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