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