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