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: }