Actual source code: svd.c
petsc-3.4.5 2014-06-29
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: PetscInt M;
55: MatGetSize(pc->pmat,&M,NULL);
56: MatGetRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,M,MAT_INITIAL_MATRIX,&redmat);
57: MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
58: MatDestroy(&redmat);
59: } else {
60: MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
61: }
62: if (!jac->diag) { /* assume square matrices */
63: MatGetVecs(jac->A,&jac->diag,&jac->work);
64: }
65: if (!jac->U) {
66: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
67: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);
68: }
69: MatGetSize(pc->pmat,&n,NULL);
70: PetscBLASIntCast(n,&nb);
71: lwork = 5*nb;
72: PetscMalloc(lwork*sizeof(PetscScalar),&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: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex");
87: #endif
88: MatDenseRestoreArray(jac->A,&a);
89: MatDenseRestoreArray(jac->U,&u);
90: MatDenseRestoreArray(jac->Vt,&v);
91: for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
92: jac->nzero = n-1-i;
93: if (jac->monitor) {
94: PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
95: 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);
96: if (n >= 10) { /* print 5 smallest and 5 largest */
97: 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]));
98: 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]));
99: } else { /* print all singular values */
100: char buf[256],*p;
101: size_t left = sizeof(buf),used;
102: PetscInt thisline;
103: for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
104: PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));
105: left -= used;
106: p += used;
107: if (thisline > 4 || i==0) {
108: PetscViewerASCIIPrintf(jac->monitor," SVD: singular values:%s\n",buf);
109: p = buf;
110: thisline = 0;
111: }
112: }
113: }
114: PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
115: }
116: PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));
117: for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
118: for (; i<n; i++) d[i] = 0.0;
119: if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
120: PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
121: VecRestoreArray(jac->diag,&d);
122: #if defined(foo)
123: {
124: PetscViewer viewer;
125: PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);
126: MatView(jac->A,viewer);
127: MatView(jac->U,viewer);
128: MatView(jac->Vt,viewer);
129: VecView(jac->diag,viewer);
130: PetscViewerDestroy(viewer);
131: }
132: #endif
133: PetscFree(work);
134: return(0);
135: #endif
136: }
140: static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
141: {
142: PC_SVD *jac = (PC_SVD*)pc->data;
144: PetscMPIInt size;
147: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
148: *xred = NULL;
149: switch (side) {
150: case PC_LEFT:
151: if (size == 1) *xred = x;
152: else {
153: if (!jac->left2red) {VecScatterCreateToAll(x,&jac->left2red,&jac->leftred);}
154: if (amode & READ) {
155: VecScatterBegin(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
156: VecScatterEnd(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
157: }
158: *xred = jac->leftred;
159: }
160: break;
161: case PC_RIGHT:
162: if (size == 1) *xred = x;
163: else {
164: if (!jac->right2red) {VecScatterCreateToAll(x,&jac->right2red,&jac->rightred);}
165: if (amode & READ) {
166: VecScatterBegin(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
167: VecScatterEnd(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
168: }
169: *xred = jac->rightred;
170: }
171: break;
172: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
173: }
174: return(0);
175: }
179: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
180: {
181: PC_SVD *jac = (PC_SVD*)pc->data;
183: PetscMPIInt size;
186: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
187: switch (side) {
188: case PC_LEFT:
189: if (size != 1 && amode & WRITE) {
190: VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
191: VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
192: }
193: break;
194: case PC_RIGHT:
195: if (size != 1 && amode & WRITE) {
196: VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
197: VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
198: }
199: break;
200: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
201: }
202: *xred = NULL;
203: return(0);
204: }
206: /* -------------------------------------------------------------------------- */
207: /*
208: PCApply_SVD - Applies the SVD preconditioner to a vector.
210: Input Parameters:
211: . pc - the preconditioner context
212: . x - input vector
214: Output Parameter:
215: . y - output vector
217: Application Interface Routine: PCApply()
218: */
221: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
222: {
223: PC_SVD *jac = (PC_SVD*)pc->data;
224: Vec work = jac->work,xred,yred;
228: PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
229: PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
230: MatMultTranspose(jac->U,xred,work);
231: VecPointwiseMult(work,work,jac->diag);
232: MatMultTranspose(jac->Vt,work,yred);
233: PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
234: PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
235: return(0);
236: }
240: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
241: {
242: PC_SVD *jac = (PC_SVD*)pc->data;
243: Vec work = jac->work,xred,yred;
247: PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);
248: PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);
249: MatMultTranspose(jac->Vt,work,yred);
250: VecPointwiseMult(work,work,jac->diag);
251: MatMultTranspose(jac->U,xred,work);
252: PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
253: PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
254: return(0);
255: }
259: static PetscErrorCode PCReset_SVD(PC pc)
260: {
261: PC_SVD *jac = (PC_SVD*)pc->data;
265: MatDestroy(&jac->A);
266: MatDestroy(&jac->U);
267: MatDestroy(&jac->Vt);
268: VecDestroy(&jac->diag);
269: VecDestroy(&jac->work);
270: VecScatterDestroy(&jac->right2red);
271: VecScatterDestroy(&jac->left2red);
272: VecDestroy(&jac->rightred);
273: VecDestroy(&jac->leftred);
274: return(0);
275: }
277: /* -------------------------------------------------------------------------- */
278: /*
279: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
280: that was created with PCCreate_SVD().
282: Input Parameter:
283: . pc - the preconditioner context
285: Application Interface Routine: PCDestroy()
286: */
289: static PetscErrorCode PCDestroy_SVD(PC pc)
290: {
291: PC_SVD *jac = (PC_SVD*)pc->data;
295: PCReset_SVD(pc);
296: PetscViewerDestroy(&jac->monitor);
297: PetscFree(pc->data);
298: return(0);
299: }
303: static PetscErrorCode PCSetFromOptions_SVD(PC pc)
304: {
306: PC_SVD *jac = (PC_SVD*)pc->data;
307: PetscBool flg,set;
310: PetscOptionsHead("SVD options");
311: PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
312: PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
313: PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
314: if (set) { /* Should make PCSVDSetMonitor() */
315: if (flg && !jac->monitor) {
316: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);
317: } else if (!flg) {
318: PetscViewerDestroy(&jac->monitor);
319: }
320: }
321: PetscOptionsTail();
322: return(0);
323: }
325: /* -------------------------------------------------------------------------- */
326: /*
327: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
328: and sets this as the private data within the generic preconditioning
329: context, PC, that was created within PCCreate().
331: Input Parameter:
332: . pc - the preconditioner context
334: Application Interface Routine: PCCreate()
335: */
337: /*MC
338: PCSVD - Use pseudo inverse defined by SVD of operator
340: Level: advanced
342: Concepts: SVD
344: Options Database:
345: - -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero
346: + -pc_svd_monitor Print information on the extreme singular values of the operator
348: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
349: M*/
353: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
354: {
355: PC_SVD *jac;
359: /*
360: Creates the private data structure for this preconditioner and
361: attach it to the PC object.
362: */
363: PetscNewLog(pc,PC_SVD,&jac);
364: jac->zerosing = 1.e-12;
365: pc->data = (void*)jac;
367: /*
368: Set the pointers for the functions that are provided above.
369: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
370: are called, they will automatically call these functions. Note we
371: choose not to provide a couple of these functions since they are
372: not needed.
373: */
374: pc->ops->apply = PCApply_SVD;
375: pc->ops->applytranspose = PCApplyTranspose_SVD;
376: pc->ops->setup = PCSetUp_SVD;
377: pc->ops->reset = PCReset_SVD;
378: pc->ops->destroy = PCDestroy_SVD;
379: pc->ops->setfromoptions = PCSetFromOptions_SVD;
380: pc->ops->view = 0;
381: pc->ops->applyrichardson = 0;
382: return(0);
383: }