Actual source code: cross.c
slepc-3.17.2 2022-08-09
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;
156: EPSProblemType ptype;
158: if (!cross->eps) SVDCrossGetEPS(svd,&cross->eps);
159: MatDestroy(&cross->C);
160: MatDestroy(&cross->D);
161: SVDCrossGetProductMat(svd,svd->A,svd->AT,&cross->C);
162: if (svd->isgeneralized) {
163: SVDCrossGetProductMat(svd,svd->B,svd->BT,&cross->D);
164: EPSSetOperators(cross->eps,cross->C,cross->D);
165: EPSGetProblemType(cross->eps,&ptype);
166: if (!ptype) EPSSetProblemType(cross->eps,EPS_GHEP);
167: } else {
168: EPSSetOperators(cross->eps,cross->C,NULL);
169: EPSSetProblemType(cross->eps,EPS_HEP);
170: }
171: if (!cross->usereps) {
172: EPSGetST(cross->eps,&st);
173: if (svd->isgeneralized && svd->which==SVD_SMALLEST) {
174: STSetType(st,STSINVERT);
175: EPSSetTarget(cross->eps,0.0);
176: EPSSetWhichEigenpairs(cross->eps,EPS_TARGET_REAL);
177: } else {
178: PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv);
179: if (issinv) EPSSetWhichEigenpairs(cross->eps,EPS_TARGET_MAGNITUDE);
180: else EPSSetWhichEigenpairs(cross->eps,svd->which==SVD_LARGEST?EPS_LARGEST_REAL:EPS_SMALLEST_REAL);
181: }
182: EPSSetDimensions(cross->eps,svd->nsv,svd->ncv,svd->mpd);
183: EPSSetTolerances(cross->eps,svd->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it);
184: switch (svd->conv) {
185: case SVD_CONV_ABS:
186: EPSSetConvergenceTest(cross->eps,EPS_CONV_ABS);break;
187: case SVD_CONV_REL:
188: EPSSetConvergenceTest(cross->eps,EPS_CONV_REL);break;
189: case SVD_CONV_NORM:
190: if (svd->isgeneralized) {
191: if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
192: if (!svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
193: EPSSetConvergenceTestFunction(cross->eps,EPSConv_Cross,svd,NULL);
194: } else {
195: EPSSetConvergenceTest(cross->eps,EPS_CONV_NORM);break;
196: }
197: break;
198: case SVD_CONV_MAXIT:
199: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
200: case SVD_CONV_USER:
201: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
202: }
203: }
204: SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
205: /* Transfer the trackall option from svd to eps */
206: SVDGetTrackAll(svd,&trackall);
207: EPSSetTrackAll(cross->eps,trackall);
208: /* Transfer the initial space from svd to eps */
209: if (svd->nini<0) {
210: EPSSetInitialSpace(cross->eps,-svd->nini,svd->IS);
211: SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
212: }
213: EPSSetUp(cross->eps);
214: EPSGetDimensions(cross->eps,NULL,&svd->ncv,&svd->mpd);
215: EPSGetTolerances(cross->eps,NULL,&svd->max_it);
216: if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
218: svd->leftbasis = PETSC_FALSE;
219: SVDAllocateSolution(svd,0);
220: PetscFunctionReturn(0);
221: }
223: PetscErrorCode SVDSolve_Cross(SVD svd)
224: {
225: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
226: PetscInt i;
227: PetscScalar lambda;
228: PetscReal sigma;
230: EPSSolve(cross->eps);
231: EPSGetConverged(cross->eps,&svd->nconv);
232: EPSGetIterationNumber(cross->eps,&svd->its);
233: EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason);
234: for (i=0;i<svd->nconv;i++) {
235: EPSGetEigenvalue(cross->eps,i,&lambda,NULL);
236: sigma = PetscRealPart(lambda);
238: if (sigma<0.0) {
239: PetscInfo(svd,"Negative eigenvalue computed by EPS: %g, resetting to 0\n",(double)sigma);
240: sigma = 0.0;
241: }
242: svd->sigma[i] = PetscSqrtReal(sigma);
243: }
244: PetscFunctionReturn(0);
245: }
247: PetscErrorCode SVDComputeVectors_Cross(SVD svd)
248: {
249: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
250: PetscInt i,mloc,ploc;
251: Vec u,v,x,uv;
252: PetscScalar *dst,alpha,lambda;
253: const PetscScalar *src;
254: PetscReal nrm;
256: if (svd->isgeneralized) {
257: MatCreateVecs(svd->A,NULL,&u);
258: VecGetLocalSize(u,&mloc);
259: MatCreateVecs(svd->B,NULL,&v);
260: VecGetLocalSize(v,&ploc);
261: for (i=0;i<svd->nconv;i++) {
262: BVGetColumn(svd->V,i,&x);
263: EPSGetEigenpair(cross->eps,i,&lambda,NULL,x,NULL);
264: MatMult(svd->A,x,u); /* u_i*c_i/alpha = A*x_i */
265: VecNormalize(u,NULL);
266: MatMult(svd->B,x,v); /* v_i*s_i/alpha = B*x_i */
267: VecNormalize(v,&nrm); /* ||v||_2 = s_i/alpha */
268: alpha = 1.0/(PetscSqrtReal(1.0+PetscRealPart(lambda))*nrm); /* alpha=s_i/||v||_2 */
269: VecScale(x,alpha);
270: BVRestoreColumn(svd->V,i,&x);
271: /* copy [u;v] to U[i] */
272: BVGetColumn(svd->U,i,&uv);
273: VecGetArrayWrite(uv,&dst);
274: VecGetArrayRead(u,&src);
275: PetscArraycpy(dst,src,mloc);
276: VecRestoreArrayRead(u,&src);
277: VecGetArrayRead(v,&src);
278: PetscArraycpy(dst+mloc,src,ploc);
279: VecRestoreArrayRead(v,&src);
280: VecRestoreArrayWrite(uv,&dst);
281: BVRestoreColumn(svd->U,i,&uv);
282: }
283: VecDestroy(&v);
284: VecDestroy(&u);
285: } else {
286: for (i=0;i<svd->nconv;i++) {
287: BVGetColumn(svd->V,i,&v);
288: EPSGetEigenvector(cross->eps,i,v,NULL);
289: BVRestoreColumn(svd->V,i,&v);
290: }
291: SVDComputeVectors_Left(svd);
292: }
293: PetscFunctionReturn(0);
294: }
296: static PetscErrorCode EPSMonitor_Cross(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
297: {
298: PetscInt i;
299: SVD svd = (SVD)ctx;
300: PetscScalar er,ei;
302: for (i=0;i<PetscMin(nest,svd->ncv);i++) {
303: er = eigr[i]; ei = eigi[i];
304: STBackTransform(eps->st,1,&er,&ei);
305: svd->sigma[i] = PetscSqrtReal(PetscAbsReal(PetscRealPart(er)));
306: svd->errest[i] = errest[i];
307: }
308: SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest);
309: PetscFunctionReturn(0);
310: }
312: PetscErrorCode SVDSetFromOptions_Cross(PetscOptionItems *PetscOptionsObject,SVD svd)
313: {
314: PetscBool set,val;
315: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
316: ST st;
318: PetscOptionsHead(PetscOptionsObject,"SVD Cross Options");
320: PetscOptionsBool("-svd_cross_explicitmatrix","Use cross explicit matrix","SVDCrossSetExplicitMatrix",cross->explicitmatrix,&val,&set);
321: if (set) SVDCrossSetExplicitMatrix(svd,val);
323: PetscOptionsTail();
325: if (!cross->eps) SVDCrossGetEPS(svd,&cross->eps);
326: if (!cross->explicitmatrix && !cross->usereps) {
327: /* use as default an ST with shell matrix and Jacobi */
328: EPSGetST(cross->eps,&st);
329: STSetMatMode(st,ST_MATMODE_SHELL);
330: }
331: EPSSetFromOptions(cross->eps);
332: PetscFunctionReturn(0);
333: }
335: static PetscErrorCode SVDCrossSetExplicitMatrix_Cross(SVD svd,PetscBool explicitmatrix)
336: {
337: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
339: if (cross->explicitmatrix != explicitmatrix) {
340: cross->explicitmatrix = explicitmatrix;
341: svd->state = SVD_STATE_INITIAL;
342: }
343: PetscFunctionReturn(0);
344: }
346: /*@
347: SVDCrossSetExplicitMatrix - Indicate if the eigensolver operator A^T*A must
348: be computed explicitly.
350: Logically Collective on svd
352: Input Parameters:
353: + svd - singular value solver
354: - explicitmat - boolean flag indicating if A^T*A is built explicitly
356: Options Database Key:
357: . -svd_cross_explicitmatrix <boolean> - Indicates the boolean flag
359: Level: advanced
361: .seealso: SVDCrossGetExplicitMatrix()
362: @*/
363: PetscErrorCode SVDCrossSetExplicitMatrix(SVD svd,PetscBool explicitmat)
364: {
367: PetscTryMethod(svd,"SVDCrossSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
368: PetscFunctionReturn(0);
369: }
371: static PetscErrorCode SVDCrossGetExplicitMatrix_Cross(SVD svd,PetscBool *explicitmat)
372: {
373: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
375: *explicitmat = cross->explicitmatrix;
376: PetscFunctionReturn(0);
377: }
379: /*@
380: SVDCrossGetExplicitMatrix - Returns the flag indicating if A^T*A is built explicitly.
382: Not Collective
384: Input Parameter:
385: . svd - singular value solver
387: Output Parameter:
388: . explicitmat - the mode flag
390: Level: advanced
392: .seealso: SVDCrossSetExplicitMatrix()
393: @*/
394: PetscErrorCode SVDCrossGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
395: {
398: PetscUseMethod(svd,"SVDCrossGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
399: PetscFunctionReturn(0);
400: }
402: static PetscErrorCode SVDCrossSetEPS_Cross(SVD svd,EPS eps)
403: {
404: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
406: PetscObjectReference((PetscObject)eps);
407: EPSDestroy(&cross->eps);
408: cross->eps = eps;
409: cross->usereps = PETSC_TRUE;
410: PetscLogObjectParent((PetscObject)svd,(PetscObject)cross->eps);
411: svd->state = SVD_STATE_INITIAL;
412: PetscFunctionReturn(0);
413: }
415: /*@
416: SVDCrossSetEPS - Associate an eigensolver object (EPS) to the
417: singular value solver.
419: Collective on svd
421: Input Parameters:
422: + svd - singular value solver
423: - eps - the eigensolver object
425: Level: advanced
427: .seealso: SVDCrossGetEPS()
428: @*/
429: PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)
430: {
434: PetscTryMethod(svd,"SVDCrossSetEPS_C",(SVD,EPS),(svd,eps));
435: PetscFunctionReturn(0);
436: }
438: static PetscErrorCode SVDCrossGetEPS_Cross(SVD svd,EPS *eps)
439: {
440: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
442: if (!cross->eps) {
443: EPSCreate(PetscObjectComm((PetscObject)svd),&cross->eps);
444: PetscObjectIncrementTabLevel((PetscObject)cross->eps,(PetscObject)svd,1);
445: EPSSetOptionsPrefix(cross->eps,((PetscObject)svd)->prefix);
446: EPSAppendOptionsPrefix(cross->eps,"svd_cross_");
447: PetscLogObjectParent((PetscObject)svd,(PetscObject)cross->eps);
448: PetscObjectSetOptions((PetscObject)cross->eps,((PetscObject)svd)->options);
449: EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_REAL);
450: EPSMonitorSet(cross->eps,EPSMonitor_Cross,svd,NULL);
451: }
452: *eps = cross->eps;
453: PetscFunctionReturn(0);
454: }
456: /*@
457: SVDCrossGetEPS - Retrieve the eigensolver object (EPS) associated
458: to the singular value solver.
460: Not Collective
462: Input Parameter:
463: . svd - singular value solver
465: Output Parameter:
466: . eps - the eigensolver object
468: Level: advanced
470: .seealso: SVDCrossSetEPS()
471: @*/
472: PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)
473: {
476: PetscUseMethod(svd,"SVDCrossGetEPS_C",(SVD,EPS*),(svd,eps));
477: PetscFunctionReturn(0);
478: }
480: PetscErrorCode SVDView_Cross(SVD svd,PetscViewer viewer)
481: {
482: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
483: PetscBool isascii;
485: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
486: if (isascii) {
487: if (!cross->eps) SVDCrossGetEPS(svd,&cross->eps);
488: PetscViewerASCIIPrintf(viewer," %s matrix\n",cross->explicitmatrix?"explicit":"implicit");
489: PetscViewerASCIIPushTab(viewer);
490: EPSView(cross->eps,viewer);
491: PetscViewerASCIIPopTab(viewer);
492: }
493: PetscFunctionReturn(0);
494: }
496: PetscErrorCode SVDReset_Cross(SVD svd)
497: {
498: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
500: EPSReset(cross->eps);
501: MatDestroy(&cross->C);
502: MatDestroy(&cross->D);
503: PetscFunctionReturn(0);
504: }
506: PetscErrorCode SVDDestroy_Cross(SVD svd)
507: {
508: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
510: EPSDestroy(&cross->eps);
511: PetscFree(svd->data);
512: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",NULL);
513: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",NULL);
514: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",NULL);
515: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",NULL);
516: PetscFunctionReturn(0);
517: }
519: SLEPC_EXTERN PetscErrorCode SVDCreate_Cross(SVD svd)
520: {
521: SVD_CROSS *cross;
523: PetscNewLog(svd,&cross);
524: svd->data = (void*)cross;
526: svd->ops->solve = SVDSolve_Cross;
527: svd->ops->solveg = SVDSolve_Cross;
528: svd->ops->setup = SVDSetUp_Cross;
529: svd->ops->setfromoptions = SVDSetFromOptions_Cross;
530: svd->ops->destroy = SVDDestroy_Cross;
531: svd->ops->reset = SVDReset_Cross;
532: svd->ops->view = SVDView_Cross;
533: svd->ops->computevectors = SVDComputeVectors_Cross;
534: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",SVDCrossSetEPS_Cross);
535: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",SVDCrossGetEPS_Cross);
536: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",SVDCrossSetExplicitMatrix_Cross);
537: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",SVDCrossGetExplicitMatrix_Cross);
538: PetscFunctionReturn(0);
539: }