Actual source code: svd.c
petsc-3.12.5 2020-03-29
2: #include <petsc/private/pcimpl.h>
3: #include <petscblaslapack.h>
5: /*
6: Private context (data structure) for the SVD preconditioner.
7: */
8: typedef struct {
9: Vec diag,work;
10: Mat A,U,Vt;
11: PetscInt nzero;
12: PetscReal zerosing; /* measure of smallest singular value treated as nonzero */
13: PetscInt essrank; /* essential rank of operator */
14: VecScatter left2red,right2red;
15: Vec leftred,rightred;
16: PetscViewer monitor;
17: } PC_SVD;
19: typedef enum {READ=1, WRITE=2, READ_WRITE=3} AccessMode;
21: /* -------------------------------------------------------------------------- */
22: /*
23: PCSetUp_SVD - Prepares for the use of the SVD preconditioner
24: by setting data structures and options.
26: Input Parameter:
27: . pc - the preconditioner context
29: Application Interface Routine: PCSetUp()
31: Notes:
32: The interface routine PCSetUp() is not usually called directly by
33: the user, but instead is called by PCApply() if necessary.
34: */
35: static PetscErrorCode PCSetUp_SVD(PC pc)
36: {
37: #if defined(PETSC_MISSING_LAPACK_GESVD)
38: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable\nNot able to provide singular value estimates.");
39: #else
40: PC_SVD *jac = (PC_SVD*)pc->data;
42: PetscScalar *a,*u,*v,*d,*work;
43: PetscBLASInt nb,lwork;
44: PetscInt i,n;
45: PetscMPIInt size;
48: MatDestroy(&jac->A);
49: MPI_Comm_size(((PetscObject)pc->pmat)->comm,&size);
50: if (size > 1) {
51: Mat redmat;
53: MatCreateRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&redmat);
54: MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
55: MatDestroy(&redmat);
56: } else {
57: MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
58: }
59: if (!jac->diag) { /* assume square matrices */
60: MatCreateVecs(jac->A,&jac->diag,&jac->work);
61: }
62: if (!jac->U) {
63: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
64: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);
65: }
66: MatGetSize(jac->A,&n,NULL);
67: if (!n) {
68: PetscInfo(pc,"Matrix has zero rows, skipping svd\n");
69: return(0);
70: }
71: PetscBLASIntCast(n,&nb);
72: lwork = 5*nb;
73: PetscMalloc1(lwork,&work);
74: MatDenseGetArray(jac->A,&a);
75: MatDenseGetArray(jac->U,&u);
76: MatDenseGetArray(jac->Vt,&v);
77: VecGetArray(jac->diag,&d);
78: #if !defined(PETSC_USE_COMPLEX)
79: {
80: PetscBLASInt lierr;
81: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
82: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
83: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
84: PetscFPTrapPop();
85: }
86: #else
87: {
88: PetscBLASInt lierr;
89: PetscReal *rwork,*dd;
90: PetscMalloc1(5*nb,&rwork);
91: PetscMalloc1(nb,&dd);
92: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
93: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,dd,u,&nb,v,&nb,work,&lwork,rwork,&lierr));
94: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
95: PetscFree(rwork);
96: for (i=0; i<n; i++) d[i] = dd[i];
97: PetscFree(dd);
98: PetscFPTrapPop();
99: }
100: #endif
101: MatDenseRestoreArray(jac->A,&a);
102: MatDenseRestoreArray(jac->U,&u);
103: MatDenseRestoreArray(jac->Vt,&v);
104: for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
105: jac->nzero = n-1-i;
106: if (jac->monitor) {
107: PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
108: PetscViewerASCIIPrintf(jac->monitor," SVD: condition number %14.12e, %D of %D singular values are (nearly) zero\n",(double)PetscRealPart(d[0]/d[n-1]),jac->nzero,n);
109: if (n >= 10) { /* print 5 smallest and 5 largest */
110: PetscViewerASCIIPrintf(jac->monitor," SVD: smallest singular values: %14.12e %14.12e %14.12e %14.12e %14.12e\n",(double)PetscRealPart(d[n-1]),(double)PetscRealPart(d[n-2]),(double)PetscRealPart(d[n-3]),(double)PetscRealPart(d[n-4]),(double)PetscRealPart(d[n-5]));
111: PetscViewerASCIIPrintf(jac->monitor," SVD: largest singular values : %14.12e %14.12e %14.12e %14.12e %14.12e\n",(double)PetscRealPart(d[4]),(double)PetscRealPart(d[3]),(double)PetscRealPart(d[2]),(double)PetscRealPart(d[1]),(double)PetscRealPart(d[0]));
112: } else { /* print all singular values */
113: char buf[256],*p;
114: size_t left = sizeof(buf),used;
115: PetscInt thisline;
116: for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
117: PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));
118: left -= used;
119: p += used;
120: if (thisline > 4 || i==0) {
121: PetscViewerASCIIPrintf(jac->monitor," SVD: singular values:%s\n",buf);
122: p = buf;
123: thisline = 0;
124: }
125: }
126: }
127: PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
128: }
129: PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));
130: for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
131: for (; i<n; i++) d[i] = 0.0;
132: if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
133: PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
134: VecRestoreArray(jac->diag,&d);
135: #if defined(foo)
136: {
137: PetscViewer viewer;
138: PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);
139: MatView(jac->A,viewer);
140: MatView(jac->U,viewer);
141: MatView(jac->Vt,viewer);
142: VecView(jac->diag,viewer);
143: PetscViewerDestroy(viewer);
144: }
145: #endif
146: PetscFree(work);
147: return(0);
148: #endif
149: }
151: static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
152: {
153: PC_SVD *jac = (PC_SVD*)pc->data;
155: PetscMPIInt size;
158: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
159: *xred = NULL;
160: switch (side) {
161: case PC_LEFT:
162: if (size == 1) *xred = x;
163: else {
164: if (!jac->left2red) {VecScatterCreateToAll(x,&jac->left2red,&jac->leftred);}
165: if (amode & READ) {
166: VecScatterBegin(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
167: VecScatterEnd(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
168: }
169: *xred = jac->leftred;
170: }
171: break;
172: case PC_RIGHT:
173: if (size == 1) *xred = x;
174: else {
175: if (!jac->right2red) {VecScatterCreateToAll(x,&jac->right2red,&jac->rightred);}
176: if (amode & READ) {
177: VecScatterBegin(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
178: VecScatterEnd(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
179: }
180: *xred = jac->rightred;
181: }
182: break;
183: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
184: }
185: return(0);
186: }
188: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
189: {
190: PC_SVD *jac = (PC_SVD*)pc->data;
192: PetscMPIInt size;
195: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
196: switch (side) {
197: case PC_LEFT:
198: if (size != 1 && amode & WRITE) {
199: VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
200: VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
201: }
202: break;
203: case PC_RIGHT:
204: if (size != 1 && amode & WRITE) {
205: VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
206: VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
207: }
208: break;
209: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
210: }
211: *xred = NULL;
212: return(0);
213: }
215: /* -------------------------------------------------------------------------- */
216: /*
217: PCApply_SVD - Applies the SVD preconditioner to a vector.
219: Input Parameters:
220: . pc - the preconditioner context
221: . x - input vector
223: Output Parameter:
224: . y - output vector
226: Application Interface Routine: PCApply()
227: */
228: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
229: {
230: PC_SVD *jac = (PC_SVD*)pc->data;
231: Vec work = jac->work,xred,yred;
235: PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
236: PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
237: #if !defined(PETSC_USE_COMPLEX)
238: MatMultTranspose(jac->U,xred,work);
239: #else
240: MatMultHermitianTranspose(jac->U,xred,work);
241: #endif
242: VecPointwiseMult(work,work,jac->diag);
243: #if !defined(PETSC_USE_COMPLEX)
244: MatMultTranspose(jac->Vt,work,yred);
245: #else
246: MatMultHermitianTranspose(jac->Vt,work,yred);
247: #endif
248: PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
249: PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
250: return(0);
251: }
253: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
254: {
255: PC_SVD *jac = (PC_SVD*)pc->data;
256: Vec work = jac->work,xred,yred;
260: PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);
261: PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);
262: MatMult(jac->Vt,xred,work);
263: VecPointwiseMult(work,work,jac->diag);
264: MatMult(jac->U,work,yred);
265: PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
266: PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
267: return(0);
268: }
270: static PetscErrorCode PCReset_SVD(PC pc)
271: {
272: PC_SVD *jac = (PC_SVD*)pc->data;
276: MatDestroy(&jac->A);
277: MatDestroy(&jac->U);
278: MatDestroy(&jac->Vt);
279: VecDestroy(&jac->diag);
280: VecDestroy(&jac->work);
281: VecScatterDestroy(&jac->right2red);
282: VecScatterDestroy(&jac->left2red);
283: VecDestroy(&jac->rightred);
284: VecDestroy(&jac->leftred);
285: return(0);
286: }
288: /* -------------------------------------------------------------------------- */
289: /*
290: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
291: that was created with PCCreate_SVD().
293: Input Parameter:
294: . pc - the preconditioner context
296: Application Interface Routine: PCDestroy()
297: */
298: static PetscErrorCode PCDestroy_SVD(PC pc)
299: {
300: PC_SVD *jac = (PC_SVD*)pc->data;
304: PCReset_SVD(pc);
305: PetscViewerDestroy(&jac->monitor);
306: PetscFree(pc->data);
307: return(0);
308: }
310: static PetscErrorCode PCSetFromOptions_SVD(PetscOptionItems *PetscOptionsObject,PC pc)
311: {
313: PC_SVD *jac = (PC_SVD*)pc->data;
314: PetscBool flg,set;
317: PetscOptionsHead(PetscOptionsObject,"SVD options");
318: PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
319: PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
320: PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
321: if (set) { /* Should make PCSVDSetMonitor() */
322: if (flg && !jac->monitor) {
323: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);
324: } else if (!flg) {
325: PetscViewerDestroy(&jac->monitor);
326: }
327: }
328: PetscOptionsTail();
329: return(0);
330: }
332: static PetscErrorCode PCView_SVD(PC pc,PetscViewer viewer)
333: {
334: PC_SVD *svd = (PC_SVD*)pc->data;
336: PetscBool iascii;
339: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
340: if (iascii) {
341: PetscViewerASCIIPrintf(viewer," All singular values smaller than %g treated as zero\n",(double)svd->zerosing);
342: PetscViewerASCIIPrintf(viewer," Provided essential rank of the matrix %D (all other eigenvalues are zeroed)\n",svd->essrank);
343: }
344: return(0);
345: }
346: /* -------------------------------------------------------------------------- */
347: /*
348: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
349: and sets this as the private data within the generic preconditioning
350: context, PC, that was created within PCCreate().
352: Input Parameter:
353: . pc - the preconditioner context
355: Application Interface Routine: PCCreate()
356: */
358: /*MC
359: PCSVD - Use pseudo inverse defined by SVD of operator
361: Level: advanced
363: Options Database:
364: + -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero
365: - -pc_svd_monitor Print information on the extreme singular values of the operator
367: Developer Note:
368: This implementation automatically creates a redundant copy of the
369: matrix on each process and uses a sequential SVD solve. Why does it do this instead
370: of using the composable PCREDUNDANT object?
372: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
373: M*/
375: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
376: {
377: PC_SVD *jac;
381: /*
382: Creates the private data structure for this preconditioner and
383: attach it to the PC object.
384: */
385: PetscNewLog(pc,&jac);
386: jac->zerosing = 1.e-12;
387: pc->data = (void*)jac;
389: /*
390: Set the pointers for the functions that are provided above.
391: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
392: are called, they will automatically call these functions. Note we
393: choose not to provide a couple of these functions since they are
394: not needed.
395: */
396: pc->ops->apply = PCApply_SVD;
397: pc->ops->applytranspose = PCApplyTranspose_SVD;
398: pc->ops->setup = PCSetUp_SVD;
399: pc->ops->reset = PCReset_SVD;
400: pc->ops->destroy = PCDestroy_SVD;
401: pc->ops->setfromoptions = PCSetFromOptions_SVD;
402: pc->ops->view = PCView_SVD;
403: pc->ops->applyrichardson = 0;
404: return(0);
405: }