Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : /*
      11             :    SLEPc singular value solver: "randomized"
      12             : 
      13             :    Method: RSVD
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Randomized singular value decomposition.
      18             : 
      19             :    References:
      20             : 
      21             :        [1] N. Halko, P.-G. Martinsson, and J. A. Tropp, "Finding
      22             :            structure with randomness: Probabilistic algorithms for
      23             :            constructing approximate matrix decompositions", SIAM Rev.,
      24             :            53(2):217-288, 2011.
      25             : */
      26             : 
      27             : #include <slepc/private/svdimpl.h>                /*I "slepcsvd.h" I*/
      28             : 
      29          10 : static PetscErrorCode SVDSetUp_Randomized(SVD svd)
      30             : {
      31          10 :   PetscInt       N;
      32             : 
      33          10 :   PetscFunctionBegin;
      34          10 :   SVDCheckStandard(svd);
      35          10 :   SVDCheckDefinite(svd);
      36          10 :   PetscCheck(svd->which==SVD_LARGEST,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"This solver supports only largest singular values");
      37          10 :   PetscCall(MatGetSize(svd->A,NULL,&N));
      38          10 :   PetscCall(SVDSetDimensions_Default(svd));
      39          10 :   PetscCheck(svd->ncv>=svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must not be smaller than nsv");
      40          10 :   if (svd->max_it==PETSC_DETERMINE) svd->max_it = PetscMax(N/svd->ncv,100);
      41          10 :   svd->leftbasis = PETSC_TRUE;
      42          10 :   svd->mpd = svd->ncv;
      43          10 :   PetscCall(SVDAllocateSolution(svd,0));
      44          10 :   PetscCall(DSSetType(svd->ds,DSSVD));
      45          10 :   PetscCall(DSAllocate(svd->ds,svd->ncv));
      46          10 :   PetscCall(SVDSetWorkVecs(svd,1,1));
      47          10 :   PetscFunctionReturn(PETSC_SUCCESS);
      48             : }
      49             : 
      50         313 : static PetscErrorCode SVDRandomizedResidualNorm(SVD svd,PetscInt i,PetscScalar sigma,PetscReal *res)
      51             : {
      52         313 :   PetscReal      norm1,norm2;
      53         313 :   Vec            u,v,wu,wv;
      54             : 
      55         313 :   PetscFunctionBegin;
      56         313 :   *res = 1.0;
      57         313 :   if (svd->conv!=SVD_CONV_MAXIT) {
      58         311 :     wu = svd->swapped? svd->workr[0]: svd->workl[0];
      59         311 :     wv = svd->swapped? svd->workl[0]: svd->workr[0];
      60         311 :     PetscCall(BVGetColumn(svd->V,i,&v));
      61         311 :     PetscCall(BVGetColumn(svd->U,i,&u));
      62             :     /* norm1 = ||A*v-sigma*u||_2 */
      63         311 :     PetscCall(MatMult(svd->A,v,wu));
      64         311 :     PetscCall(VecAXPY(wu,-sigma,u));
      65         311 :     PetscCall(VecNorm(wu,NORM_2,&norm1));
      66             :     /* norm2 = ||A^T*u-sigma*v||_2 */
      67         311 :     PetscCall(MatMult(svd->AT,u,wv));
      68         311 :     PetscCall(VecAXPY(wv,-sigma,v));
      69         311 :     PetscCall(VecNorm(wv,NORM_2,&norm2));
      70         311 :     PetscCall(BVRestoreColumn(svd->V,i,&v));
      71         311 :     PetscCall(BVRestoreColumn(svd->U,i,&u));
      72         311 :     *res = SlepcAbs(norm1,norm2);
      73             :   }
      74         313 :   PetscFunctionReturn(PETSC_SUCCESS);
      75             : }
      76             : 
      77             : /* If A is a virtual Hermitian transpose, then BVMatMult will fail if PRODUCT_AhB is not implemented */
      78         544 : static PetscErrorCode BlockMatMult(BV V,Mat A,BV Y,Mat AT)
      79             : {
      80         544 :   PetscFunctionBegin;
      81         544 :   if (!PetscDefined(USE_COMPLEX)) PetscCall(BVMatMult(V,A,Y));
      82             :   else {
      83         544 :     PetscBool flg=PETSC_FALSE;
      84         544 :     PetscCall(PetscObjectTypeCompare((PetscObject)A,MATHERMITIANTRANSPOSEVIRTUAL,&flg));
      85         544 :     if (flg) PetscCall(BVMatMultHermitianTranspose(V,AT,Y));
      86         523 :     else PetscCall(BVMatMult(V,A,Y));
      87             :   }
      88         544 :   PetscFunctionReturn(PETSC_SUCCESS);
      89             : }
      90             : 
      91          10 : static PetscErrorCode SVDSolve_Randomized(SVD svd)
      92             : {
      93          10 :   PetscScalar    *w;
      94          10 :   PetscReal      res=1.0;
      95          10 :   PetscInt       i,k=0;
      96          10 :   Mat            A,U,V;
      97             : 
      98          10 :   PetscFunctionBegin;
      99             :   /* Form random matrix, G. Complete the initial basis with random vectors */
     100          10 :   PetscCall(BVSetActiveColumns(svd->V,svd->nini,svd->ncv));
     101          10 :   PetscCall(BVSetRandomNormal(svd->V));
     102          10 :   PetscCall(PetscCalloc1(svd->ncv,&w));
     103             : 
     104             :   /* Subspace Iteration */
     105         272 :   do {
     106         272 :     svd->its++;
     107         272 :     PetscCall(BVSetActiveColumns(svd->V,svd->nconv,svd->ncv));
     108         272 :     PetscCall(BVSetActiveColumns(svd->U,svd->nconv,svd->ncv));
     109             :     /* Form AG */
     110         272 :     PetscCall(BlockMatMult(svd->V,svd->A,svd->U,svd->AT));
     111             :     /* Orthogonalization Q=qr(AG)*/
     112         272 :     PetscCall(BVOrthogonalize(svd->U,NULL));
     113             :     /* Form B^*= AQ */
     114         272 :     PetscCall(BlockMatMult(svd->U,svd->AT,svd->V,svd->A));
     115             : 
     116         272 :     PetscCall(DSSetDimensions(svd->ds,svd->ncv,svd->nconv,svd->ncv));
     117         272 :     PetscCall(DSSVDSetDimensions(svd->ds,svd->ncv));
     118         272 :     PetscCall(DSGetMat(svd->ds,DS_MAT_A,&A));
     119         272 :     PetscCall(MatZeroEntries(A));
     120         272 :     PetscCall(BVOrthogonalize(svd->V,A));
     121         272 :     PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&A));
     122         272 :     PetscCall(DSSetState(svd->ds,DS_STATE_RAW));
     123         272 :     PetscCall(DSSolve(svd->ds,w,NULL));
     124         272 :     PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
     125         272 :     PetscCall(DSSynchronize(svd->ds,w,NULL));
     126         272 :     PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
     127         272 :     PetscCall(DSGetMat(svd->ds,DS_MAT_V,&V));
     128         272 :     PetscCall(BVMultInPlace(svd->U,V,svd->nconv,svd->ncv));
     129         272 :     PetscCall(BVMultInPlace(svd->V,U,svd->nconv,svd->ncv));
     130         272 :     PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
     131         272 :     PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&V));
     132             :     /* Check convergence */
     133         272 :     k = 0;
     134         314 :     for (i=svd->nconv;i<svd->ncv;i++) {
     135         313 :       PetscCall(SVDRandomizedResidualNorm(svd,i,w[i],&res));
     136         313 :       svd->sigma[i] = PetscRealPart(w[i]);
     137         313 :       PetscCall((*svd->converged)(svd,svd->sigma[i],res,&svd->errest[i],svd->convergedctx));
     138         313 :       if (svd->errest[i] < svd->tol) k++;
     139             :       else break;
     140             :     }
     141         272 :     if (svd->conv == SVD_CONV_MAXIT && svd->its >= svd->max_it) {
     142           2 :       k = svd->nsv;
     143          22 :       for (i=0;i<svd->ncv;i++) svd->sigma[i] = PetscRealPart(w[i]);
     144             :     }
     145         272 :     PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,svd->nconv+k,svd->nsv,&svd->reason,svd->stoppingctx));
     146         272 :     svd->nconv += k;
     147         272 :     PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,svd->ncv));
     148         272 :   } while (svd->reason == SVD_CONVERGED_ITERATING);
     149          10 :   PetscCall(PetscFree(w));
     150          10 :   PetscFunctionReturn(PETSC_SUCCESS);
     151             : }
     152             : 
     153           8 : SLEPC_EXTERN PetscErrorCode SVDCreate_Randomized(SVD svd)
     154             : {
     155           8 :   PetscFunctionBegin;
     156           8 :   svd->ops->setup          = SVDSetUp_Randomized;
     157           8 :   svd->ops->solve          = SVDSolve_Randomized;
     158           8 :   PetscFunctionReturn(PETSC_SUCCESS);
     159             : }