Actual source code: svdksvd.c
slepc-main 2024-11-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_DETERMINE) PetscCall(PetscInfo(svd,"Warning: parameter mpd ignored\n"));
36: if (svd->max_it==PETSC_DETERMINE) 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: const 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(MatAssemblyBegin(Z,MAT_FINAL_ASSEMBLY));
73: PetscCall(MatAssemblyEnd(Z,MAT_FINAL_ASSEMBLY));
74: z = (Mat_ScaLAPACK*)Z->data;
75: PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&Q));
76: PetscCall(MatSetSizes(Q,mn,n,PETSC_DECIDE,PETSC_DECIDE));
77: PetscCall(MatSetType(Q,MATSCALAPACK));
78: PetscCall(MatAssemblyBegin(Q,MAT_FINAL_ASSEMBLY));
79: PetscCall(MatAssemblyEnd(Q,MAT_FINAL_ASSEMBLY));
80: q = (Mat_ScaLAPACK*)Q->data;
82: /* configure solver */
83: setPolarAlgorithm((int)ctx->polar);
84: switch (ctx->eigen) {
85: case SVD_KSVD_EIGEN_MRRR: eigen = "r"; break;
86: case SVD_KSVD_EIGEN_DC: eigen = "d"; break;
87: case SVD_KSVD_EIGEN_ELPA: eigen = "e"; break;
88: }
90: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
91: /* allocate workspace */
92: 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));
93: PetscCheck(!info,PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"Error in KSVD subroutine pdgeqsvd: info=%d",(int)info);
94: PetscCall(PetscBLASIntCast((PetscInt)minlwork,&lwork));
95: PetscCall(PetscBLASIntCast(minliwork,&liwork));
96: PetscCall(PetscMalloc2(lwork,&work,liwork,&iwork));
97: /* call computational routine */
98: 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));
99: PetscCheck(!info,PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"Error in KSVD subroutine pdgeqsvd: info=%d",(int)info);
100: PetscCall(PetscFree2(work,iwork));
101: PetscCall(PetscFPTrapPop());
103: PetscCall(BVGetMat(svd->U,&U));
104: PetscCall(BVGetMat(svd->V,&V));
105: if (M>=N) {
106: PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&U));
107: PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
108: } else {
109: PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&U));
110: PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&V));
111: }
112: PetscCall(BVRestoreMat(svd->U,&U));
113: PetscCall(BVRestoreMat(svd->V,&V));
114: PetscCall(MatDestroy(&Z));
115: PetscCall(MatDestroy(&Q));
117: svd->nconv = svd->ncv;
118: svd->its = 1;
119: svd->reason = SVD_CONVERGED_TOL;
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: static PetscErrorCode SVDDestroy_KSVD(SVD svd)
124: {
125: PetscFunctionBegin;
126: PetscCall(PetscFree(svd->data));
127: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetEigenMethod_C",NULL));
128: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetEigenMethod_C",NULL));
129: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetPolarMethod_C",NULL));
130: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetPolarMethod_C",NULL));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: static PetscErrorCode SVDReset_KSVD(SVD svd)
135: {
136: SVD_KSVD *ctx = (SVD_KSVD*)svd->data;
138: PetscFunctionBegin;
139: PetscCall(MatDestroy(&ctx->As));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: static PetscErrorCode SVDView_KSVD(SVD svd,PetscViewer viewer)
144: {
145: PetscBool isascii;
146: SVD_KSVD *ctx = (SVD_KSVD*)svd->data;
148: PetscFunctionBegin;
149: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
150: if (isascii) {
151: PetscCall(PetscViewerASCIIPrintf(viewer," eigensolver method: %s\n",SVDKSVDEigenMethods[ctx->eigen]));
152: PetscCall(PetscViewerASCIIPrintf(viewer," polar decomposition method: %s\n",SVDKSVDPolarMethods[ctx->polar]));
153: }
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode SVDSetFromOptions_KSVD(SVD svd,PetscOptionItems *PetscOptionsObject)
158: {
159: SVD_KSVD *ctx = (SVD_KSVD*)svd->data;
160: SVDKSVDEigenMethod eigen;
161: SVDKSVDPolarMethod polar;
162: PetscBool flg;
164: PetscFunctionBegin;
165: PetscOptionsHeadBegin(PetscOptionsObject,"SVD KSVD Options");
167: PetscCall(PetscOptionsEnum("-svd_ksvd_eigen_method","Method for solving the internal eigenvalue problem","SVDKSVDSetEigenMethod",SVDKSVDEigenMethods,(PetscEnum)ctx->eigen,(PetscEnum*)&eigen,&flg));
168: if (flg) PetscCall(SVDKSVDSetEigenMethod(svd,eigen));
169: PetscCall(PetscOptionsEnum("-svd_ksvd_polar_method","Method for computing the polar decomposition","SVDKSVDSetPolarMethod",SVDKSVDPolarMethods,(PetscEnum)ctx->polar,(PetscEnum*)&polar,&flg));
170: if (flg) PetscCall(SVDKSVDSetPolarMethod(svd,polar));
172: PetscOptionsHeadEnd();
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: static PetscErrorCode SVDKSVDSetEigenMethod_KSVD(SVD svd,SVDKSVDEigenMethod eigen)
177: {
178: SVD_KSVD *ctx = (SVD_KSVD*)svd->data;
180: PetscFunctionBegin;
181: ctx->eigen = eigen;
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: /*@
186: SVDKSVDSetEigenMethod - Sets the method to solve the internal eigenproblem within the KSVD solver.
188: Logically Collective
190: Input Parameters:
191: + svd - the singular value solver context
192: - eigen - method that will be used by KSVD for the eigenproblem
194: Options Database Key:
195: . -svd_ksvd_eigen_method - Sets the method for the KSVD eigensolver
197: If not set, the method defaults to SVD_KSVD_EIGEN_MRRR.
199: Level: advanced
201: .seealso: SVDKSVDGetEigenMethod(), SVDKSVDEigenMethod
202: @*/
203: PetscErrorCode SVDKSVDSetEigenMethod(SVD svd,SVDKSVDEigenMethod eigen)
204: {
205: PetscFunctionBegin;
208: PetscTryMethod(svd,"SVDKSVDSetEigenMethod_C",(SVD,SVDKSVDEigenMethod),(svd,eigen));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: static PetscErrorCode SVDKSVDGetEigenMethod_KSVD(SVD svd,SVDKSVDEigenMethod *eigen)
213: {
214: SVD_KSVD *ctx = (SVD_KSVD*)svd->data;
216: PetscFunctionBegin;
217: *eigen = ctx->eigen;
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: /*@
222: SVDKSVDGetEigenMethod - Gets the method for the KSVD eigensolver.
224: Not Collective
226: Input Parameter:
227: . svd - the singular value solver context
229: Output Parameter:
230: . eigen - method that will be used by KSVD for the eigenproblem
232: Level: advanced
234: .seealso: SVDKSVDSetEigenMethod(), SVDKSVDEigenMethod
235: @*/
236: PetscErrorCode SVDKSVDGetEigenMethod(SVD svd,SVDKSVDEigenMethod *eigen)
237: {
238: PetscFunctionBegin;
240: PetscAssertPointer(eigen,2);
241: PetscUseMethod(svd,"SVDKSVDGetEigenMethod_C",(SVD,SVDKSVDEigenMethod*),(svd,eigen));
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: static PetscErrorCode SVDKSVDSetPolarMethod_KSVD(SVD svd,SVDKSVDPolarMethod polar)
246: {
247: SVD_KSVD *ctx = (SVD_KSVD*)svd->data;
249: PetscFunctionBegin;
250: ctx->polar = polar;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: /*@
255: SVDKSVDSetPolarMethod - Sets the method to compute the polar decomposition within the KSVD solver.
257: Logically Collective
259: Input Parameters:
260: + svd - the singular value solver context
261: - polar - method that will be used by KSVD for the polar decomposition
263: Options Database Key:
264: . -svd_ksvd_polar_method - Sets the method for the KSVD polar decomposition
266: If not set, the method defaults to SVD_KSVD_POLAR_QDWH.
268: Level: advanced
270: .seealso: SVDKSVDGetPolarMethod(), SVDKSVDPolarMethod
271: @*/
272: PetscErrorCode SVDKSVDSetPolarMethod(SVD svd,SVDKSVDPolarMethod polar)
273: {
274: PetscFunctionBegin;
277: PetscTryMethod(svd,"SVDKSVDSetPolarMethod_C",(SVD,SVDKSVDPolarMethod),(svd,polar));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: static PetscErrorCode SVDKSVDGetPolarMethod_KSVD(SVD svd,SVDKSVDPolarMethod *polar)
282: {
283: SVD_KSVD *ctx = (SVD_KSVD*)svd->data;
285: PetscFunctionBegin;
286: *polar = ctx->polar;
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: /*@
291: SVDKSVDGetPolarMethod - Gets the method for the KSVD polar decomposition.
293: Not Collective
295: Input Parameter:
296: . svd - the singular value solver context
298: Output Parameter:
299: . polar - method that will be used by KSVD for the polar decomposition
301: Level: advanced
303: .seealso: SVDKSVDSetPolarMethod(), SVDKSVDPolarMethod
304: @*/
305: PetscErrorCode SVDKSVDGetPolarMethod(SVD svd,SVDKSVDPolarMethod *polar)
306: {
307: PetscFunctionBegin;
309: PetscAssertPointer(polar,2);
310: PetscUseMethod(svd,"SVDKSVDGetPolarMethod_C",(SVD,SVDKSVDPolarMethod*),(svd,polar));
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: SLEPC_EXTERN PetscErrorCode SVDCreate_KSVD(SVD svd)
315: {
316: SVD_KSVD *ctx;
318: PetscFunctionBegin;
319: PetscCall(PetscNew(&ctx));
320: svd->data = (void*)ctx;
322: svd->ops->solve = SVDSolve_KSVD;
323: svd->ops->setup = SVDSetUp_KSVD;
324: svd->ops->setfromoptions = SVDSetFromOptions_KSVD;
325: svd->ops->destroy = SVDDestroy_KSVD;
326: svd->ops->reset = SVDReset_KSVD;
327: svd->ops->view = SVDView_KSVD;
329: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetEigenMethod_C",SVDKSVDSetEigenMethod_KSVD));
330: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetEigenMethod_C",SVDKSVDGetEigenMethod_KSVD));
331: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDSetPolarMethod_C",SVDKSVDSetPolarMethod_KSVD));
332: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDKSVDGetPolarMethod_C",SVDKSVDGetPolarMethod_KSVD));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }