Actual source code: svd.c
petsc-3.11.4 2019-09-28
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;
52: MatCreateRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&redmat);
53: MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
54: MatDestroy(&redmat);
55: } else {
56: MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
57: }
58: if (!jac->diag) { /* assume square matrices */
59: MatCreateVecs(jac->A,&jac->diag,&jac->work);
60: }
61: if (!jac->U) {
62: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
63: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);
64: }
65: MatGetSize(jac->A,&n,NULL);
66: if (!n) {
67: PetscInfo(pc,"Matrix has zero rows, skipping svd\n");
68: return(0);
69: }
70: PetscBLASIntCast(n,&nb);
71: lwork = 5*nb;
72: PetscMalloc1(lwork,&work);
73: MatDenseGetArray(jac->A,&a);
74: MatDenseGetArray(jac->U,&u);
75: MatDenseGetArray(jac->Vt,&v);
76: VecGetArray(jac->diag,&d);
77: #if !defined(PETSC_USE_COMPLEX)
78: {
79: PetscBLASInt lierr;
80: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
81: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
82: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
83: PetscFPTrapPop();
84: }
85: #else
86: {
87: PetscBLASInt lierr;
88: PetscReal *rwork,*dd;
89: PetscMalloc1(5*nb,&rwork);
90: PetscMalloc1(nb,&dd);
91: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
92: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,dd,u,&nb,v,&nb,work,&lwork,rwork,&lierr));
93: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
94: PetscFree(rwork);
95: for (i=0; i<n; i++) d[i] = dd[i];
96: PetscFree(dd);
97: PetscFPTrapPop();
98: }
99: #endif
100: MatDenseRestoreArray(jac->A,&a);
101: MatDenseRestoreArray(jac->U,&u);
102: MatDenseRestoreArray(jac->Vt,&v);
103: for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
104: jac->nzero = n-1-i;
105: if (jac->monitor) {
106: PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
107: 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);
108: if (n >= 10) { /* print 5 smallest and 5 largest */
109: 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]));
110: 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]));
111: } else { /* print all singular values */
112: char buf[256],*p;
113: size_t left = sizeof(buf),used;
114: PetscInt thisline;
115: for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
116: PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));
117: left -= used;
118: p += used;
119: if (thisline > 4 || i==0) {
120: PetscViewerASCIIPrintf(jac->monitor," SVD: singular values:%s\n",buf);
121: p = buf;
122: thisline = 0;
123: }
124: }
125: }
126: PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
127: }
128: PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));
129: for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
130: for (; i<n; i++) d[i] = 0.0;
131: if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
132: PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
133: VecRestoreArray(jac->diag,&d);
134: #if defined(foo)
135: {
136: PetscViewer viewer;
137: PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);
138: MatView(jac->A,viewer);
139: MatView(jac->U,viewer);
140: MatView(jac->Vt,viewer);
141: VecView(jac->diag,viewer);
142: PetscViewerDestroy(viewer);
143: }
144: #endif
145: PetscFree(work);
146: return(0);
147: #endif
148: }
150: static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
151: {
152: PC_SVD *jac = (PC_SVD*)pc->data;
154: PetscMPIInt size;
157: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
158: *xred = NULL;
159: switch (side) {
160: case PC_LEFT:
161: if (size == 1) *xred = x;
162: else {
163: if (!jac->left2red) {VecScatterCreateToAll(x,&jac->left2red,&jac->leftred);}
164: if (amode & READ) {
165: VecScatterBegin(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
166: VecScatterEnd(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
167: }
168: *xred = jac->leftred;
169: }
170: break;
171: case PC_RIGHT:
172: if (size == 1) *xred = x;
173: else {
174: if (!jac->right2red) {VecScatterCreateToAll(x,&jac->right2red,&jac->rightred);}
175: if (amode & READ) {
176: VecScatterBegin(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
177: VecScatterEnd(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
178: }
179: *xred = jac->rightred;
180: }
181: break;
182: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
183: }
184: return(0);
185: }
187: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
188: {
189: PC_SVD *jac = (PC_SVD*)pc->data;
191: PetscMPIInt size;
194: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
195: switch (side) {
196: case PC_LEFT:
197: if (size != 1 && amode & WRITE) {
198: VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
199: VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
200: }
201: break;
202: case PC_RIGHT:
203: if (size != 1 && amode & WRITE) {
204: VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
205: VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
206: }
207: break;
208: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
209: }
210: *xred = NULL;
211: return(0);
212: }
214: /* -------------------------------------------------------------------------- */
215: /*
216: PCApply_SVD - Applies the SVD preconditioner to a vector.
218: Input Parameters:
219: . pc - the preconditioner context
220: . x - input vector
222: Output Parameter:
223: . y - output vector
225: Application Interface Routine: PCApply()
226: */
227: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
228: {
229: PC_SVD *jac = (PC_SVD*)pc->data;
230: Vec work = jac->work,xred,yred;
234: PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
235: PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
236: #if !defined(PETSC_USE_COMPLEX)
237: MatMultTranspose(jac->U,xred,work);
238: #else
239: MatMultHermitianTranspose(jac->U,xred,work);
240: #endif
241: VecPointwiseMult(work,work,jac->diag);
242: #if !defined(PETSC_USE_COMPLEX)
243: MatMultTranspose(jac->Vt,work,yred);
244: #else
245: MatMultHermitianTranspose(jac->Vt,work,yred);
246: #endif
247: PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
248: PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
249: return(0);
250: }
252: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
253: {
254: PC_SVD *jac = (PC_SVD*)pc->data;
255: Vec work = jac->work,xred,yred;
259: PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);
260: PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);
261: MatMult(jac->Vt,xred,work);
262: VecPointwiseMult(work,work,jac->diag);
263: MatMult(jac->U,work,yred);
264: PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
265: PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
266: return(0);
267: }
269: static PetscErrorCode PCReset_SVD(PC pc)
270: {
271: PC_SVD *jac = (PC_SVD*)pc->data;
275: MatDestroy(&jac->A);
276: MatDestroy(&jac->U);
277: MatDestroy(&jac->Vt);
278: VecDestroy(&jac->diag);
279: VecDestroy(&jac->work);
280: VecScatterDestroy(&jac->right2red);
281: VecScatterDestroy(&jac->left2red);
282: VecDestroy(&jac->rightred);
283: VecDestroy(&jac->leftred);
284: return(0);
285: }
287: /* -------------------------------------------------------------------------- */
288: /*
289: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
290: that was created with PCCreate_SVD().
292: Input Parameter:
293: . pc - the preconditioner context
295: Application Interface Routine: PCDestroy()
296: */
297: static PetscErrorCode PCDestroy_SVD(PC pc)
298: {
299: PC_SVD *jac = (PC_SVD*)pc->data;
303: PCReset_SVD(pc);
304: PetscViewerDestroy(&jac->monitor);
305: PetscFree(pc->data);
306: return(0);
307: }
309: static PetscErrorCode PCSetFromOptions_SVD(PetscOptionItems *PetscOptionsObject,PC pc)
310: {
312: PC_SVD *jac = (PC_SVD*)pc->data;
313: PetscBool flg,set;
316: PetscOptionsHead(PetscOptionsObject,"SVD options");
317: PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
318: PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
319: PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
320: if (set) { /* Should make PCSVDSetMonitor() */
321: if (flg && !jac->monitor) {
322: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);
323: } else if (!flg) {
324: PetscViewerDestroy(&jac->monitor);
325: }
326: }
327: PetscOptionsTail();
328: return(0);
329: }
331: static PetscErrorCode PCView_SVD(PC pc,PetscViewer viewer)
332: {
333: PC_SVD *svd = (PC_SVD*)pc->data;
335: PetscBool iascii;
338: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
339: if (iascii) {
340: PetscViewerASCIIPrintf(viewer," All singular values smaller than %g treated as zero\n",(double)svd->zerosing);
341: PetscViewerASCIIPrintf(viewer," Provided essential rank of the matrix %D (all other eigenvalues are zeroed)\n",svd->essrank);
342: }
343: return(0);
344: }
345: /* -------------------------------------------------------------------------- */
346: /*
347: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
348: and sets this as the private data within the generic preconditioning
349: context, PC, that was created within PCCreate().
351: Input Parameter:
352: . pc - the preconditioner context
354: Application Interface Routine: PCCreate()
355: */
357: /*MC
358: PCSVD - Use pseudo inverse defined by SVD of operator
360: Level: advanced
362: Concepts: SVD
364: Options Database:
365: - -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero
366: + -pc_svd_monitor Print information on the extreme singular values of the operator
368: Developer Note: 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: }