Actual source code: svd.c
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: PetscViewerFormat monitorformat;
18: } PC_SVD;
20: typedef enum {READ=1, WRITE=2, READ_WRITE=3} AccessMode;
22: /* -------------------------------------------------------------------------- */
23: /*
24: PCSetUp_SVD - Prepares for the use of the SVD preconditioner
25: by setting data structures and options.
27: Input Parameter:
28: . pc - the preconditioner context
30: Application Interface Routine: PCSetUp()
32: Notes:
33: The interface routine PCSetUp() is not usually called directly by
34: the user, but instead is called by PCApply() if necessary.
35: */
36: static PetscErrorCode PCSetUp_SVD(PC pc)
37: {
38: PC_SVD *jac = (PC_SVD*)pc->data;
39: PetscScalar *a,*u,*v,*d,*work;
40: PetscBLASInt nb,lwork;
41: PetscInt i,n;
42: PetscMPIInt size;
44: MatDestroy(&jac->A);
45: MPI_Comm_size(((PetscObject)pc->pmat)->comm,&size);
46: if (size > 1) {
47: Mat redmat;
49: MatCreateRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&redmat);
50: MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
51: MatDestroy(&redmat);
52: } else {
53: MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
54: }
55: if (!jac->diag) { /* assume square matrices */
56: MatCreateVecs(jac->A,&jac->diag,&jac->work);
57: }
58: if (!jac->U) {
59: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
60: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);
61: }
62: MatGetSize(jac->A,&n,NULL);
63: if (!n) {
64: PetscInfo(pc,"Matrix has zero rows, skipping svd\n");
65: return 0;
66: }
67: PetscBLASIntCast(n,&nb);
68: lwork = 5*nb;
69: PetscMalloc1(lwork,&work);
70: MatDenseGetArray(jac->A,&a);
71: MatDenseGetArray(jac->U,&u);
72: MatDenseGetArray(jac->Vt,&v);
73: VecGetArray(jac->diag,&d);
74: #if !defined(PETSC_USE_COMPLEX)
75: {
76: PetscBLASInt lierr;
77: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
78: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
80: PetscFPTrapPop();
81: }
82: #else
83: {
84: PetscBLASInt lierr;
85: PetscReal *rwork,*dd;
86: PetscMalloc1(5*nb,&rwork);
87: PetscMalloc1(nb,&dd);
88: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
89: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,dd,u,&nb,v,&nb,work,&lwork,rwork,&lierr));
91: PetscFree(rwork);
92: for (i=0; i<n; i++) d[i] = dd[i];
93: PetscFree(dd);
94: PetscFPTrapPop();
95: }
96: #endif
97: MatDenseRestoreArray(jac->A,&a);
98: MatDenseRestoreArray(jac->U,&u);
99: MatDenseRestoreArray(jac->Vt,&v);
100: for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
101: jac->nzero = n-1-i;
102: if (jac->monitor) {
103: PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
104: 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);
105: if (n < 10 || jac->monitorformat == PETSC_VIEWER_ALL) {
106: PetscViewerASCIIPrintf(jac->monitor," SVD: singular values:\n");
107: for (i=0; i<n; i++) {
108: if (i%5 == 0) {
109: if (i != 0) {
110: PetscViewerASCIIPrintf(jac->monitor,"\n");
111: }
112: PetscViewerASCIIPrintf(jac->monitor," ");
113: }
114: PetscViewerASCIIPrintf(jac->monitor," %14.12e",(double)PetscRealPart(d[i]));
115: }
116: PetscViewerASCIIPrintf(jac->monitor,"\n");
117: } else { /* print 5 smallest and 5 largest */
118: 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]));
119: 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]));
120: }
121: PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
122: }
123: PetscInfo(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));
124: for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
125: for (; i<n; i++) d[i] = 0.0;
126: if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
127: PetscInfo(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
128: VecRestoreArray(jac->diag,&d);
129: PetscFree(work);
130: return 0;
131: }
133: static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
134: {
135: PC_SVD *jac = (PC_SVD*)pc->data;
136: PetscMPIInt size;
138: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
139: *xred = NULL;
140: switch (side) {
141: case PC_LEFT:
142: if (size == 1) *xred = x;
143: else {
144: if (!jac->left2red) VecScatterCreateToAll(x,&jac->left2red,&jac->leftred);
145: if (amode & READ) {
146: VecScatterBegin(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
147: VecScatterEnd(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
148: }
149: *xred = jac->leftred;
150: }
151: break;
152: case PC_RIGHT:
153: if (size == 1) *xred = x;
154: else {
155: if (!jac->right2red) VecScatterCreateToAll(x,&jac->right2red,&jac->rightred);
156: if (amode & READ) {
157: VecScatterBegin(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
158: VecScatterEnd(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
159: }
160: *xred = jac->rightred;
161: }
162: break;
163: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
164: }
165: return 0;
166: }
168: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
169: {
170: PC_SVD *jac = (PC_SVD*)pc->data;
171: PetscMPIInt size;
173: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
174: switch (side) {
175: case PC_LEFT:
176: if (size != 1 && amode & WRITE) {
177: VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
178: VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
179: }
180: break;
181: case PC_RIGHT:
182: if (size != 1 && amode & WRITE) {
183: VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
184: VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
185: }
186: break;
187: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
188: }
189: *xred = NULL;
190: return 0;
191: }
193: /* -------------------------------------------------------------------------- */
194: /*
195: PCApply_SVD - Applies the SVD preconditioner to a vector.
197: Input Parameters:
198: . pc - the preconditioner context
199: . x - input vector
201: Output Parameter:
202: . y - output vector
204: Application Interface Routine: PCApply()
205: */
206: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
207: {
208: PC_SVD *jac = (PC_SVD*)pc->data;
209: Vec work = jac->work,xred,yred;
211: PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
212: PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
213: #if !defined(PETSC_USE_COMPLEX)
214: MatMultTranspose(jac->U,xred,work);
215: #else
216: MatMultHermitianTranspose(jac->U,xred,work);
217: #endif
218: VecPointwiseMult(work,work,jac->diag);
219: #if !defined(PETSC_USE_COMPLEX)
220: MatMultTranspose(jac->Vt,work,yred);
221: #else
222: MatMultHermitianTranspose(jac->Vt,work,yred);
223: #endif
224: PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
225: PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
226: return 0;
227: }
229: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
230: {
231: PC_SVD *jac = (PC_SVD*)pc->data;
232: Vec work = jac->work,xred,yred;
234: PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);
235: PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);
236: MatMult(jac->Vt,xred,work);
237: VecPointwiseMult(work,work,jac->diag);
238: MatMult(jac->U,work,yred);
239: PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
240: PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
241: return 0;
242: }
244: static PetscErrorCode PCReset_SVD(PC pc)
245: {
246: PC_SVD *jac = (PC_SVD*)pc->data;
248: MatDestroy(&jac->A);
249: MatDestroy(&jac->U);
250: MatDestroy(&jac->Vt);
251: VecDestroy(&jac->diag);
252: VecDestroy(&jac->work);
253: VecScatterDestroy(&jac->right2red);
254: VecScatterDestroy(&jac->left2red);
255: VecDestroy(&jac->rightred);
256: VecDestroy(&jac->leftred);
257: return 0;
258: }
260: /* -------------------------------------------------------------------------- */
261: /*
262: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
263: that was created with PCCreate_SVD().
265: Input Parameter:
266: . pc - the preconditioner context
268: Application Interface Routine: PCDestroy()
269: */
270: static PetscErrorCode PCDestroy_SVD(PC pc)
271: {
272: PC_SVD *jac = (PC_SVD*)pc->data;
274: PCReset_SVD(pc);
275: PetscViewerDestroy(&jac->monitor);
276: PetscFree(pc->data);
277: return 0;
278: }
280: static PetscErrorCode PCSetFromOptions_SVD(PetscOptionItems *PetscOptionsObject,PC pc)
281: {
282: PC_SVD *jac = (PC_SVD*)pc->data;
283: PetscBool flg;
285: PetscOptionsHead(PetscOptionsObject,"SVD options");
286: PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
287: PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
288: PetscOptionsViewer("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",&jac->monitor,&jac->monitorformat,&flg);
289: PetscOptionsTail();
290: return 0;
291: }
293: static PetscErrorCode PCView_SVD(PC pc,PetscViewer viewer)
294: {
295: PC_SVD *svd = (PC_SVD*)pc->data;
296: PetscBool iascii;
298: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
299: if (iascii) {
300: PetscViewerASCIIPrintf(viewer," All singular values smaller than %g treated as zero\n",(double)svd->zerosing);
301: PetscViewerASCIIPrintf(viewer," Provided essential rank of the matrix %D (all other eigenvalues are zeroed)\n",svd->essrank);
302: }
303: return 0;
304: }
305: /* -------------------------------------------------------------------------- */
306: /*
307: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
308: and sets this as the private data within the generic preconditioning
309: context, PC, that was created within PCCreate().
311: Input Parameter:
312: . pc - the preconditioner context
314: Application Interface Routine: PCCreate()
315: */
317: /*MC
318: PCSVD - Use pseudo inverse defined by SVD of operator
320: Level: advanced
322: Options Database:
323: + -pc_svd_zero_sing <rtol> - Singular values smaller than this are treated as zero
324: - -pc_svd_monitor - Print information on the extreme singular values of the operator
326: Developer Note:
327: This implementation automatically creates a redundant copy of the
328: matrix on each process and uses a sequential SVD solve. Why does it do this instead
329: of using the composable PCREDUNDANT object?
331: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
332: M*/
334: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
335: {
336: PC_SVD *jac;
338: /*
339: Creates the private data structure for this preconditioner and
340: attach it to the PC object.
341: */
342: PetscNewLog(pc,&jac);
343: jac->zerosing = 1.e-12;
344: pc->data = (void*)jac;
346: /*
347: Set the pointers for the functions that are provided above.
348: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
349: are called, they will automatically call these functions. Note we
350: choose not to provide a couple of these functions since they are
351: not needed.
352: */
353: pc->ops->apply = PCApply_SVD;
354: pc->ops->applytranspose = PCApplyTranspose_SVD;
355: pc->ops->setup = PCSetUp_SVD;
356: pc->ops->reset = PCReset_SVD;
357: pc->ops->destroy = PCDestroy_SVD;
358: pc->ops->setfromoptions = PCSetFromOptions_SVD;
359: pc->ops->view = PCView_SVD;
360: pc->ops->applyrichardson = NULL;
361: return 0;
362: }