Actual source code: cross.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcsvd.h" I*/
17: #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
19: typedef struct {
20: PetscBool explicitmatrix;
21: EPS eps;
22: PetscBool usereps;
23: Mat mat;
24: Vec w,diag;
25: } SVD_CROSS;
27: static PetscErrorCode MatMult_Cross(Mat B,Vec x,Vec y)
28: {
30: SVD svd;
31: SVD_CROSS *cross;
34: MatShellGetContext(B,(void**)&svd);
35: cross = (SVD_CROSS*)svd->data;
36: SVDMatMult(svd,PETSC_FALSE,x,cross->w);
37: SVDMatMult(svd,PETSC_TRUE,cross->w,y);
38: return(0);
39: }
41: static PetscErrorCode MatCreateVecs_Cross(Mat B,Vec *right,Vec *left)
42: {
44: SVD svd;
47: MatShellGetContext(B,(void**)&svd);
48: if (right) {
49: SVDMatCreateVecs(svd,right,NULL);
50: if (left) { VecDuplicate(*right,left); }
51: } else {
52: SVDMatCreateVecs(svd,left,NULL);
53: }
54: return(0);
55: }
57: static PetscErrorCode MatGetDiagonal_Cross(Mat B,Vec d)
58: {
59: PetscErrorCode ierr;
60: SVD svd;
61: SVD_CROSS *cross;
62: PetscMPIInt len;
63: PetscInt N,n,i,j,start,end,ncols;
64: PetscScalar *work1,*work2,*diag;
65: const PetscInt *cols;
66: const PetscScalar *vals;
69: MatShellGetContext(B,(void**)&svd);
70: cross = (SVD_CROSS*)svd->data;
71: if (!cross->diag) {
72: /* compute diagonal from rows and store in cross->diag */
73: VecDuplicate(d,&cross->diag);
74: SVDMatGetSize(svd,NULL,&N);
75: SVDMatGetLocalSize(svd,NULL,&n);
76: PetscCalloc2(N,&work1,N,&work2);
77: if (svd->AT) {
78: MatGetOwnershipRange(svd->AT,&start,&end);
79: for (i=start;i<end;i++) {
80: MatGetRow(svd->AT,i,&ncols,NULL,&vals);
81: for (j=0;j<ncols;j++)
82: work1[i] += vals[j]*vals[j];
83: MatRestoreRow(svd->AT,i,&ncols,NULL,&vals);
84: }
85: } else {
86: MatGetOwnershipRange(svd->A,&start,&end);
87: for (i=start;i<end;i++) {
88: MatGetRow(svd->A,i,&ncols,&cols,&vals);
89: for (j=0;j<ncols;j++)
90: work1[cols[j]] += vals[j]*vals[j];
91: MatRestoreRow(svd->A,i,&ncols,&cols,&vals);
92: }
93: }
94: PetscMPIIntCast(N,&len);
95: MPI_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)svd));
96: VecGetOwnershipRange(cross->diag,&start,&end);
97: VecGetArray(cross->diag,&diag);
98: for (i=start;i<end;i++) diag[i-start] = work2[i];
99: VecRestoreArray(cross->diag,&diag);
100: PetscFree2(work1,work2);
101: }
102: VecCopy(cross->diag,d);
103: return(0);
104: }
106: PetscErrorCode SVDSetUp_Cross(SVD svd)
107: {
109: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
110: PetscInt n;
111: PetscBool trackall;
114: if (!cross->mat) {
115: if (cross->explicitmatrix) {
116: if (svd->A && svd->AT) { /* explicit transpose */
117: MatMatMult(svd->AT,svd->A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&cross->mat);
118: } else { /* implicit transpose */
119: #if defined(PETSC_USE_COMPLEX)
120: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Must use explicit transpose with complex scalars");
121: #else
122: if (svd->A) {
123: MatTransposeMatMult(svd->A,svd->A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&cross->mat);
124: } else {
125: MatMatTransposeMult(svd->AT,svd->AT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&cross->mat);
126: }
127: #endif
128: }
129: } else {
130: SVDMatGetLocalSize(svd,NULL,&n);
131: MatCreateShell(PetscObjectComm((PetscObject)svd),n,n,PETSC_DETERMINE,PETSC_DETERMINE,svd,&cross->mat);
132: MatShellSetOperation(cross->mat,MATOP_MULT,(void(*)(void))MatMult_Cross);
133: MatShellSetOperation(cross->mat,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_Cross);
134: MatShellSetOperation(cross->mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cross);
135: SVDMatCreateVecs(svd,NULL,&cross->w);
136: PetscLogObjectParent((PetscObject)svd,(PetscObject)cross->w);
137: }
138: PetscLogObjectParent((PetscObject)svd,(PetscObject)cross->mat);
139: }
141: if (!cross->eps) { SVDCrossGetEPS(svd,&cross->eps); }
142: EPSSetOperators(cross->eps,cross->mat,NULL);
143: EPSSetProblemType(cross->eps,EPS_HEP);
144: if (!cross->usereps) {
145: EPSSetWhichEigenpairs(cross->eps,svd->which==SVD_LARGEST?EPS_LARGEST_REAL:EPS_SMALLEST_REAL);
146: EPSSetDimensions(cross->eps,svd->nsv,svd->ncv?svd->ncv:PETSC_DEFAULT,svd->mpd?svd->mpd:PETSC_DEFAULT);
147: EPSSetTolerances(cross->eps,svd->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it?svd->max_it:PETSC_DEFAULT);
148: switch (svd->conv) {
149: case SVD_CONV_ABS:
150: EPSSetConvergenceTest(cross->eps,EPS_CONV_ABS);break;
151: case SVD_CONV_REL:
152: EPSSetConvergenceTest(cross->eps,EPS_CONV_REL);break;
153: case SVD_CONV_USER:
154: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
155: }
156: }
157: if (svd->stop!=SVD_STOP_BASIC) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined stopping test not supported in this solver");
158: /* Transfer the trackall option from svd to eps */
159: SVDGetTrackAll(svd,&trackall);
160: EPSSetTrackAll(cross->eps,trackall);
161: /* Transfer the initial space from svd to eps */
162: if (svd->nini<0) {
163: EPSSetInitialSpace(cross->eps,-svd->nini,svd->IS);
164: SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
165: }
166: EPSSetUp(cross->eps);
167: EPSGetDimensions(cross->eps,NULL,&svd->ncv,&svd->mpd);
168: EPSGetTolerances(cross->eps,NULL,&svd->max_it);
169: if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
171: svd->leftbasis = PETSC_FALSE;
172: SVDAllocateSolution(svd,0);
173: return(0);
174: }
176: PetscErrorCode SVDSolve_Cross(SVD svd)
177: {
179: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
180: PetscInt i;
181: PetscScalar lambda;
182: PetscReal sigma;
183: Vec v;
186: EPSSolve(cross->eps);
187: EPSGetConverged(cross->eps,&svd->nconv);
188: EPSGetIterationNumber(cross->eps,&svd->its);
189: EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason);
190: for (i=0;i<svd->nconv;i++) {
191: BVGetColumn(svd->V,i,&v);
192: EPSGetEigenpair(cross->eps,i,&lambda,NULL,v,NULL);
193: BVRestoreColumn(svd->V,i,&v);
194: sigma = PetscRealPart(lambda);
195: if (sigma<-10*PETSC_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject)svd),1,"Negative eigenvalue computed by EPS: %g",sigma);
196: if (sigma<0.0) {
197: PetscInfo1(svd,"Negative eigenvalue computed by EPS: %g, resetting to 0\n",sigma);
198: sigma = 0.0;
199: }
200: svd->sigma[i] = PetscSqrtReal(sigma);
201: }
202: return(0);
203: }
205: static PetscErrorCode EPSMonitor_Cross(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
206: {
207: PetscInt i;
208: SVD svd = (SVD)ctx;
209: PetscScalar er,ei;
213: for (i=0;i<PetscMin(nest,svd->ncv);i++) {
214: er = eigr[i]; ei = eigi[i];
215: STBackTransform(eps->st,1,&er,&ei);
216: svd->sigma[i] = PetscSqrtReal(PetscAbsReal(PetscRealPart(er)));
217: svd->errest[i] = errest[i];
218: }
219: SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest);
220: return(0);
221: }
223: PetscErrorCode SVDSetFromOptions_Cross(PetscOptionItems *PetscOptionsObject,SVD svd)
224: {
226: PetscBool set,val;
227: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
228: ST st;
231: PetscOptionsHead(PetscOptionsObject,"SVD Cross Options");
233: PetscOptionsBool("-svd_cross_explicitmatrix","Use cross explicit matrix","SVDCrossSetExplicitMatrix",cross->explicitmatrix,&val,&set);
234: if (set) { SVDCrossSetExplicitMatrix(svd,val); }
236: PetscOptionsTail();
238: if (!cross->eps) { SVDCrossGetEPS(svd,&cross->eps); }
239: if (!cross->explicitmatrix && !cross->usereps) {
240: /* use as default an ST with shell matrix and Jacobi */
241: EPSGetST(cross->eps,&st);
242: STSetMatMode(st,ST_MATMODE_SHELL);
243: }
244: EPSSetFromOptions(cross->eps);
245: return(0);
246: }
248: static PetscErrorCode SVDCrossSetExplicitMatrix_Cross(SVD svd,PetscBool explicitmatrix)
249: {
250: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
253: if (cross->explicitmatrix != explicitmatrix) {
254: cross->explicitmatrix = explicitmatrix;
255: svd->state = SVD_STATE_INITIAL;
256: }
257: return(0);
258: }
260: /*@
261: SVDCrossSetExplicitMatrix - Indicate if the eigensolver operator A^T*A must
262: be computed explicitly.
264: Logically Collective on SVD
266: Input Parameters:
267: + svd - singular value solver
268: - explicit - boolean flag indicating if A^T*A is built explicitly
270: Options Database Key:
271: . -svd_cross_explicitmatrix <boolean> - Indicates the boolean flag
273: Level: advanced
275: .seealso: SVDCrossGetExplicitMatrix()
276: @*/
277: PetscErrorCode SVDCrossSetExplicitMatrix(SVD svd,PetscBool explicitmatrix)
278: {
284: PetscTryMethod(svd,"SVDCrossSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmatrix));
285: return(0);
286: }
288: static PetscErrorCode SVDCrossGetExplicitMatrix_Cross(SVD svd,PetscBool *explicitmatrix)
289: {
290: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
293: *explicitmatrix = cross->explicitmatrix;
294: return(0);
295: }
297: /*@
298: SVDCrossGetExplicitMatrix - Returns the flag indicating if A^T*A is built explicitly.
300: Not Collective
302: Input Parameter:
303: . svd - singular value solver
305: Output Parameter:
306: . explicit - the mode flag
308: Level: advanced
310: .seealso: SVDCrossSetExplicitMatrix()
311: @*/
312: PetscErrorCode SVDCrossGetExplicitMatrix(SVD svd,PetscBool *explicitmatrix)
313: {
319: PetscUseMethod(svd,"SVDCrossGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmatrix));
320: return(0);
321: }
323: static PetscErrorCode SVDCrossSetEPS_Cross(SVD svd,EPS eps)
324: {
326: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
329: PetscObjectReference((PetscObject)eps);
330: EPSDestroy(&cross->eps);
331: cross->eps = eps;
332: cross->usereps = PETSC_TRUE;
333: PetscLogObjectParent((PetscObject)svd,(PetscObject)cross->eps);
334: svd->state = SVD_STATE_INITIAL;
335: return(0);
336: }
338: /*@
339: SVDCrossSetEPS - Associate an eigensolver object (EPS) to the
340: singular value solver.
342: Collective on SVD
344: Input Parameters:
345: + svd - singular value solver
346: - eps - the eigensolver object
348: Level: advanced
350: .seealso: SVDCrossGetEPS()
351: @*/
352: PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)
353: {
360: PetscTryMethod(svd,"SVDCrossSetEPS_C",(SVD,EPS),(svd,eps));
361: return(0);
362: }
364: static PetscErrorCode SVDCrossGetEPS_Cross(SVD svd,EPS *eps)
365: {
366: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
370: if (!cross->eps) {
371: EPSCreate(PetscObjectComm((PetscObject)svd),&cross->eps);
372: PetscObjectIncrementTabLevel((PetscObject)cross->eps,(PetscObject)svd,1);
373: EPSSetOptionsPrefix(cross->eps,((PetscObject)svd)->prefix);
374: EPSAppendOptionsPrefix(cross->eps,"svd_cross_");
375: PetscLogObjectParent((PetscObject)svd,(PetscObject)cross->eps);
376: PetscObjectSetOptions((PetscObject)cross->eps,((PetscObject)svd)->options);
377: EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_REAL);
378: EPSMonitorSet(cross->eps,EPSMonitor_Cross,svd,NULL);
379: }
380: *eps = cross->eps;
381: return(0);
382: }
384: /*@
385: SVDCrossGetEPS - Retrieve the eigensolver object (EPS) associated
386: to the singular value solver.
388: Not Collective
390: Input Parameter:
391: . svd - singular value solver
393: Output Parameter:
394: . eps - the eigensolver object
396: Level: advanced
398: .seealso: SVDCrossSetEPS()
399: @*/
400: PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)
401: {
407: PetscUseMethod(svd,"SVDCrossGetEPS_C",(SVD,EPS*),(svd,eps));
408: return(0);
409: }
411: PetscErrorCode SVDView_Cross(SVD svd,PetscViewer viewer)
412: {
414: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
415: PetscBool isascii;
418: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
419: if (isascii) {
420: if (!cross->eps) { SVDCrossGetEPS(svd,&cross->eps); }
421: PetscViewerASCIIPrintf(viewer," %s matrix\n",cross->explicitmatrix?"explicit":"implicit");
422: PetscViewerASCIIPushTab(viewer);
423: EPSView(cross->eps,viewer);
424: PetscViewerASCIIPopTab(viewer);
425: }
426: return(0);
427: }
429: PetscErrorCode SVDReset_Cross(SVD svd)
430: {
432: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
435: EPSReset(cross->eps);
436: MatDestroy(&cross->mat);
437: VecDestroy(&cross->w);
438: VecDestroy(&cross->diag);
439: return(0);
440: }
442: PetscErrorCode SVDDestroy_Cross(SVD svd)
443: {
445: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
448: EPSDestroy(&cross->eps);
449: PetscFree(svd->data);
450: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",NULL);
451: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",NULL);
452: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",NULL);
453: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",NULL);
454: return(0);
455: }
457: SLEPC_EXTERN PetscErrorCode SVDCreate_Cross(SVD svd)
458: {
460: SVD_CROSS *cross;
463: PetscNewLog(svd,&cross);
464: svd->data = (void*)cross;
466: svd->ops->solve = SVDSolve_Cross;
467: svd->ops->setup = SVDSetUp_Cross;
468: svd->ops->setfromoptions = SVDSetFromOptions_Cross;
469: svd->ops->destroy = SVDDestroy_Cross;
470: svd->ops->reset = SVDReset_Cross;
471: svd->ops->view = SVDView_Cross;
472: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",SVDCrossSetEPS_Cross);
473: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",SVDCrossGetEPS_Cross);
474: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",SVDCrossSetExplicitMatrix_Cross);
475: PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",SVDCrossGetExplicitMatrix_Cross);
476: return(0);
477: }