Actual source code: svd.c
petsc-3.14.6 2021-03-30
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: 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;
45: MatDestroy(&jac->A);
46: MPI_Comm_size(((PetscObject)pc->pmat)->comm,&size);
47: if (size > 1) {
48: Mat redmat;
50: MatCreateRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&redmat);
51: MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
52: MatDestroy(&redmat);
53: } else {
54: MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
55: }
56: if (!jac->diag) { /* assume square matrices */
57: MatCreateVecs(jac->A,&jac->diag,&jac->work);
58: }
59: if (!jac->U) {
60: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
61: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);
62: }
63: MatGetSize(jac->A,&n,NULL);
64: if (!n) {
65: PetscInfo(pc,"Matrix has zero rows, skipping svd\n");
66: return(0);
67: }
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<n; 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: }
147: static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
148: {
149: PC_SVD *jac = (PC_SVD*)pc->data;
151: PetscMPIInt size;
154: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
155: *xred = NULL;
156: switch (side) {
157: case PC_LEFT:
158: if (size == 1) *xred = x;
159: else {
160: if (!jac->left2red) {VecScatterCreateToAll(x,&jac->left2red,&jac->leftred);}
161: if (amode & READ) {
162: VecScatterBegin(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
163: VecScatterEnd(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
164: }
165: *xred = jac->leftred;
166: }
167: break;
168: case PC_RIGHT:
169: if (size == 1) *xred = x;
170: else {
171: if (!jac->right2red) {VecScatterCreateToAll(x,&jac->right2red,&jac->rightred);}
172: if (amode & READ) {
173: VecScatterBegin(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
174: VecScatterEnd(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
175: }
176: *xred = jac->rightred;
177: }
178: break;
179: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
180: }
181: return(0);
182: }
184: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
185: {
186: PC_SVD *jac = (PC_SVD*)pc->data;
188: PetscMPIInt size;
191: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
192: switch (side) {
193: case PC_LEFT:
194: if (size != 1 && amode & WRITE) {
195: VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
196: VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
197: }
198: break;
199: case PC_RIGHT:
200: if (size != 1 && amode & WRITE) {
201: VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
202: VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
203: }
204: break;
205: default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
206: }
207: *xred = NULL;
208: return(0);
209: }
211: /* -------------------------------------------------------------------------- */
212: /*
213: PCApply_SVD - Applies the SVD preconditioner to a vector.
215: Input Parameters:
216: . pc - the preconditioner context
217: . x - input vector
219: Output Parameter:
220: . y - output vector
222: Application Interface Routine: PCApply()
223: */
224: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
225: {
226: PC_SVD *jac = (PC_SVD*)pc->data;
227: Vec work = jac->work,xred,yred;
231: PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
232: PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
233: #if !defined(PETSC_USE_COMPLEX)
234: MatMultTranspose(jac->U,xred,work);
235: #else
236: MatMultHermitianTranspose(jac->U,xred,work);
237: #endif
238: VecPointwiseMult(work,work,jac->diag);
239: #if !defined(PETSC_USE_COMPLEX)
240: MatMultTranspose(jac->Vt,work,yred);
241: #else
242: MatMultHermitianTranspose(jac->Vt,work,yred);
243: #endif
244: PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
245: PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
246: return(0);
247: }
249: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
250: {
251: PC_SVD *jac = (PC_SVD*)pc->data;
252: Vec work = jac->work,xred,yred;
256: PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);
257: PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);
258: MatMult(jac->Vt,xred,work);
259: VecPointwiseMult(work,work,jac->diag);
260: MatMult(jac->U,work,yred);
261: PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
262: PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
263: return(0);
264: }
266: static PetscErrorCode PCReset_SVD(PC pc)
267: {
268: PC_SVD *jac = (PC_SVD*)pc->data;
272: MatDestroy(&jac->A);
273: MatDestroy(&jac->U);
274: MatDestroy(&jac->Vt);
275: VecDestroy(&jac->diag);
276: VecDestroy(&jac->work);
277: VecScatterDestroy(&jac->right2red);
278: VecScatterDestroy(&jac->left2red);
279: VecDestroy(&jac->rightred);
280: VecDestroy(&jac->leftred);
281: return(0);
282: }
284: /* -------------------------------------------------------------------------- */
285: /*
286: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
287: that was created with PCCreate_SVD().
289: Input Parameter:
290: . pc - the preconditioner context
292: Application Interface Routine: PCDestroy()
293: */
294: static PetscErrorCode PCDestroy_SVD(PC pc)
295: {
296: PC_SVD *jac = (PC_SVD*)pc->data;
300: PCReset_SVD(pc);
301: PetscViewerDestroy(&jac->monitor);
302: PetscFree(pc->data);
303: return(0);
304: }
306: static PetscErrorCode PCSetFromOptions_SVD(PetscOptionItems *PetscOptionsObject,PC pc)
307: {
309: PC_SVD *jac = (PC_SVD*)pc->data;
310: PetscBool flg,set;
313: PetscOptionsHead(PetscOptionsObject,"SVD options");
314: PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
315: PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
316: PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
317: if (set) { /* Should make PCSVDSetMonitor() */
318: if (flg && !jac->monitor) {
319: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);
320: } else if (!flg) {
321: PetscViewerDestroy(&jac->monitor);
322: }
323: }
324: PetscOptionsTail();
325: return(0);
326: }
328: static PetscErrorCode PCView_SVD(PC pc,PetscViewer viewer)
329: {
330: PC_SVD *svd = (PC_SVD*)pc->data;
332: PetscBool iascii;
335: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
336: if (iascii) {
337: PetscViewerASCIIPrintf(viewer," All singular values smaller than %g treated as zero\n",(double)svd->zerosing);
338: PetscViewerASCIIPrintf(viewer," Provided essential rank of the matrix %D (all other eigenvalues are zeroed)\n",svd->essrank);
339: }
340: return(0);
341: }
342: /* -------------------------------------------------------------------------- */
343: /*
344: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
345: and sets this as the private data within the generic preconditioning
346: context, PC, that was created within PCCreate().
348: Input Parameter:
349: . pc - the preconditioner context
351: Application Interface Routine: PCCreate()
352: */
354: /*MC
355: PCSVD - Use pseudo inverse defined by SVD of operator
357: Level: advanced
359: Options Database:
360: + -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero
361: - -pc_svd_monitor Print information on the extreme singular values of the operator
363: Developer Note:
364: This implementation automatically creates a redundant copy of the
365: matrix on each process and uses a sequential SVD solve. Why does it do this instead
366: of using the composable PCREDUNDANT object?
368: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
369: M*/
371: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
372: {
373: PC_SVD *jac;
377: /*
378: Creates the private data structure for this preconditioner and
379: attach it to the PC object.
380: */
381: PetscNewLog(pc,&jac);
382: jac->zerosing = 1.e-12;
383: pc->data = (void*)jac;
385: /*
386: Set the pointers for the functions that are provided above.
387: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
388: are called, they will automatically call these functions. Note we
389: choose not to provide a couple of these functions since they are
390: not needed.
391: */
392: pc->ops->apply = PCApply_SVD;
393: pc->ops->applytranspose = PCApplyTranspose_SVD;
394: pc->ops->setup = PCSetUp_SVD;
395: pc->ops->reset = PCReset_SVD;
396: pc->ops->destroy = PCDestroy_SVD;
397: pc->ops->setfromoptions = PCSetFromOptions_SVD;
398: pc->ops->view = PCView_SVD;
399: pc->ops->applyrichardson = NULL;
400: return(0);
401: }