Actual source code: svd.c
petsc-3.8.4 2018-03-24
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: #if !defined(PETSC_USE_COMPLEX)
262: MatMultTranspose(jac->Vt,work,yred);
263: #else
264: MatMultHermitianTranspose(jac->Vt,work,yred);
265: #endif
266: VecPointwiseMult(work,work,jac->diag);
267: #if !defined(PETSC_USE_COMPLEX)
268: MatMultTranspose(jac->U,xred,work);
269: #else
270: MatMultHermitianTranspose(jac->U,xred,work);
271: #endif
272: PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
273: PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
274: return(0);
275: }
277: static PetscErrorCode PCReset_SVD(PC pc)
278: {
279: PC_SVD *jac = (PC_SVD*)pc->data;
283: MatDestroy(&jac->A);
284: MatDestroy(&jac->U);
285: MatDestroy(&jac->Vt);
286: VecDestroy(&jac->diag);
287: VecDestroy(&jac->work);
288: VecScatterDestroy(&jac->right2red);
289: VecScatterDestroy(&jac->left2red);
290: VecDestroy(&jac->rightred);
291: VecDestroy(&jac->leftred);
292: return(0);
293: }
295: /* -------------------------------------------------------------------------- */
296: /*
297: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
298: that was created with PCCreate_SVD().
300: Input Parameter:
301: . pc - the preconditioner context
303: Application Interface Routine: PCDestroy()
304: */
305: static PetscErrorCode PCDestroy_SVD(PC pc)
306: {
307: PC_SVD *jac = (PC_SVD*)pc->data;
311: PCReset_SVD(pc);
312: PetscViewerDestroy(&jac->monitor);
313: PetscFree(pc->data);
314: return(0);
315: }
317: static PetscErrorCode PCSetFromOptions_SVD(PetscOptionItems *PetscOptionsObject,PC pc)
318: {
320: PC_SVD *jac = (PC_SVD*)pc->data;
321: PetscBool flg,set;
324: PetscOptionsHead(PetscOptionsObject,"SVD options");
325: PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
326: PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
327: PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
328: if (set) { /* Should make PCSVDSetMonitor() */
329: if (flg && !jac->monitor) {
330: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);
331: } else if (!flg) {
332: PetscViewerDestroy(&jac->monitor);
333: }
334: }
335: PetscOptionsTail();
336: return(0);
337: }
339: static PetscErrorCode PCView_SVD(PC pc,PetscViewer viewer)
340: {
341: PC_SVD *svd = (PC_SVD*)pc->data;
343: PetscBool iascii;
346: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
347: if (iascii) {
348: PetscViewerASCIIPrintf(viewer," All singular values smaller than %g treated as zero\n",(double)svd->zerosing);
349: PetscViewerASCIIPrintf(viewer," Provided essential rank of the matrix %D (all other eigenvalues are zeroed)\n",svd->essrank);
350: }
351: return(0);
352: }
353: /* -------------------------------------------------------------------------- */
354: /*
355: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
356: and sets this as the private data within the generic preconditioning
357: context, PC, that was created within PCCreate().
359: Input Parameter:
360: . pc - the preconditioner context
362: Application Interface Routine: PCCreate()
363: */
365: /*MC
366: PCSVD - Use pseudo inverse defined by SVD of operator
368: Level: advanced
370: Concepts: SVD
372: Options Database:
373: - -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero
374: + -pc_svd_monitor Print information on the extreme singular values of the operator
376: Developer Note: This implementation automatically creates a redundant copy of the
377: matrix on each process and uses a sequential SVD solve. Why does it do this instead
378: of using the composable PCREDUNDANT object?
380: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
381: M*/
383: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
384: {
385: PC_SVD *jac;
389: /*
390: Creates the private data structure for this preconditioner and
391: attach it to the PC object.
392: */
393: PetscNewLog(pc,&jac);
394: jac->zerosing = 1.e-12;
395: pc->data = (void*)jac;
397: /*
398: Set the pointers for the functions that are provided above.
399: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
400: are called, they will automatically call these functions. Note we
401: choose not to provide a couple of these functions since they are
402: not needed.
403: */
404: pc->ops->apply = PCApply_SVD;
405: pc->ops->applytranspose = PCApplyTranspose_SVD;
406: pc->ops->setup = PCSetUp_SVD;
407: pc->ops->reset = PCReset_SVD;
408: pc->ops->destroy = PCDestroy_SVD;
409: pc->ops->setfromoptions = PCSetFromOptions_SVD;
410: pc->ops->view = PCView_SVD;
411: pc->ops->applyrichardson = 0;
412: return(0);
413: }