Actual source code: cross.c
slepc-3.17.0 2022-03-31
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: SLEPc singular value solver: "cross"
13: Method: Uses a Hermitian eigensolver for A^T*A
14: */
16: #include <slepc/private/svdimpl.h>
17: #include <slepc/private/epsimpl.h>
19: typedef struct {
20: PetscBool explicitmatrix;
21: EPS eps;
22: PetscBool usereps;
23: Mat C,D;
24: } SVD_CROSS;
26: typedef struct {
27: Mat A,AT;
28: Vec w,diag;
29: PetscBool swapped;
30: } SVD_CROSS_SHELL;
32: static PetscErrorCode MatMult_Cross(Mat B,Vec x,Vec y)
33: {
34: SVD_CROSS_SHELL *ctx;
36: MatShellGetContext(B,&ctx);
37: MatMult(ctx->A,x,ctx->w);
38: MatMult(ctx->AT,ctx->w,y);
39: PetscFunctionReturn(0);
40: }
42: static PetscErrorCode MatGetDiagonal_Cross(Mat B,Vec d)
43: {
44: SVD_CROSS_SHELL *ctx;
45: PetscMPIInt len;
46: PetscInt N,n,i,j,start,end,ncols;
47: PetscScalar *work1,*work2,*diag;
48: const PetscInt *cols;
49: const PetscScalar *vals;
51: MatShellGetContext(B,&ctx);
52: if (!ctx->diag) {
53: /* compute diagonal from rows and store in ctx->diag */
54: VecDuplicate(d,&ctx->diag);
55: MatGetSize(ctx->A,NULL,&N);
56: MatGetLocalSize(ctx->A,NULL,&n);
57: PetscCalloc2(N,&work1,N,&work2);
58: if (ctx->swapped) {
59: MatGetOwnershipRange(ctx->AT,&start,&end);
60: for (i=start;i<end;i++) {
61: MatGetRow(ctx->AT,i,&ncols,NULL,&vals);
62: for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
63: MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals);
64: }
65: } else {
66: MatGetOwnershipRange(ctx->A,&start,&end);
67: for (i=start;i<end;i++) {
68: MatGetRow(ctx->A,i,&ncols,&cols,&vals);
69: for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
70: MatRestoreRow(ctx->A,i,&ncols,&cols,&vals);
71: }
72: }
73: PetscMPIIntCast(N,&len);
74: MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B));
75: VecGetOwnershipRange(ctx->diag,&start,&end);
76: VecGetArrayWrite(ctx->diag,&diag);
77: for (i=start;i<end;i++) diag[i-start] = work2[i];
78: VecRestoreArrayWrite(ctx->diag,&diag);
79: PetscFree2(work1,work2);
80: }
81: VecCopy(ctx->diag,d);
82: PetscFunctionReturn(0);
83: }
85: static PetscErrorCode MatDestroy_Cross(Mat B)
86: {
87: SVD_CROSS_SHELL *ctx;
89: MatShellGetContext(B,&ctx);
90: VecDestroy(&ctx->w);
91: VecDestroy(&ctx->diag);
92: PetscFree(ctx);
93: PetscFunctionReturn(0);
94: }
96: static PetscErrorCode SVDCrossGetProductMat(SVD svd,Mat A,Mat AT,Mat *C)
97: {
98: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
99: SVD_CROSS_SHELL *ctx;
100: PetscInt n;
101: VecType vtype;
103: if (cross->explicitmatrix) {
104: if (svd->expltrans) { /* explicit transpose */
105: MatProductCreate(AT,A,NULL,C);
106: MatProductSetType(*C,MATPRODUCT_AB);
107: } else { /* implicit transpose */
108: #if defined(PETSC_USE_COMPLEX)
109: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Must use explicit transpose with complex scalars");
110: #else
111: if (!svd->swapped) {
112: MatProductCreate(A,A,NULL,C);
113: MatProductSetType(*C,MATPRODUCT_AtB);
114: } else {
115: MatProductCreate(AT,AT,NULL,C);
116: MatProductSetType(*C,MATPRODUCT_ABt);
117: }
118: #endif
119: }
120: MatProductSetFromOptions(*C);
121: MatProductSymbolic(*C);
122: MatProductNumeric(*C);
123: } else {
124: PetscNew(&ctx);
125: ctx->A = A;
126: ctx->AT = AT;
127: ctx->swapped = svd->swapped;
128: MatCreateVecs(A,NULL,&ctx->w);
129: PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->w);
130: MatGetLocalSize(A,NULL,&n);
131: MatCreateShell(PetscObjectComm((PetscObject)svd),n,n,PETSC_DETERMINE,PETSC_DETERMINE,(void*)ctx,C);
132: MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cross);
133: MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cross);
134: MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cross);
135: MatGetVecType(A,&vtype);
136: MatSetVecType(*C,vtype);
137: }
138: PetscLogObjectParent((PetscObject)svd,(PetscObject)*C);
139: PetscFunctionReturn(0);
140: }
142: /* Convergence test relative to the norm of R (used in GSVD only) */
143: static PetscErrorCode EPSConv_Cross(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
144: {
145: SVD svd = (SVD)ctx;
147: *errest = res/PetscMax(svd->nrma,svd->nrmb);
148: PetscFunctionReturn(0);
149: }
151: PetscErrorCode SVDSetUp_Cross(SVD svd)
152: {
153: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
154: ST st;
155: PetscBool trackall,issinv;
157: if (!cross->eps) SVDCrossGetEPS(svd,&cross->eps);
158: MatDestroy(&cross->C);
159: MatDestroy(&cross->D);
160: SVDCrossGetProductMat(svd,svd->A,svd->AT,&cross->C);
161: if (svd->isgeneralized) {
162: SVDCrossGetProductMat(svd,svd->B,svd->BT,&cross->D);
163: EPSSetOperators(cross->eps,cross->C,cross->D);
164: EPSSetProblemType(cross->eps,EPS_GHEP);
165: } else {
166: EPSSetOperators(cross->eps,cross->C,NULL);
167: EPSSetProblemType(cross->eps,EPS_HEP);
168: }
169: if (!cross->usereps) {
170: EPSGetST(cross->eps,&st);
171: if (svd->isgeneralized && svd->which==SVD_SMALLEST) {
172: STSetType(st,STSINVERT);
173: EPSSetTarget(cross->eps,0.0);
174: EPSSetWhichEigenpairs(cross->eps,EPS_TARGET_REAL);
175: } else {
176: PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv);
177: if (issinv) EPSSetWhichEigenpairs(cross->eps,EPS_TARGET_MAGNITUDE);
178: else EPSSetWhichEigenpairs(cross->eps,svd->which==SVD_LARGEST?EPS_LARGEST_REAL:EPS_SMALLEST_REAL);
179: }
180: EPSSetDimensions(cross->eps,svd->nsv,svd->ncv,svd->mpd);
181: EPSSetTolerances(cross->eps,svd->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it);
182: switch (svd->conv) {
183: case SVD_CONV_ABS:
184: EPSSetConvergenceTest(cross->eps,EPS_CONV_ABS);break;
185: case SVD_CONV_REL:
186: EPSSetConvergenceTest(cross->eps,EPS_CONV_REL);break;
187: case SVD_CONV_NORM:
188: if (svd->isgeneralized) {
189: if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
190: if (!svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
191: EPSSetConvergenceTestFunction(cross->eps,EPSConv_Cross,svd,NULL);
192: } else {
193: EPSSetConvergenceTest(cross->eps,EPS_CONV_NORM);break;
194: }
195: break;
196: case SVD_CONV_MAXIT:
197: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
198: case SVD_CONV_USER:
199: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
200: }
201: }
202: SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
203: /* Transfer the trackall option from svd to eps */
204: SVDGetTrackAll(svd,&trackall);
205: EPSSetTrackAll(cross->eps,trackall);
206: /* Transfer the initial space from svd to eps */
207: if (svd->nini<0) {
208: EPSSetInitialSpace(cross->eps,-svd->nini,svd->IS);
209: SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
210: }
211: EPSSetUp(cross->eps);
212: EPSGetDimensions(cross->eps,NULL,&svd->ncv,&svd->mpd);
213: EPSGetTolerances(cross->eps,NULL,&svd->max_it);
214: if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
216: svd->leftbasis = PETSC_FALSE;
217: SVDAllocateSolution(svd,0);
218: PetscFunctionReturn(0);
219: }
221: PetscErrorCode SVDSolve_Cross(SVD svd)
222: {
223: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
224: PetscInt i;
225: PetscScalar lambda;
226: PetscReal sigma;
228: EPSSolve(cross->eps);
229: EPSGetConverged(cross->eps,&svd->nconv);
230: EPSGetIterationNumber(cross->eps,&svd->its);
231: EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason);
232: for (i=0;i<svd->nconv;i++) {
233: EPSGetEigenvalue(cross->eps,i,&lambda,NULL);
234: sigma = PetscRealPart(lambda);
236: if (sigma<0.0) {
237: PetscInfo(svd,"Negative eigenvalue computed by EPS: %g, resetting to 0\n",(double)sigma);
238: sigma = 0.0;
239: }
240: svd->sigma[i] = PetscSqrtReal(sigma);
241: }
242: PetscFunctionReturn(0);
243: }
245: PetscErrorCode SVDComputeVectors_Cross(SVD svd)
246: {
247: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
248: PetscInt i,mloc,ploc;
249: Vec u,v,x,uv;
250: PetscScalar *dst,alpha,lambda;
251: const PetscScalar *src;
252: PetscReal nrm;
254: if (svd->isgeneralized) {
255: MatCreateVecs(svd->A,NULL,&u);
256: VecGetLocalSize(u,&mloc);
257: MatCreateVecs(svd->B,NULL,&v);
258: VecGetLocalSize(v,&ploc);
259: for (i=0;i<svd->nconv;i++) {
260: BVGetColumn(svd->V,i,&x);
261: EPSGetEigenpair(cross->eps,i,&lambda,NULL,x,NULL);
262: MatMult(svd->A,x,u); /* u_i*c_i/alpha = A*x_i */
263: VecNormalize(u,NULL);
264: MatMult(svd->B,x,v); /* v_i*s_i/alpha = B*x_i */
265: VecNormalize(v,&nrm); /* ||v||_2 = s_i/alpha */
266: alpha = 1.0/(PetscSqrtReal(1.0+PetscRealPart(lambda))*nrm); /* alpha=s_i/||v||_2 */
267: VecScale(x,alpha);
268: BVRestoreColumn(svd->V,i,&x);
269: /* copy [u;v] to U[i] */
270: BVGetColumn(svd->U,i,&uv);
271: VecGetArrayWrite(uv,&dst);
272: VecGetArrayRead(u,&src);
273: PetscArraycpy(dst,src,mloc);
274: VecRestoreArrayRead(u,&src);
275: VecGetArrayRead(v,&src);
276: PetscArraycpy(dst+mloc,src,ploc);
277: VecRestoreArrayRead(v,&src);
278: VecRestoreArrayWrite(uv,&dst);
279: BVRestoreColumn(svd->U,i,&uv);
280: }
281: VecDestroy(&v);
282: VecDestroy(&u);
283: } else {
284: for (i=0;i<svd->nconv;i++) {
285: BVGetColumn(svd->V,i,&v);
286: EPSGetEigenvector(cross->eps,i,v,NULL);
287: BVRestoreColumn(svd->V,i,&v);
288: }
289: SVDComputeVectors_Left(svd);
290: }
291: PetscFunctionReturn(0);
292: }
294: static PetscErrorCode EPSMonitor_Cross(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
295: {
296: PetscInt i;
297: SVD svd = (SVD)ctx;
298: PetscScalar er,ei;
300: for (i=0;i<PetscMin(nest,svd->ncv);i++) {
301: er = eigr[i]; ei = eigi[i];
302: STBackTransform(eps->st,1,&er,&ei);
303: svd->sigma[i] = PetscSqrtReal(PetscAbsReal(PetscRealPart(er)));
304: svd->errest[i] = errest[i];
305: }
306: SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest);
307: PetscFunctionReturn(0);
308: }
310: PetscErrorCode SVDSetFromOptions_Cross(PetscOptionItems *PetscOptionsObject,SVD svd)
311: {
312: PetscBool set,val;
313: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
314: ST st;
316: PetscOptionsHead(PetscOptionsObject,"SVD Cross Options");
318: PetscOptionsBool("-svd_cross_explicitmatrix","Use cross explicit matrix","SVDCrossSetExplicitMatrix",cross->explicitmatrix,&val,&set);
319: if (set) SVDCrossSetExplicitMatrix(svd,val);
321: PetscOptionsTail();
323: if (!cross->eps) SVDCrossGetEPS(svd,&cross->eps);
324: if (!cross->explicitmatrix && !cross->usereps) {
325: /* use as default an ST with shell matrix and Jacobi */
326: EPSGetST(cross->eps,&st);
327: STSetMatMode(st,ST_MATMODE_SHELL);
328: }
329: EPSSetFromOptions(cross->eps);
330: PetscFunctionReturn(0);
331: }
333: static PetscErrorCode SVDCrossSetExplicitMatrix_Cross(SVD svd,PetscBool explicitmatrix)
334: {
335: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
337: if (cross->explicitmatrix != explicitmatrix) {
338: cross->explicitmatrix = explicitmatrix;
339: svd->state = SVD_STATE_INITIAL;
340: }
341: PetscFunctionReturn(0);
342: }
344: /*@
345: SVDCrossSetExplicitMatrix - Indicate if the eigensolver operator A^T*A must
346: be computed explicitly.
348: Logically Collective on svd
350: Input Parameters:
351: + svd - singular value solver
352: - explicitmat - boolean flag indicating if A^T*A is built explicitly
354: Options Database Key:
355: . -svd_cross_explicitmatrix <boolean> - Indicates the boolean flag
357: Level: advanced
359: .seealso: SVDCrossGetExplicitMatrix()
360: @*/
361: PetscErrorCode SVDCrossSetExplicitMatrix(SVD svd,PetscBool explicitmat)
362: {
365: PetscTryMethod(svd,"SVDCrossSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
366: PetscFunctionReturn(0);
367: }
369: static PetscErrorCode SVDCrossGetExplicitMatrix_Cross(SVD svd,PetscBool *explicitmat)
370: {
371: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
373: *explicitmat = cross->explicitmatrix;
374: PetscFunctionReturn(0);
375: }
377: /*@
378: SVDCrossGetExplicitMatrix - Returns the flag indicating if A^T*A is built explicitly.
380: Not Collective
382: Input Parameter:
383: . svd - singular value solver
385: Output Parameter:
386: . explicitmat - the mode flag
388: Level: advanced
390: .seealso: SVDCrossSetExplicitMatrix()
391: @*/
392: PetscErrorCode SVDCrossGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
393: {
396: PetscUseMethod(svd,"SVDCrossGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
397: PetscFunctionReturn(0);
398: }
400: static PetscErrorCode SVDCrossSetEPS_Cross(SVD svd,EPS eps)
401: {
402: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
404: PetscObjectReference((PetscObject)eps);
405: EPSDestroy(&cross->eps);
406: cross->eps = eps;
407: cross->usereps = PETSC_TRUE;
408: PetscLogObjectParent((PetscObject)svd,(PetscObject)cross->eps);
409: svd->state = SVD_STATE_INITIAL;
410: PetscFunctionReturn(0);
411: }
413: /*@
414: SVDCrossSetEPS - Associate an eigensolver object (EPS) to the
415: singular value solver.
417: Collective on svd
419: Input Parameters:
420: + svd - singular value solver
421: - eps - the eigensolver object
423: Level: advanced
425: .seealso: SVDCrossGetEPS()
426: @*/
427: PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)
428: {
432: PetscTryMethod(svd,"SVDCrossSetEPS_C",(SVD,EPS),(svd,eps));
433: PetscFunctionReturn(0);
434: }
436: static PetscErrorCode SVDCrossGetEPS_Cross(SVD svd,EPS *eps)
437: {
438: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
440: if (!cross->eps) {
441: EPSCreate(PetscObjectComm((PetscObject)svd),&cross->eps);
442: PetscObjectIncrementTabLevel((PetscObject)cross->eps,(PetscObject)svd,1);
443: EPSSetOptionsPrefix(cross->eps,((PetscObject)svd)->prefix);
444: EPSAppendOptionsPrefix(cross->eps,"svd_cross_");
445: PetscLogObjectParent((PetscObject)svd,(PetscObject)cross->eps);
446: PetscObjectSetOptions((PetscObject)cross->eps,((PetscObject)svd)->options);
447: EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_REAL);
448: EPSMonitorSet(cross->eps,EPSMonitor_Cross,svd,NULL);
449: }
450: *eps = cross->eps;
451: PetscFunctionReturn(0);
452: }
454: /*@
455: SVDCrossGetEPS - Retrieve the eigensolver object (EPS) associated
456: to the singular value solver.
458: Not Collective
460: Input Parameter:
461: . svd - singular value solver
463: Output Parameter:
464: . eps - the eigensolver object
466: Level: advanced
468: .seealso: SVDCrossSetEPS()
469: @*/
470: PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)
471: {
474: PetscUseMethod(svd,"SVDCrossGetEPS_C",(SVD,EPS*),(svd,eps));
475: PetscFunctionReturn(0);
476: }
478: PetscErrorCode SVDView_Cross(SVD svd,PetscViewer viewer)
479: {
480: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
481: PetscBool isascii;
483: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
484: if (isascii) {
485: if (!cross->eps) SVDCrossGetEPS(svd,&cross->eps);
486: PetscViewerASCIIPrintf(viewer," %s matrix\n",cross->explicitmatrix?"explicit":"implicit");
487: PetscViewerASCIIPushTab(viewer);
488: EPSView(cross->eps,viewer);
489: PetscViewerASCIIPopTab(viewer);
490: }
491: PetscFunctionReturn(0);
492: }
494: PetscErrorCode SVDReset_Cross(SVD svd)
495: {
496: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
498: EPSReset(cross->eps);
499: MatDestroy(&cross->C);
500: MatDestroy(&cross->D);
501: PetscFunctionReturn(0);
502: }
504: PetscErrorCode SVDDestroy_Cross(SVD svd)
505: {
506: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
508: EPSDestroy(&cross->eps);
509: PetscFree(svd->data);
510: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",NULL);
511: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",NULL);
512: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",NULL);
513: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",NULL);
514: PetscFunctionReturn(0);
515: }
517: SLEPC_EXTERN PetscErrorCode SVDCreate_Cross(SVD svd)
518: {
519: SVD_CROSS *cross;
521: PetscNewLog(svd,&cross);
522: svd->data = (void*)cross;
524: svd->ops->solve = SVDSolve_Cross;
525: svd->ops->solveg = SVDSolve_Cross;
526: svd->ops->setup = SVDSetUp_Cross;
527: svd->ops->setfromoptions = SVDSetFromOptions_Cross;
528: svd->ops->destroy = SVDDestroy_Cross;
529: svd->ops->reset = SVDReset_Cross;
530: svd->ops->view = SVDView_Cross;
531: svd->ops->computevectors = SVDComputeVectors_Cross;
532: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",SVDCrossSetEPS_Cross);
533: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",SVDCrossGetEPS_Cross);
534: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",SVDCrossSetExplicitMatrix_Cross);
535: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",SVDCrossGetExplicitMatrix_Cross);
536: PetscFunctionReturn(0);
537: }