1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
2 #include <petscblaslapack.h>
3
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;
18
19 typedef enum {
20 READ = 1,
21 WRITE = 2,
22 READ_WRITE = 3
23 } AccessMode;
24
25 /*
26 PCSetUp_SVD - Prepares for the use of the SVD preconditioner
27 by setting data structures and options.
28
29 Input Parameter:
30 . pc - the preconditioner context
31
32 Application Interface Routine: PCSetUp()
33
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 */
PCSetUp_SVD(PC pc)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;
45
46 PetscFunctionBegin;
47 PetscCall(MatDestroy(&jac->A));
48 PetscCallMPI(MPI_Comm_size(((PetscObject)pc->pmat)->comm, &size));
49 if (size > 1) {
50 Mat redmat;
51
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);
67
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 }
133
PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec * xred)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;
138
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 }
170
PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec * xred)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;
175
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 }
197
198 /*
199 PCApply_SVD - Applies the SVD preconditioner to a vector.
200
201 Input Parameters:
202 . pc - the preconditioner context
203 . x - input vector
204
205 Output Parameter:
206 . y - output vector
207
208 Application Interface Routine: PCApply()
209 */
PCApply_SVD(PC pc,Vec x,Vec y)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;
214
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 }
233
PCMatApply_SVD(PC pc,Mat X,Mat Y)234 static PetscErrorCode PCMatApply_SVD(PC pc, Mat X, Mat Y)
235 {
236 PC_SVD *jac = (PC_SVD *)pc->data;
237 Mat W;
238
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 }
246
PCApplyTranspose_SVD(PC pc,Vec x,Vec y)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;
251
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 }
262
PCReset_SVD(PC pc)263 static PetscErrorCode PCReset_SVD(PC pc)
264 {
265 PC_SVD *jac = (PC_SVD *)pc->data;
266
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 }
279
280 /*
281 PCDestroy_SVD - Destroys the private context for the SVD preconditioner
282 that was created with PCCreate_SVD().
283
284 Input Parameter:
285 . pc - the preconditioner context
286
287 Application Interface Routine: PCDestroy()
288 */
PCDestroy_SVD(PC pc)289 static PetscErrorCode PCDestroy_SVD(PC pc)
290 {
291 PC_SVD *jac = (PC_SVD *)pc->data;
292
293 PetscFunctionBegin;
294 PetscCall(PCReset_SVD(pc));
295 PetscCall(PetscViewerDestroy(&jac->monitor));
296 PetscCall(PetscFree(pc->data));
297 PetscFunctionReturn(PETSC_SUCCESS);
298 }
299
PCSetFromOptions_SVD(PC pc,PetscOptionItems PetscOptionsObject)300 static PetscErrorCode PCSetFromOptions_SVD(PC pc, PetscOptionItems PetscOptionsObject)
301 {
302 PC_SVD *jac = (PC_SVD *)pc->data;
303 PetscBool flg;
304
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 }
313
PCView_SVD(PC pc,PetscViewer viewer)314 static PetscErrorCode PCView_SVD(PC pc, PetscViewer viewer)
315 {
316 PC_SVD *svd = (PC_SVD *)pc->data;
317 PetscBool isascii;
318
319 PetscFunctionBegin;
320 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
321 if (isascii) {
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 }
327
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().
332
333 Input Parameter:
334 . pc - the preconditioner context
335
336 Application Interface Routine: PCCreate()
337 */
338
339 /*MC
340 PCSVD - Use pseudo inverse defined by SVD of operator
341
342 Level: advanced
343
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
347
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?
352
353 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCREDUNDANT`
354 M*/
355
PCCreate_SVD(PC pc)356 PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
357 {
358 PC_SVD *jac;
359 PetscMPIInt size = 0;
360
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;
369
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 */
377
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 }
392