Actual source code: svd.c
petsc-3.6.1 2015-08-06
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(pc->pmat,&n,NULL);
68: PetscBLASIntCast(n,&nb);
69: lwork = 5*nb;
70: PetscMalloc1(lwork,&work);
71: MatDenseGetArray(jac->A,&a);
72: MatDenseGetArray(jac->U,&u);
73: MatDenseGetArray(jac->Vt,&v);
74: VecGetArray(jac->diag,&d);
75: #if !defined(PETSC_USE_COMPLEX)
76: {
77: PetscBLASInt lierr;
78: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
79: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
80: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
81: PetscFPTrapPop();
82: }
83: #else
84: {
85: PetscBLASInt lierr;
86: PetscReal *rwork,*dd;
87: PetscMalloc1(5*nb,&rwork);
88: PetscMalloc1(nb,&dd);
89: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
90: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,dd,u,&nb,v,&nb,work,&lwork,rwork,&lierr));
91: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
92: PetscFree(rwork);
93: for (i=0; i<nb; i++) d[i] = dd[i];
94: PetscFree(dd);
95: PetscFPTrapPop();
96: }
97: #endif
98: MatDenseRestoreArray(jac->A,&a);
99: MatDenseRestoreArray(jac->U,&u);
100: MatDenseRestoreArray(jac->Vt,&v);
101: for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
102: jac->nzero = n-1-i;
103: if (jac->monitor) {
104: PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
105: 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);
106: if (n >= 10) { /* print 5 smallest and 5 largest */
107: 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]));
108: 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]));
109: } else { /* print all singular values */
110: char buf[256],*p;
111: size_t left = sizeof(buf),used;
112: PetscInt thisline;
113: for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
114: PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));
115: left -= used;
116: p += used;
117: if (thisline > 4 || i==0) {
118: PetscViewerASCIIPrintf(jac->monitor," SVD: singular values:%s\n",buf);
119: p = buf;
120: thisline = 0;
121: }
122: }
123: }
124: PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
125: }
126: PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));
127: for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
128: for (; i<n; i++) d[i] = 0.0;
129: if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
130: PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
131: VecRestoreArray(jac->diag,&d);
132: #if defined(foo)
133: {
134: PetscViewer viewer;
135: PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);
136: MatView(jac->A,viewer);
137: MatView(jac->U,viewer);
138: MatView(jac->Vt,viewer);
139: VecView(jac->diag,viewer);
140: PetscViewerDestroy(viewer);
141: }
142: #endif
143: PetscFree(work);
144: return(0);
145: #endif
146: }
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: }
189: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
190: {
191: PC_SVD *jac = (PC_SVD*)pc->data;
193: PetscMPIInt size;
196: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
197: switch (side) {
198: case PC_LEFT:
199: if (size != 1 && amode & WRITE) {
200: VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
201: VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
202: }
203: break;
204: case PC_RIGHT:
205: if (size != 1 && amode & WRITE) {
206: VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
207: VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
208: }
209: break;
210: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
211: }
212: *xred = NULL;
213: return(0);
214: }
216: /* -------------------------------------------------------------------------- */
217: /*
218: PCApply_SVD - Applies the SVD preconditioner to a vector.
220: Input Parameters:
221: . pc - the preconditioner context
222: . x - input vector
224: Output Parameter:
225: . y - output vector
227: Application Interface Routine: PCApply()
228: */
231: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
232: {
233: PC_SVD *jac = (PC_SVD*)pc->data;
234: Vec work = jac->work,xred,yred;
238: PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
239: PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
240: #if !defined(PETSC_USE_COMPLEX)
241: MatMultTranspose(jac->U,xred,work);
242: #else
243: MatMultHermitianTranspose(jac->U,xred,work);
244: #endif
245: VecPointwiseMult(work,work,jac->diag);
246: #if !defined(PETSC_USE_COMPLEX)
247: MatMultTranspose(jac->Vt,work,yred);
248: #else
249: MatMultHermitianTranspose(jac->Vt,work,yred);
250: #endif
251: PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
252: PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
253: return(0);
254: }
258: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
259: {
260: PC_SVD *jac = (PC_SVD*)pc->data;
261: Vec work = jac->work,xred,yred;
265: PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);
266: PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);
267: #if !defined(PETSC_USE_COMPLEX)
268: MatMultTranspose(jac->Vt,work,yred);
269: #else
270: MatMultHermitianTranspose(jac->Vt,work,yred);
271: #endif
272: VecPointwiseMult(work,work,jac->diag);
273: #if !defined(PETSC_USE_COMPLEX)
274: MatMultTranspose(jac->U,xred,work);
275: #else
276: MatMultHermitianTranspose(jac->U,xred,work);
277: #endif
278: PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
279: PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
280: return(0);
281: }
285: static PetscErrorCode PCReset_SVD(PC pc)
286: {
287: PC_SVD *jac = (PC_SVD*)pc->data;
291: MatDestroy(&jac->A);
292: MatDestroy(&jac->U);
293: MatDestroy(&jac->Vt);
294: VecDestroy(&jac->diag);
295: VecDestroy(&jac->work);
296: VecScatterDestroy(&jac->right2red);
297: VecScatterDestroy(&jac->left2red);
298: VecDestroy(&jac->rightred);
299: VecDestroy(&jac->leftred);
300: return(0);
301: }
303: /* -------------------------------------------------------------------------- */
304: /*
305: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
306: that was created with PCCreate_SVD().
308: Input Parameter:
309: . pc - the preconditioner context
311: Application Interface Routine: PCDestroy()
312: */
315: static PetscErrorCode PCDestroy_SVD(PC pc)
316: {
317: PC_SVD *jac = (PC_SVD*)pc->data;
321: PCReset_SVD(pc);
322: PetscViewerDestroy(&jac->monitor);
323: PetscFree(pc->data);
324: return(0);
325: }
329: static PetscErrorCode PCSetFromOptions_SVD(PetscOptions *PetscOptionsObject,PC pc)
330: {
332: PC_SVD *jac = (PC_SVD*)pc->data;
333: PetscBool flg,set;
336: PetscOptionsHead(PetscOptionsObject,"SVD options");
337: PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
338: PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
339: PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
340: if (set) { /* Should make PCSVDSetMonitor() */
341: if (flg && !jac->monitor) {
342: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);
343: } else if (!flg) {
344: PetscViewerDestroy(&jac->monitor);
345: }
346: }
347: PetscOptionsTail();
348: return(0);
349: }
351: /* -------------------------------------------------------------------------- */
352: /*
353: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
354: and sets this as the private data within the generic preconditioning
355: context, PC, that was created within PCCreate().
357: Input Parameter:
358: . pc - the preconditioner context
360: Application Interface Routine: PCCreate()
361: */
363: /*MC
364: PCSVD - Use pseudo inverse defined by SVD of operator
366: Level: advanced
368: Concepts: SVD
370: Options Database:
371: - -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero
372: + -pc_svd_monitor Print information on the extreme singular values of the operator
374: Developer Note: This implementation automatically creates a redundant copy of the
375: matrix on each process and uses a sequential SVD solve. Why does it do this instead
376: of using the composable PCREDUNDANT object?
378: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
379: 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 = 0;
411: pc->ops->applyrichardson = 0;
412: return(0);
413: }