Actual source code: svd.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
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: */
37: static PetscErrorCode PCSetUp_SVD(PC pc)
38: {
39: #if defined(PETSC_MISSING_LAPACK_GESVD)
40: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable\nNot able to provide singular value estimates.");
41: #else
42: PC_SVD *jac = (PC_SVD*)pc->data;
44: PetscScalar *a,*u,*v,*d,*work;
45: PetscBLASInt nb,lwork;
46: PetscInt i,n;
47: PetscMPIInt size;
50: MatDestroy(&jac->A);
51: MPI_Comm_size(((PetscObject)pc->pmat)->comm,&size);
52: if (size > 1) {
53: Mat redmat;
54: MatCreateRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&redmat);
55: MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
56: MatDestroy(&redmat);
57: } else {
58: MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
59: }
60: if (!jac->diag) { /* assume square matrices */
61: MatCreateVecs(jac->A,&jac->diag,&jac->work);
62: }
63: if (!jac->U) {
64: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
65: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);
66: }
67: MatGetSize(jac->A,&n,NULL);
68: if (!n) {
69: PetscInfo(pc,"Matrix has zero rows, skipping svd\n");
70: return(0);
71: }
72: PetscBLASIntCast(n,&nb);
73: lwork = 5*nb;
74: PetscMalloc1(lwork,&work);
75: MatDenseGetArray(jac->A,&a);
76: MatDenseGetArray(jac->U,&u);
77: MatDenseGetArray(jac->Vt,&v);
78: VecGetArray(jac->diag,&d);
79: #if !defined(PETSC_USE_COMPLEX)
80: {
81: PetscBLASInt lierr;
82: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
83: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
84: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
85: PetscFPTrapPop();
86: }
87: #else
88: {
89: PetscBLASInt lierr;
90: PetscReal *rwork,*dd;
91: PetscMalloc1(5*nb,&rwork);
92: PetscMalloc1(nb,&dd);
93: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
94: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,dd,u,&nb,v,&nb,work,&lwork,rwork,&lierr));
95: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
96: PetscFree(rwork);
97: for (i=0; i<n; i++) d[i] = dd[i];
98: PetscFree(dd);
99: PetscFPTrapPop();
100: }
101: #endif
102: MatDenseRestoreArray(jac->A,&a);
103: MatDenseRestoreArray(jac->U,&u);
104: MatDenseRestoreArray(jac->Vt,&v);
105: for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
106: jac->nzero = n-1-i;
107: if (jac->monitor) {
108: PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
109: 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);
110: if (n >= 10) { /* print 5 smallest and 5 largest */
111: 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]));
112: 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]));
113: } else { /* print all singular values */
114: char buf[256],*p;
115: size_t left = sizeof(buf),used;
116: PetscInt thisline;
117: for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
118: PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));
119: left -= used;
120: p += used;
121: if (thisline > 4 || i==0) {
122: PetscViewerASCIIPrintf(jac->monitor," SVD: singular values:%s\n",buf);
123: p = buf;
124: thisline = 0;
125: }
126: }
127: }
128: PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
129: }
130: PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));
131: for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
132: for (; i<n; i++) d[i] = 0.0;
133: if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
134: PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
135: VecRestoreArray(jac->diag,&d);
136: #if defined(foo)
137: {
138: PetscViewer viewer;
139: PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);
140: MatView(jac->A,viewer);
141: MatView(jac->U,viewer);
142: MatView(jac->Vt,viewer);
143: VecView(jac->diag,viewer);
144: PetscViewerDestroy(viewer);
145: }
146: #endif
147: PetscFree(work);
148: return(0);
149: #endif
150: }
154: static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
155: {
156: PC_SVD *jac = (PC_SVD*)pc->data;
158: PetscMPIInt size;
161: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
162: *xred = NULL;
163: switch (side) {
164: case PC_LEFT:
165: if (size == 1) *xred = x;
166: else {
167: if (!jac->left2red) {VecScatterCreateToAll(x,&jac->left2red,&jac->leftred);}
168: if (amode & READ) {
169: VecScatterBegin(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
170: VecScatterEnd(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
171: }
172: *xred = jac->leftred;
173: }
174: break;
175: case PC_RIGHT:
176: if (size == 1) *xred = x;
177: else {
178: if (!jac->right2red) {VecScatterCreateToAll(x,&jac->right2red,&jac->rightred);}
179: if (amode & READ) {
180: VecScatterBegin(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
181: VecScatterEnd(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
182: }
183: *xred = jac->rightred;
184: }
185: break;
186: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
187: }
188: return(0);
189: }
193: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
194: {
195: PC_SVD *jac = (PC_SVD*)pc->data;
197: PetscMPIInt size;
200: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
201: switch (side) {
202: case PC_LEFT:
203: if (size != 1 && amode & WRITE) {
204: VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
205: VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
206: }
207: break;
208: case PC_RIGHT:
209: if (size != 1 && amode & WRITE) {
210: VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
211: VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
212: }
213: break;
214: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
215: }
216: *xred = NULL;
217: return(0);
218: }
220: /* -------------------------------------------------------------------------- */
221: /*
222: PCApply_SVD - Applies the SVD preconditioner to a vector.
224: Input Parameters:
225: . pc - the preconditioner context
226: . x - input vector
228: Output Parameter:
229: . y - output vector
231: Application Interface Routine: PCApply()
232: */
235: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
236: {
237: PC_SVD *jac = (PC_SVD*)pc->data;
238: Vec work = jac->work,xred,yred;
242: PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
243: PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
244: #if !defined(PETSC_USE_COMPLEX)
245: MatMultTranspose(jac->U,xred,work);
246: #else
247: MatMultHermitianTranspose(jac->U,xred,work);
248: #endif
249: VecPointwiseMult(work,work,jac->diag);
250: #if !defined(PETSC_USE_COMPLEX)
251: MatMultTranspose(jac->Vt,work,yred);
252: #else
253: MatMultHermitianTranspose(jac->Vt,work,yred);
254: #endif
255: PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
256: PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
257: return(0);
258: }
262: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
263: {
264: PC_SVD *jac = (PC_SVD*)pc->data;
265: Vec work = jac->work,xred,yred;
269: PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);
270: PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);
271: #if !defined(PETSC_USE_COMPLEX)
272: MatMultTranspose(jac->Vt,work,yred);
273: #else
274: MatMultHermitianTranspose(jac->Vt,work,yred);
275: #endif
276: VecPointwiseMult(work,work,jac->diag);
277: #if !defined(PETSC_USE_COMPLEX)
278: MatMultTranspose(jac->U,xred,work);
279: #else
280: MatMultHermitianTranspose(jac->U,xred,work);
281: #endif
282: PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
283: PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
284: return(0);
285: }
289: static PetscErrorCode PCReset_SVD(PC pc)
290: {
291: PC_SVD *jac = (PC_SVD*)pc->data;
295: MatDestroy(&jac->A);
296: MatDestroy(&jac->U);
297: MatDestroy(&jac->Vt);
298: VecDestroy(&jac->diag);
299: VecDestroy(&jac->work);
300: VecScatterDestroy(&jac->right2red);
301: VecScatterDestroy(&jac->left2red);
302: VecDestroy(&jac->rightred);
303: VecDestroy(&jac->leftred);
304: return(0);
305: }
307: /* -------------------------------------------------------------------------- */
308: /*
309: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
310: that was created with PCCreate_SVD().
312: Input Parameter:
313: . pc - the preconditioner context
315: Application Interface Routine: PCDestroy()
316: */
319: static PetscErrorCode PCDestroy_SVD(PC pc)
320: {
321: PC_SVD *jac = (PC_SVD*)pc->data;
325: PCReset_SVD(pc);
326: PetscViewerDestroy(&jac->monitor);
327: PetscFree(pc->data);
328: return(0);
329: }
333: static PetscErrorCode PCSetFromOptions_SVD(PetscOptionItems *PetscOptionsObject,PC pc)
334: {
336: PC_SVD *jac = (PC_SVD*)pc->data;
337: PetscBool flg,set;
340: PetscOptionsHead(PetscOptionsObject,"SVD options");
341: PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
342: PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
343: PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
344: if (set) { /* Should make PCSVDSetMonitor() */
345: if (flg && !jac->monitor) {
346: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);
347: } else if (!flg) {
348: PetscViewerDestroy(&jac->monitor);
349: }
350: }
351: PetscOptionsTail();
352: return(0);
353: }
355: /* -------------------------------------------------------------------------- */
356: /*
357: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
358: and sets this as the private data within the generic preconditioning
359: context, PC, that was created within PCCreate().
361: Input Parameter:
362: . pc - the preconditioner context
364: Application Interface Routine: PCCreate()
365: */
367: /*MC
368: PCSVD - Use pseudo inverse defined by SVD of operator
370: Level: advanced
372: Concepts: SVD
374: Options Database:
375: - -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero
376: + -pc_svd_monitor Print information on the extreme singular values of the operator
378: Developer Note: This implementation automatically creates a redundant copy of the
379: matrix on each process and uses a sequential SVD solve. Why does it do this instead
380: of using the composable PCREDUNDANT object?
382: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
383: M*/
387: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
388: {
389: PC_SVD *jac;
393: /*
394: Creates the private data structure for this preconditioner and
395: attach it to the PC object.
396: */
397: PetscNewLog(pc,&jac);
398: jac->zerosing = 1.e-12;
399: pc->data = (void*)jac;
401: /*
402: Set the pointers for the functions that are provided above.
403: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
404: are called, they will automatically call these functions. Note we
405: choose not to provide a couple of these functions since they are
406: not needed.
407: */
408: pc->ops->apply = PCApply_SVD;
409: pc->ops->applytranspose = PCApplyTranspose_SVD;
410: pc->ops->setup = PCSetUp_SVD;
411: pc->ops->reset = PCReset_SVD;
412: pc->ops->destroy = PCDestroy_SVD;
413: pc->ops->setfromoptions = PCSetFromOptions_SVD;
414: pc->ops->view = 0;
415: pc->ops->applyrichardson = 0;
416: return(0);
417: }