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