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