Actual source code: svdksvd.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */
     10: /*
     11:    This file implements a wrapper to the KSVD SVD solver
     12: */

     14: #include <slepc/private/svdimpl.h>
     15: #include <slepc/private/slepcscalapack.h>
     16: #include <ksvd.h>

     18: typedef struct {
     19:   Mat                As;        /* converted matrix */
     20:   SVDKSVDEigenMethod eigen;
     21:   SVDKSVDPolarMethod polar;
     22: } SVD_KSVD;

     24: static PetscErrorCode SVDSetUp_KSVD(SVD svd)
     25: {
     26:   SVD_KSVD  *ctx = (SVD_KSVD*)svd->data;
     27:   PetscInt  M,N;

     29:   PetscFunctionBegin;
     30:   SVDCheckStandard(svd);
     31:   SVDCheckDefinite(svd);
     32:   if (svd->nsv==0) svd->nsv = 1;
     33:   PetscCall(MatGetSize(svd->A,&M,&N));
     34:   PetscCheck(M==N,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The interface to KSVD does not support rectangular matrices");
     35:   svd->ncv = N;
     36:   if (svd->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(svd,"Warning: parameter mpd ignored\n"));
     37:   if (svd->max_it==PETSC_DETERMINE) svd->max_it = 1;
     38:   svd->leftbasis = PETSC_TRUE;
     39:   SVDCheckIgnored(svd,SVD_FEATURE_STOPPING);
     40:   PetscCall(SVDAllocateSolution(svd,0));

     42:   /* default methods */
     43:   if (!ctx->eigen) ctx->eigen = SVD_KSVD_EIGEN_MRRR;
     44:   if (!ctx->polar) ctx->polar = SVD_KSVD_POLAR_QDWH;

     46:   /* convert matrix */
     47:   PetscCall(MatDestroy(&ctx->As));
     48:   PetscCall(MatConvert(svd->OP,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->As));
     49:   PetscFunctionReturn(PETSC_SUCCESS);
     50: }

     52: static PetscErrorCode SVDSolve_KSVD(SVD svd)
     53: {
     54:   SVD_KSVD       *ctx = (SVD_KSVD*)svd->data;
     55:   Mat            A = ctx->As,Z,Q,U,V;
     56:   Mat_ScaLAPACK  *a = (Mat_ScaLAPACK*)A->data,*q,*z;
     57:   PetscScalar    *work,minlwork;
     58:   PetscBLASInt   info,lwork=-1,*iwork,liwork=-1,minliwork,one=1;
     59:   PetscInt       M,N,m,n,mn;
     60:   const char     *eigen;
     61: #if defined(PETSC_USE_COMPLEX)
     62:   PetscBLASInt   lrwork;
     63:   PetscReal      *rwork,dummy;
     64: #endif

     66:   PetscFunctionBegin;
     67:   PetscCall(MatGetSize(A,&M,&N));
     68:   PetscCall(MatGetLocalSize(A,&m,&n));
     69:   mn = (M>=N)? n: m;
     70:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&Z));
     71:   PetscCall(MatSetSizes(Z,m,mn,PETSC_DECIDE,PETSC_DECIDE));
     72:   PetscCall(MatSetType(Z,MATSCALAPACK));
     73:   PetscCall(MatAssemblyBegin(Z,MAT_FINAL_ASSEMBLY));
     74:   PetscCall(MatAssemblyEnd(Z,MAT_FINAL_ASSEMBLY));
     75:   z = (Mat_ScaLAPACK*)Z->data;
     76:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&Q));
     77:   PetscCall(MatSetSizes(Q,mn,n,PETSC_DECIDE,PETSC_DECIDE));
     78:   PetscCall(MatSetType(Q,MATSCALAPACK));
     79:   PetscCall(MatAssemblyBegin(Q,MAT_FINAL_ASSEMBLY));
     80:   PetscCall(MatAssemblyEnd(Q,MAT_FINAL_ASSEMBLY));
     81:   q = (Mat_ScaLAPACK*)Q->data;

     83:   /* configure solver */
     84:   setPolarAlgorithm((int)ctx->polar);
     85:   switch (ctx->eigen) {
     86:     case SVD_KSVD_EIGEN_MRRR: eigen = "r"; break;
     87:     case SVD_KSVD_EIGEN_DC:   eigen = "d"; break;
     88:     case SVD_KSVD_EIGEN_ELPA: eigen = "e"; break;
     89:   }

     91:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     92:   /* allocate workspace */
     93:   PetscCallExternalVoid("pdgeqsvd",pdgeqsvd("V","V",eigen,a->M,a->N,a->loc,one,one,a->desc,svd->sigma,z->loc,one,one,z->desc,q->loc,one,one,q->desc,&minlwork,lwork,&minliwork,liwork,&info));
     94:   PetscCheck(!info,PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"Error in KSVD subroutine pdgeqsvd: info=%d",(int)info);
     95:   PetscCall(PetscBLASIntCast((PetscInt)minlwork,&lwork));
     96:   PetscCall(PetscBLASIntCast(minliwork,&liwork));
     97:   PetscCall(PetscMalloc2(lwork,&work,liwork,&iwork));
     98:   /* call computational routine */
     99:   PetscCallExternalVoid("pdgeqsvd",pdgeqsvd("V","V",eigen,a->M,a->N,a->loc,one,one,a->desc,svd->sigma,z->loc,one,one,z->desc,q->loc,one,one,q->desc,work,lwork,iwork,liwork,&info));
    100:   PetscCheck(!info,PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"Error in KSVD subroutine pdgeqsvd: info=%d",(int)info);
    101:   PetscCall(PetscFree2(work,iwork));
    102:   PetscCall(PetscFPTrapPop());

    104:   PetscCall(BVGetMat(svd->U,&U));
    105:   PetscCall(BVGetMat(svd->V,&V));
    106:   if (M>=N) {
    107:     PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&U));
    108:     PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
    109:   } else {
    110:     PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&U));
    111:     PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&V));
    112:   }
    113:   PetscCall(BVRestoreMat(svd->U,&U));
    114:   PetscCall(BVRestoreMat(svd->V,&V));
    115:   PetscCall(MatDestroy(&Z));
    116:   PetscCall(MatDestroy(&Q));

    118:   svd->nconv  = svd->ncv;
    119:   svd->its    = 1;
    120:   svd->reason = SVD_CONVERGED_TOL;
    121:   PetscFunctionReturn(PETSC_SUCCESS);
    122: }

    124: static PetscErrorCode SVDDestroy_KSVD(SVD svd)
    125: {
    126:   PetscFunctionBegin;
    127:   PetscCall(PetscFree(svd->data));
    128:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetEigenMethod_C",NULL));
    129:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetEigenMethod_C",NULL));
    130:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetPolarMethod_C",NULL));
    131:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetPolarMethod_C",NULL));
    132:   PetscFunctionReturn(PETSC_SUCCESS);
    133: }

    135: static PetscErrorCode SVDReset_KSVD(SVD svd)
    136: {
    137:   SVD_KSVD  *ctx = (SVD_KSVD*)svd->data;

    139:   PetscFunctionBegin;
    140:   PetscCall(MatDestroy(&ctx->As));
    141:   PetscFunctionReturn(PETSC_SUCCESS);
    142: }

    144: static PetscErrorCode SVDView_KSVD(SVD svd,PetscViewer viewer)
    145: {
    146:   PetscBool      isascii;
    147:   SVD_KSVD       *ctx = (SVD_KSVD*)svd->data;

    149:   PetscFunctionBegin;
    150:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    151:   if (isascii) {
    152:     PetscCall(PetscViewerASCIIPrintf(viewer,"  eigensolver method: %s\n",SVDKSVDEigenMethods[ctx->eigen]));
    153:     PetscCall(PetscViewerASCIIPrintf(viewer,"  polar decomposition method: %s\n",SVDKSVDPolarMethods[ctx->polar]));
    154:   }
    155:   PetscFunctionReturn(PETSC_SUCCESS);
    156: }

    158: static PetscErrorCode SVDSetFromOptions_KSVD(SVD svd,PetscOptionItems PetscOptionsObject)
    159: {
    160:   SVD_KSVD           *ctx = (SVD_KSVD*)svd->data;
    161:   SVDKSVDEigenMethod eigen;
    162:   SVDKSVDPolarMethod polar;
    163:   PetscBool          flg;

    165:   PetscFunctionBegin;
    166:   PetscOptionsHeadBegin(PetscOptionsObject,"SVD KSVD Options");

    168:     PetscCall(PetscOptionsEnum("-svd_ksvd_eigen_method","Method for solving the internal eigenvalue problem","SVDKSVDSetEigenMethod",SVDKSVDEigenMethods,(PetscEnum)ctx->eigen,(PetscEnum*)&eigen,&flg));
    169:     if (flg) PetscCall(SVDKSVDSetEigenMethod(svd,eigen));
    170:     PetscCall(PetscOptionsEnum("-svd_ksvd_polar_method","Method for computing the polar decomposition","SVDKSVDSetPolarMethod",SVDKSVDPolarMethods,(PetscEnum)ctx->polar,(PetscEnum*)&polar,&flg));
    171:     if (flg) PetscCall(SVDKSVDSetPolarMethod(svd,polar));

    173:   PetscOptionsHeadEnd();
    174:   PetscFunctionReturn(PETSC_SUCCESS);
    175: }

    177: static PetscErrorCode SVDKSVDSetEigenMethod_KSVD(SVD svd,SVDKSVDEigenMethod eigen)
    178: {
    179:   SVD_KSVD *ctx = (SVD_KSVD*)svd->data;

    181:   PetscFunctionBegin;
    182:   ctx->eigen = eigen;
    183:   PetscFunctionReturn(PETSC_SUCCESS);
    184: }

    186: /*@
    187:    SVDKSVDSetEigenMethod - Sets the method to solve the internal eigenproblem within the KSVD solver.

    189:    Logically Collective

    191:    Input Parameters:
    192: +  svd - the singular value solver context
    193: -  eigen - method that will be used by KSVD for the eigenproblem

    195:    Options Database Key:
    196: .  -svd_ksvd_eigen_method (mrrr|dc|elpa) - sets the method for the KSVD eigensolver

    198:    If not set, the method defaults to `SVD_KSVD_EIGEN_MRRR`.

    200:    Level: advanced

    202: .seealso: [](ch:svd), `SVDKSVD`, `SVDKSVDGetEigenMethod()`, `SVDKSVDEigenMethod`
    203: @*/
    204: PetscErrorCode SVDKSVDSetEigenMethod(SVD svd,SVDKSVDEigenMethod eigen)
    205: {
    206:   PetscFunctionBegin;
    209:   PetscTryMethod(svd,"SVDKSVDSetEigenMethod_C",(SVD,SVDKSVDEigenMethod),(svd,eigen));
    210:   PetscFunctionReturn(PETSC_SUCCESS);
    211: }

    213: static PetscErrorCode SVDKSVDGetEigenMethod_KSVD(SVD svd,SVDKSVDEigenMethod *eigen)
    214: {
    215:   SVD_KSVD *ctx = (SVD_KSVD*)svd->data;

    217:   PetscFunctionBegin;
    218:   *eigen = ctx->eigen;
    219:   PetscFunctionReturn(PETSC_SUCCESS);
    220: }

    222: /*@
    223:    SVDKSVDGetEigenMethod - Gets the method for the KSVD eigensolver.

    225:    Not Collective

    227:    Input Parameter:
    228: .  svd - the singular value solver context

    230:    Output Parameter:
    231: .  eigen - method that will be used by KSVD for the eigenproblem

    233:    Level: advanced

    235: .seealso: [](ch:svd), `SVDKSVD`, `SVDKSVDSetEigenMethod()`, `SVDKSVDEigenMethod`
    236: @*/
    237: PetscErrorCode SVDKSVDGetEigenMethod(SVD svd,SVDKSVDEigenMethod *eigen)
    238: {
    239:   PetscFunctionBegin;
    241:   PetscAssertPointer(eigen,2);
    242:   PetscUseMethod(svd,"SVDKSVDGetEigenMethod_C",(SVD,SVDKSVDEigenMethod*),(svd,eigen));
    243:   PetscFunctionReturn(PETSC_SUCCESS);
    244: }

    246: static PetscErrorCode SVDKSVDSetPolarMethod_KSVD(SVD svd,SVDKSVDPolarMethod polar)
    247: {
    248:   SVD_KSVD *ctx = (SVD_KSVD*)svd->data;

    250:   PetscFunctionBegin;
    251:   ctx->polar = polar;
    252:   PetscFunctionReturn(PETSC_SUCCESS);
    253: }

    255: /*@
    256:    SVDKSVDSetPolarMethod - Sets the method to compute the polar decomposition within the KSVD solver.

    258:    Logically Collective

    260:    Input Parameters:
    261: +  svd - the singular value solver context
    262: -  polar - method that will be used by KSVD for the polar decomposition

    264:    Options Database Key:
    265: .  -svd_ksvd_polar_method (qdwh|zolopd) - sets the method for the KSVD polar decomposition

    267:    If not set, the method defaults to `SVD_KSVD_POLAR_QDWH`.

    269:    Level: advanced

    271: .seealso: [](ch:svd), `SVDKSVD`, `SVDKSVDGetPolarMethod()`, `SVDKSVDPolarMethod`
    272: @*/
    273: PetscErrorCode SVDKSVDSetPolarMethod(SVD svd,SVDKSVDPolarMethod polar)
    274: {
    275:   PetscFunctionBegin;
    278:   PetscTryMethod(svd,"SVDKSVDSetPolarMethod_C",(SVD,SVDKSVDPolarMethod),(svd,polar));
    279:   PetscFunctionReturn(PETSC_SUCCESS);
    280: }

    282: static PetscErrorCode SVDKSVDGetPolarMethod_KSVD(SVD svd,SVDKSVDPolarMethod *polar)
    283: {
    284:   SVD_KSVD *ctx = (SVD_KSVD*)svd->data;

    286:   PetscFunctionBegin;
    287:   *polar = ctx->polar;
    288:   PetscFunctionReturn(PETSC_SUCCESS);
    289: }

    291: /*@
    292:    SVDKSVDGetPolarMethod - Gets the method for the KSVD polar decomposition.

    294:    Not Collective

    296:    Input Parameter:
    297: .  svd - the singular value solver context

    299:    Output Parameter:
    300: .  polar - method that will be used by KSVD for the polar decomposition

    302:    Level: advanced

    304: .seealso: [](ch:svd), `SVDKSVD`, `SVDKSVDSetPolarMethod()`, `SVDKSVDPolarMethod`
    305: @*/
    306: PetscErrorCode SVDKSVDGetPolarMethod(SVD svd,SVDKSVDPolarMethod *polar)
    307: {
    308:   PetscFunctionBegin;
    310:   PetscAssertPointer(polar,2);
    311:   PetscUseMethod(svd,"SVDKSVDGetPolarMethod_C",(SVD,SVDKSVDPolarMethod*),(svd,polar));
    312:   PetscFunctionReturn(PETSC_SUCCESS);
    313: }

    315: /*MC
    316:    SVDKSVD - SVDKSVD = "ksvd" - A wrapper to the KSVD singular value
    317:    solver {cite:p}`Suk19`.

    319:    Notes:
    320:    Only available for standard SVD problems.

    322:    This is a direct singular value solver, that is, the full decomposition
    323:    is computed. The computation involves redistributing the matrices from PETSc
    324:    storage to ScaLAPACK distribution, and vice versa (this is done automatically
    325:    by SLEPc). Alternatively, the user may create the problem matrices
    326:    already with type `MATSCALAPACK`.

    328:    The implemented method is supposed to be more scalable than the one
    329:    in `SVDSCALAPACK`.

    331:    Level: beginner

    333: .seealso: [](ch:svd), `SVD`, `SVDType`, `SVDSetType()`, `SVDSCALAPACK`
    334: M*/
    335: SLEPC_EXTERN PetscErrorCode SVDCreate_KSVD(SVD svd)
    336: {
    337:   SVD_KSVD  *ctx;

    339:   PetscFunctionBegin;
    340:   PetscCall(PetscNew(&ctx));
    341:   svd->data = (void*)ctx;

    343:   svd->ops->solve          = SVDSolve_KSVD;
    344:   svd->ops->setup          = SVDSetUp_KSVD;
    345:   svd->ops->setfromoptions = SVDSetFromOptions_KSVD;
    346:   svd->ops->destroy        = SVDDestroy_KSVD;
    347:   svd->ops->reset          = SVDReset_KSVD;
    348:   svd->ops->view           = SVDView_KSVD;

    350:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetEigenMethod_C",SVDKSVDSetEigenMethod_KSVD));
    351:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetEigenMethod_C",SVDKSVDGetEigenMethod_KSVD));
    352:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetPolarMethod_C",SVDKSVDSetPolarMethod_KSVD));
    353:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetPolarMethod_C",SVDKSVDGetPolarMethod_KSVD));
    354:   PetscFunctionReturn(PETSC_SUCCESS);
    355: }