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 :   SVDCheckUnsupported(svd,SVD_FEATURE_THRESHOLD);
      37          10 :   PetscCheck(svd->which==SVD_LARGEST,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"This solver supports only largest singular values");
      38          10 :   PetscCall(MatGetSize(svd->A,NULL,&N));
      39          10 :   PetscCall(SVDSetDimensions_Default(svd));
      40          10 :   PetscCheck(svd->ncv>=svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must not be smaller than nsv");
      41          10 :   if (svd->max_it==PETSC_DETERMINE) svd->max_it = PetscMax(N/svd->ncv,100);
      42          10 :   svd->leftbasis = PETSC_TRUE;
      43          10 :   svd->mpd = svd->ncv;
      44          10 :   PetscCall(SVDAllocateSolution(svd,0));
      45          10 :   PetscCall(DSSetType(svd->ds,DSSVD));
      46          10 :   PetscCall(DSAllocate(svd->ds,svd->ncv));
      47          10 :   PetscCall(SVDSetWorkVecs(svd,1,1));
      48          10 :   PetscFunctionReturn(PETSC_SUCCESS);
      49             : }
      50             : 
      51         313 : static PetscErrorCode SVDRandomizedResidualNorm(SVD svd,PetscInt i,PetscScalar sigma,PetscReal *res)
      52             : {
      53         313 :   PetscReal      norm1,norm2;
      54         313 :   Vec            u,v,wu,wv;
      55             : 
      56         313 :   PetscFunctionBegin;
      57         313 :   *res = 1.0;
      58         313 :   if (svd->conv!=SVD_CONV_MAXIT) {
      59         311 :     wu = svd->swapped? svd->workr[0]: svd->workl[0];
      60         311 :     wv = svd->swapped? svd->workl[0]: svd->workr[0];
      61         311 :     PetscCall(BVGetColumn(svd->V,i,&v));
      62         311 :     PetscCall(BVGetColumn(svd->U,i,&u));
      63             :     /* norm1 = ||A*v-sigma*u||_2 */
      64         311 :     PetscCall(MatMult(svd->A,v,wu));
      65         311 :     PetscCall(VecAXPY(wu,-sigma,u));
      66         311 :     PetscCall(VecNorm(wu,NORM_2,&norm1));
      67             :     /* norm2 = ||A^T*u-sigma*v||_2 */
      68         311 :     PetscCall(MatMult(svd->AT,u,wv));
      69         311 :     PetscCall(VecAXPY(wv,-sigma,v));
      70         311 :     PetscCall(VecNorm(wv,NORM_2,&norm2));
      71         311 :     PetscCall(BVRestoreColumn(svd->V,i,&v));
      72         311 :     PetscCall(BVRestoreColumn(svd->U,i,&u));
      73         311 :     *res = SlepcAbs(norm1,norm2);
      74             :   }
      75         313 :   PetscFunctionReturn(PETSC_SUCCESS);
      76             : }
      77             : 
      78             : /* If A is a virtual Hermitian transpose, then BVMatMult will fail if PRODUCT_AhB is not implemented */
      79         544 : static PetscErrorCode BlockMatMult(BV V,Mat A,BV Y,Mat AT)
      80             : {
      81         544 :   PetscFunctionBegin;
      82         544 :   if (!PetscDefined(USE_COMPLEX)) PetscCall(BVMatMult(V,A,Y));
      83             :   else {
      84         544 :     PetscBool flg=PETSC_FALSE;
      85         544 :     PetscCall(PetscObjectTypeCompare((PetscObject)A,MATHERMITIANTRANSPOSEVIRTUAL,&flg));
      86         544 :     if (flg) PetscCall(BVMatMultHermitianTranspose(V,AT,Y));
      87         523 :     else PetscCall(BVMatMult(V,A,Y));
      88             :   }
      89         544 :   PetscFunctionReturn(PETSC_SUCCESS);
      90             : }
      91             : 
      92          10 : static PetscErrorCode SVDSolve_Randomized(SVD svd)
      93             : {
      94          10 :   PetscScalar    *w;
      95          10 :   PetscReal      res=1.0;
      96          10 :   PetscInt       i,k=0;
      97          10 :   Mat            A,U,V;
      98             : 
      99          10 :   PetscFunctionBegin;
     100             :   /* Form random matrix, G. Complete the initial basis with random vectors */
     101          10 :   PetscCall(BVSetActiveColumns(svd->V,svd->nini,svd->ncv));
     102          10 :   PetscCall(BVSetRandomNormal(svd->V));
     103          10 :   PetscCall(PetscCalloc1(svd->ncv,&w));
     104             : 
     105             :   /* Subspace Iteration */
     106         272 :   do {
     107         272 :     svd->its++;
     108         272 :     PetscCall(BVSetActiveColumns(svd->V,svd->nconv,svd->ncv));
     109         272 :     PetscCall(BVSetActiveColumns(svd->U,svd->nconv,svd->ncv));
     110             :     /* Form AG */
     111         272 :     PetscCall(BlockMatMult(svd->V,svd->A,svd->U,svd->AT));
     112             :     /* Orthogonalization Q=qr(AG)*/
     113         272 :     PetscCall(BVOrthogonalize(svd->U,NULL));
     114             :     /* Form B^*= AQ */
     115         272 :     PetscCall(BlockMatMult(svd->U,svd->AT,svd->V,svd->A));
     116             : 
     117         272 :     PetscCall(DSSetDimensions(svd->ds,svd->ncv,svd->nconv,svd->ncv));
     118         272 :     PetscCall(DSSVDSetDimensions(svd->ds,svd->ncv));
     119         272 :     PetscCall(DSGetMat(svd->ds,DS_MAT_A,&A));
     120         272 :     PetscCall(MatZeroEntries(A));
     121         272 :     PetscCall(BVOrthogonalize(svd->V,A));
     122         272 :     PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&A));
     123         272 :     PetscCall(DSSetState(svd->ds,DS_STATE_RAW));
     124         272 :     PetscCall(DSSolve(svd->ds,w,NULL));
     125         272 :     PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
     126         272 :     PetscCall(DSSynchronize(svd->ds,w,NULL));
     127         272 :     PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
     128         272 :     PetscCall(DSGetMat(svd->ds,DS_MAT_V,&V));
     129         272 :     PetscCall(BVMultInPlace(svd->U,V,svd->nconv,svd->ncv));
     130         272 :     PetscCall(BVMultInPlace(svd->V,U,svd->nconv,svd->ncv));
     131         272 :     PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
     132         272 :     PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&V));
     133             :     /* Check convergence */
     134         272 :     k = 0;
     135         314 :     for (i=svd->nconv;i<svd->ncv;i++) {
     136         313 :       PetscCall(SVDRandomizedResidualNorm(svd,i,w[i],&res));
     137         313 :       svd->sigma[i] = PetscRealPart(w[i]);
     138         313 :       PetscCall((*svd->converged)(svd,svd->sigma[i],res,&svd->errest[i],svd->convergedctx));
     139         313 :       if (svd->errest[i] < svd->tol) k++;
     140             :       else break;
     141             :     }
     142         272 :     if (svd->conv == SVD_CONV_MAXIT && svd->its >= svd->max_it) {
     143           2 :       k = svd->nsv;
     144          22 :       for (i=0;i<svd->ncv;i++) svd->sigma[i] = PetscRealPart(w[i]);
     145             :     }
     146         272 :     PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,svd->nconv+k,svd->nsv,&svd->reason,svd->stoppingctx));
     147         272 :     svd->nconv += k;
     148         272 :     PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,svd->ncv));
     149         272 :   } while (svd->reason == SVD_CONVERGED_ITERATING);
     150          10 :   PetscCall(PetscFree(w));
     151          10 :   PetscFunctionReturn(PETSC_SUCCESS);
     152             : }
     153             : 
     154           8 : SLEPC_EXTERN PetscErrorCode SVDCreate_Randomized(SVD svd)
     155             : {
     156           8 :   PetscFunctionBegin;
     157           8 :   svd->ops->setup          = SVDSetUp_Randomized;
     158           8 :   svd->ops->solve          = SVDSolve_Randomized;
     159           8 :   PetscFunctionReturn(PETSC_SUCCESS);
     160             : }