Actual source code: svdksvd.c

slepc-3.20.2 2024-03-15
Report Typos and Errors
  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:   PetscCall(MatGetSize(svd->A,&M,&N));
 33:   PetscCheck(M==N,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The interface to KSVD does not support rectangular matrices");
 34:   svd->ncv = N;
 35:   if (svd->mpd!=PETSC_DEFAULT) PetscCall(PetscInfo(svd,"Warning: parameter mpd ignored\n"));
 36:   if (svd->max_it==PETSC_DEFAULT) svd->max_it = 1;
 37:   svd->leftbasis = PETSC_TRUE;
 38:   SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
 39:   PetscCall(SVDAllocateSolution(svd,0));

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

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

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

 65:   PetscFunctionBegin;
 66:   PetscCall(MatGetSize(A,&M,&N));
 67:   PetscCall(MatGetLocalSize(A,&m,&n));
 68:   mn = (M>=N)? n: m;
 69:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&Z));
 70:   PetscCall(MatSetSizes(Z,m,mn,PETSC_DECIDE,PETSC_DECIDE));
 71:   PetscCall(MatSetType(Z,MATSCALAPACK));
 72:   PetscCall(MatSetUp(Z));
 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(MatSetUp(Q));
 80:   PetscCall(MatAssemblyBegin(Q,MAT_FINAL_ASSEMBLY));
 81:   PetscCall(MatAssemblyEnd(Q,MAT_FINAL_ASSEMBLY));
 82:   q = (Mat_ScaLAPACK*)Q->data;

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

 92:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 93:   /* allocate workspace */
 94:   PetscStackCallExternalVoid("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));
 95:   PetscCheck(!info,PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"Error in KSVD subroutine pdgeqsvd: info=%d",(int)info);
 96:   PetscCall(PetscBLASIntCast((PetscInt)minlwork,&lwork));
 97:   PetscCall(PetscBLASIntCast(minliwork,&liwork));
 98:   PetscCall(PetscMalloc2(lwork,&work,liwork,&iwork));
 99:   /* call computational routine */
100:   PetscStackCallExternalVoid("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));
101:   PetscCheck(!info,PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"Error in KSVD subroutine pdgeqsvd: info=%d",(int)info);
102:   PetscCall(PetscFree2(work,iwork));
103:   PetscCall(PetscFPTrapPop());

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

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

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

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

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

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

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

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

167:   PetscFunctionBegin;
168:   PetscOptionsHeadBegin(PetscOptionsObject,"SVD KSVD Options");

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

175:   PetscOptionsHeadEnd();
176:   PetscFunctionReturn(PETSC_SUCCESS);
177: }

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

183:   PetscFunctionBegin;
184:   ctx->eigen = eigen;
185:   PetscFunctionReturn(PETSC_SUCCESS);
186: }

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

191:    Logically Collective

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

197:    Options Database Key:
198: .  -svd_ksvd_eigen_method - Sets the method for the KSVD eigensolver

200:    If not set, the method defaults to SVD_KSVD_EIGEN_MRRR.

202:    Level: advanced

204: .seealso: SVDKSVDGetEigenMethod(), SVDKSVDEigenMethod
205: @*/
206: PetscErrorCode SVDKSVDSetEigenMethod(SVD svd,SVDKSVDEigenMethod eigen)
207: {
208:   PetscFunctionBegin;
211:   PetscTryMethod(svd,"SVDKSVDSetEigenMethod_C",(SVD,SVDKSVDEigenMethod),(svd,eigen));
212:   PetscFunctionReturn(PETSC_SUCCESS);
213: }

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

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

224: /*@
225:    SVDKSVDGetEigenMethod - Gets the method for the KSVD eigensolver.

227:    Not Collective

229:    Input Parameter:
230: .  svd - the singular value solver context

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

235:    Level: advanced

237: .seealso: SVDKSVDSetEigenMethod(), SVDKSVDEigenMethod
238: @*/
239: PetscErrorCode SVDKSVDGetEigenMethod(SVD svd,SVDKSVDEigenMethod *eigen)
240: {
241:   PetscFunctionBegin;
243:   PetscAssertPointer(eigen,2);
244:   PetscUseMethod(svd,"SVDKSVDGetEigenMethod_C",(SVD,SVDKSVDEigenMethod*),(svd,eigen));
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

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

252:   PetscFunctionBegin;
253:   ctx->polar = polar;
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

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

260:    Logically Collective

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

266:    Options Database Key:
267: .  -svd_ksvd_polar_method - Sets the method for the KSVD polar decomposition

269:    If not set, the method defaults to SVD_KSVD_POLAR_QDWH.

271:    Level: advanced

273: .seealso: SVDKSVDGetPolarMethod(), SVDKSVDPolarMethod
274: @*/
275: PetscErrorCode SVDKSVDSetPolarMethod(SVD svd,SVDKSVDPolarMethod polar)
276: {
277:   PetscFunctionBegin;
280:   PetscTryMethod(svd,"SVDKSVDSetPolarMethod_C",(SVD,SVDKSVDPolarMethod),(svd,polar));
281:   PetscFunctionReturn(PETSC_SUCCESS);
282: }

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

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

293: /*@
294:    SVDKSVDGetPolarMethod - Gets the method for the KSVD polar decomposition.

296:    Not Collective

298:    Input Parameter:
299: .  svd - the singular value solver context

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

304:    Level: advanced

306: .seealso: SVDKSVDSetPolarMethod(), SVDKSVDPolarMethod
307: @*/
308: PetscErrorCode SVDKSVDGetPolarMethod(SVD svd,SVDKSVDPolarMethod *polar)
309: {
310:   PetscFunctionBegin;
312:   PetscAssertPointer(polar,2);
313:   PetscUseMethod(svd,"SVDKSVDGetPolarMethod_C",(SVD,SVDKSVDPolarMethod*),(svd,polar));
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

317: SLEPC_EXTERN PetscErrorCode SVDCreate_KSVD(SVD svd)
318: {
319:   SVD_KSVD  *ctx;

321:   PetscFunctionBegin;
322:   PetscCall(PetscNew(&ctx));
323:   svd->data = (void*)ctx;

325:   svd->ops->solve          = SVDSolve_KSVD;
326:   svd->ops->setup          = SVDSetUp_KSVD;
327:   svd->ops->setfromoptions = SVDSetFromOptions_KSVD;
328:   svd->ops->destroy        = SVDDestroy_KSVD;
329:   svd->ops->reset          = SVDReset_KSVD;
330:   svd->ops->view           = SVDView_KSVD;

332:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetEigenMethod_C",SVDKSVDSetEigenMethod_KSVD));
333:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetEigenMethod_C",SVDKSVDGetEigenMethod_KSVD));
334:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetPolarMethod_C",SVDKSVDSetPolarMethod_KSVD));
335:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetPolarMethod_C",SVDKSVDGetPolarMethod_KSVD));
336:   PetscFunctionReturn(PETSC_SUCCESS);
337: }