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: PetscFunctionBegin;
48: PetscCall(MatDestroy(&jac->A));
49: PetscCallMPI(MPI_Comm_size(((PetscObject)pc->pmat)->comm, &size));
50: if (size > 1) {
51: Mat redmat;
53: PetscCall(MatCreateRedundantMatrix(pc->pmat, size, PETSC_COMM_SELF, MAT_INITIAL_MATRIX, &redmat));
54: PetscCall(MatConvert(redmat, MATSEQDENSE, MAT_INITIAL_MATRIX, &jac->A));
55: PetscCall(MatDestroy(&redmat));
56: } else {
57: PetscCall(MatConvert(pc->pmat, MATSEQDENSE, MAT_INITIAL_MATRIX, &jac->A));
58: }
59: if (!jac->diag) { /* assume square matrices */
60: PetscCall(MatCreateVecs(jac->A, &jac->diag, &jac->work));
61: }
62: if (!jac->U) {
63: PetscCall(MatDuplicate(jac->A, MAT_DO_NOT_COPY_VALUES, &jac->U));
64: PetscCall(MatDuplicate(jac->A, MAT_DO_NOT_COPY_VALUES, &jac->Vt));
65: }
66: PetscCall(MatGetSize(jac->A, &n, NULL));
67: if (!n) {
68: PetscCall(PetscInfo(pc, "Matrix has zero rows, skipping svd\n"));
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
71: PetscCall(PetscBLASIntCast(n, &nb));
72: lwork = 5 * nb;
73: PetscCall(PetscMalloc1(lwork, &work));
74: PetscCall(MatDenseGetArray(jac->A, &a));
75: PetscCall(MatDenseGetArray(jac->U, &u));
76: PetscCall(MatDenseGetArray(jac->Vt, &v));
77: PetscCall(VecGetArray(jac->diag, &d));
78: #if !defined(PETSC_USE_COMPLEX)
79: {
80: PetscBLASInt lierr;
81: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
82: PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &nb, &nb, a, &nb, d, u, &nb, v, &nb, work, &lwork, &lierr));
83: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "gesv() error %" PetscBLASInt_FMT, lierr);
84: PetscCall(PetscFPTrapPop());
85: }
86: #else
87: {
88: PetscBLASInt lierr;
89: PetscReal *rwork, *dd;
90: PetscCall(PetscMalloc1(5 * nb, &rwork));
91: PetscCall(PetscMalloc1(nb, &dd));
92: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
93: PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &nb, &nb, a, &nb, dd, u, &nb, v, &nb, work, &lwork, rwork, &lierr));
94: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "gesv() error %" PetscBLASInt_FMT, lierr);
95: PetscCall(PetscFree(rwork));
96: for (i = 0; i < n; i++) d[i] = dd[i];
97: PetscCall(PetscFree(dd));
98: PetscCall(PetscFPTrapPop());
99: }
100: #endif
101: PetscCall(MatDenseRestoreArray(jac->A, &a));
102: PetscCall(MatDenseRestoreArray(jac->U, &u));
103: PetscCall(MatDenseRestoreArray(jac->Vt, &v));
104: for (i = n - 1; i >= 0; i--)
105: if (PetscRealPart(d[i]) > jac->zerosing) break;
106: jac->nzero = n - 1 - i;
107: if (jac->monitor) {
108: PetscCall(PetscViewerASCIIAddTab(jac->monitor, ((PetscObject)pc)->tablevel));
109: PetscCall(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));
110: if (n < 10 || jac->monitorformat == PETSC_VIEWER_ALL) {
111: PetscCall(PetscViewerASCIIPrintf(jac->monitor, " SVD: singular values:\n"));
112: for (i = 0; i < n; i++) {
113: if (i % 5 == 0) {
114: if (i != 0) PetscCall(PetscViewerASCIIPrintf(jac->monitor, "\n"));
115: PetscCall(PetscViewerASCIIPrintf(jac->monitor, " "));
116: }
117: PetscCall(PetscViewerASCIIPrintf(jac->monitor, " %14.12e", (double)PetscRealPart(d[i])));
118: }
119: PetscCall(PetscViewerASCIIPrintf(jac->monitor, "\n"));
120: } else { /* print 5 smallest and 5 largest */
121: PetscCall(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])));
122: PetscCall(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])));
123: }
124: PetscCall(PetscViewerASCIISubtractTab(jac->monitor, ((PetscObject)pc)->tablevel));
125: }
126: PetscCall(PetscInfo(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)
130: for (i = 0; i < n - jac->nzero - jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
131: PetscCall(PetscInfo(pc, "Number of zero or nearly singular values %" PetscInt_FMT "\n", jac->nzero));
132: PetscCall(VecRestoreArray(jac->diag, &d));
133: PetscCall(PetscFree(work));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode PCSVDGetVec(PC pc, PCSide side, AccessMode amode, Vec x, Vec *xred)
138: {
139: PC_SVD *jac = (PC_SVD *)pc->data;
140: PetscMPIInt size;
142: PetscFunctionBegin;
143: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
144: *xred = NULL;
145: switch (side) {
146: case PC_LEFT:
147: if (size == 1) *xred = x;
148: else {
149: if (!jac->left2red) PetscCall(VecScatterCreateToAll(x, &jac->left2red, &jac->leftred));
150: if (amode & READ) {
151: PetscCall(VecScatterBegin(jac->left2red, x, jac->leftred, INSERT_VALUES, SCATTER_FORWARD));
152: PetscCall(VecScatterEnd(jac->left2red, x, jac->leftred, INSERT_VALUES, SCATTER_FORWARD));
153: }
154: *xred = jac->leftred;
155: }
156: break;
157: case PC_RIGHT:
158: if (size == 1) *xred = x;
159: else {
160: if (!jac->right2red) PetscCall(VecScatterCreateToAll(x, &jac->right2red, &jac->rightred));
161: if (amode & READ) {
162: PetscCall(VecScatterBegin(jac->right2red, x, jac->rightred, INSERT_VALUES, SCATTER_FORWARD));
163: PetscCall(VecScatterEnd(jac->right2red, x, jac->rightred, INSERT_VALUES, SCATTER_FORWARD));
164: }
165: *xred = jac->rightred;
166: }
167: break;
168: default:
169: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Side must be LEFT or RIGHT");
170: }
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: static PetscErrorCode PCSVDRestoreVec(PC pc, PCSide side, AccessMode amode, Vec x, Vec *xred)
175: {
176: PC_SVD *jac = (PC_SVD *)pc->data;
177: PetscMPIInt size;
179: PetscFunctionBegin;
180: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
181: switch (side) {
182: case PC_LEFT:
183: if (size != 1 && amode & WRITE) {
184: PetscCall(VecScatterBegin(jac->left2red, jac->leftred, x, INSERT_VALUES, SCATTER_REVERSE));
185: PetscCall(VecScatterEnd(jac->left2red, jac->leftred, x, INSERT_VALUES, SCATTER_REVERSE));
186: }
187: break;
188: case PC_RIGHT:
189: if (size != 1 && amode & WRITE) {
190: PetscCall(VecScatterBegin(jac->right2red, jac->rightred, x, INSERT_VALUES, SCATTER_REVERSE));
191: PetscCall(VecScatterEnd(jac->right2red, jac->rightred, x, INSERT_VALUES, SCATTER_REVERSE));
192: }
193: break;
194: default:
195: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Side must be LEFT or RIGHT");
196: }
197: *xred = NULL;
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: /*
202: PCApply_SVD - Applies the SVD preconditioner to a vector.
204: Input Parameters:
205: . pc - the preconditioner context
206: . x - input vector
208: Output Parameter:
209: . y - output vector
211: Application Interface Routine: PCApply()
212: */
213: static PetscErrorCode PCApply_SVD(PC pc, Vec x, Vec y)
214: {
215: PC_SVD *jac = (PC_SVD *)pc->data;
216: Vec work = jac->work, xred, yred;
218: PetscFunctionBegin;
219: PetscCall(PCSVDGetVec(pc, PC_RIGHT, READ, x, &xred));
220: PetscCall(PCSVDGetVec(pc, PC_LEFT, WRITE, y, &yred));
221: #if !defined(PETSC_USE_COMPLEX)
222: PetscCall(MatMultTranspose(jac->U, xred, work));
223: #else
224: PetscCall(MatMultHermitianTranspose(jac->U, xred, work));
225: #endif
226: PetscCall(VecPointwiseMult(work, work, jac->diag));
227: #if !defined(PETSC_USE_COMPLEX)
228: PetscCall(MatMultTranspose(jac->Vt, work, yred));
229: #else
230: PetscCall(MatMultHermitianTranspose(jac->Vt, work, yred));
231: #endif
232: PetscCall(PCSVDRestoreVec(pc, PC_RIGHT, READ, x, &xred));
233: PetscCall(PCSVDRestoreVec(pc, PC_LEFT, WRITE, y, &yred));
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: static PetscErrorCode PCMatApply_SVD(PC pc, Mat X, Mat Y)
238: {
239: PC_SVD *jac = (PC_SVD *)pc->data;
240: Mat W;
242: PetscFunctionBegin;
243: PetscCall(MatTransposeMatMult(jac->U, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &W));
244: PetscCall(MatDiagonalScale(W, jac->diag, NULL));
245: PetscCall(MatTransposeMatMult(jac->Vt, W, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y));
246: PetscCall(MatDestroy(&W));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: static PetscErrorCode PCApplyTranspose_SVD(PC pc, Vec x, Vec y)
251: {
252: PC_SVD *jac = (PC_SVD *)pc->data;
253: Vec work = jac->work, xred, yred;
255: PetscFunctionBegin;
256: PetscCall(PCSVDGetVec(pc, PC_LEFT, READ, x, &xred));
257: PetscCall(PCSVDGetVec(pc, PC_RIGHT, WRITE, y, &yred));
258: PetscCall(MatMult(jac->Vt, xred, work));
259: PetscCall(VecPointwiseMult(work, work, jac->diag));
260: PetscCall(MatMult(jac->U, work, yred));
261: PetscCall(PCSVDRestoreVec(pc, PC_LEFT, READ, x, &xred));
262: PetscCall(PCSVDRestoreVec(pc, PC_RIGHT, WRITE, y, &yred));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: static PetscErrorCode PCReset_SVD(PC pc)
267: {
268: PC_SVD *jac = (PC_SVD *)pc->data;
270: PetscFunctionBegin;
271: PetscCall(MatDestroy(&jac->A));
272: PetscCall(MatDestroy(&jac->U));
273: PetscCall(MatDestroy(&jac->Vt));
274: PetscCall(VecDestroy(&jac->diag));
275: PetscCall(VecDestroy(&jac->work));
276: PetscCall(VecScatterDestroy(&jac->right2red));
277: PetscCall(VecScatterDestroy(&jac->left2red));
278: PetscCall(VecDestroy(&jac->rightred));
279: PetscCall(VecDestroy(&jac->leftred));
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: /*
284: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
285: that was created with PCCreate_SVD().
287: Input Parameter:
288: . pc - the preconditioner context
290: Application Interface Routine: PCDestroy()
291: */
292: static PetscErrorCode PCDestroy_SVD(PC pc)
293: {
294: PC_SVD *jac = (PC_SVD *)pc->data;
296: PetscFunctionBegin;
297: PetscCall(PCReset_SVD(pc));
298: PetscCall(PetscViewerDestroy(&jac->monitor));
299: PetscCall(PetscFree(pc->data));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode PCSetFromOptions_SVD(PC pc, PetscOptionItems *PetscOptionsObject)
304: {
305: PC_SVD *jac = (PC_SVD *)pc->data;
306: PetscBool flg;
308: PetscFunctionBegin;
309: PetscOptionsHeadBegin(PetscOptionsObject, "SVD options");
310: PetscCall(PetscOptionsReal("-pc_svd_zero_sing", "Singular values smaller than this treated as zero", "None", jac->zerosing, &jac->zerosing, NULL));
311: PetscCall(PetscOptionsInt("-pc_svd_ess_rank", "Essential rank of operator (0 to use entire operator)", "None", jac->essrank, &jac->essrank, NULL));
312: PetscCall(PetscOptionsViewer("-pc_svd_monitor", "Monitor the conditioning, and extremal singular values", "None", &jac->monitor, &jac->monitorformat, &flg));
313: PetscOptionsHeadEnd();
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: static PetscErrorCode PCView_SVD(PC pc, PetscViewer viewer)
318: {
319: PC_SVD *svd = (PC_SVD *)pc->data;
320: PetscBool iascii;
322: PetscFunctionBegin;
323: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
324: if (iascii) {
325: PetscCall(PetscViewerASCIIPrintf(viewer, " All singular values smaller than %g treated as zero\n", (double)svd->zerosing));
326: PetscCall(PetscViewerASCIIPrintf(viewer, " Provided essential rank of the matrix %" PetscInt_FMT " (all other eigenvalues are zeroed)\n", svd->essrank));
327: }
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: /*
332: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
333: and sets this as the private data within the generic preconditioning
334: context, PC, that was created within PCCreate().
336: Input Parameter:
337: . pc - the preconditioner context
339: Application Interface Routine: PCCreate()
340: */
342: /*MC
343: PCSVD - Use pseudo inverse defined by SVD of operator
345: Level: advanced
347: Options Database Keys:
348: + -pc_svd_zero_sing <rtol> - Singular values smaller than this are treated as zero
349: - -pc_svd_monitor - Print information on the extreme singular values of the operator
351: Developer Note:
352: This implementation automatically creates a redundant copy of the
353: matrix on each process and uses a sequential SVD solve. Why does it do this instead
354: of using the composable `PCREDUNDANT` object?
356: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCREDUNDANT`
357: M*/
359: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
360: {
361: PC_SVD *jac;
362: PetscMPIInt size = 0;
364: PetscFunctionBegin;
365: /*
366: Creates the private data structure for this preconditioner and
367: attach it to the PC object.
368: */
369: PetscCall(PetscNew(&jac));
370: jac->zerosing = 1.e-12;
371: pc->data = (void *)jac;
373: /*
374: Set the pointers for the functions that are provided above.
375: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
376: are called, they will automatically call these functions. Note we
377: choose not to provide a couple of these functions since they are
378: not needed.
379: */
381: #if defined(PETSC_HAVE_COMPLEX)
382: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
383: #endif
384: if (size == 1) pc->ops->matapply = PCMatApply_SVD;
385: pc->ops->apply = PCApply_SVD;
386: pc->ops->applytranspose = PCApplyTranspose_SVD;
387: pc->ops->setup = PCSetUp_SVD;
388: pc->ops->reset = PCReset_SVD;
389: pc->ops->destroy = PCDestroy_SVD;
390: pc->ops->setfromoptions = PCSetFromOptions_SVD;
391: pc->ops->view = PCView_SVD;
392: pc->ops->applyrichardson = NULL;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }