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: 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));
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: 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));
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 \<eigen\> - 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 \<polar\> - 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: }