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