Actual source code: svdksvd.c
slepc-3.20.2 2024-03-15
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: }