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:   MatDestroy(&jac->A);
 48:   MPI_Comm_size(((PetscObject)pc->pmat)->comm, &size);
 49:   if (size > 1) {
 50:     Mat redmat;

 52:     MatCreateRedundantMatrix(pc->pmat, size, PETSC_COMM_SELF, MAT_INITIAL_MATRIX, &redmat);
 53:     MatConvert(redmat, MATSEQDENSE, MAT_INITIAL_MATRIX, &jac->A);
 54:     MatDestroy(&redmat);
 55:   } else {
 56:     MatConvert(pc->pmat, MATSEQDENSE, MAT_INITIAL_MATRIX, &jac->A);
 57:   }
 58:   if (!jac->diag) { /* assume square matrices */
 59:     MatCreateVecs(jac->A, &jac->diag, &jac->work);
 60:   }
 61:   if (!jac->U) {
 62:     MatDuplicate(jac->A, MAT_DO_NOT_COPY_VALUES, &jac->U);
 63:     MatDuplicate(jac->A, MAT_DO_NOT_COPY_VALUES, &jac->Vt);
 64:   }
 65:   MatGetSize(jac->A, &n, NULL);
 66:   if (!n) {
 67:     PetscInfo(pc, "Matrix has zero rows, skipping svd\n");
 68:     return 0;
 69:   }
 70:   PetscBLASIntCast(n, &nb);
 71:   lwork = 5 * nb;
 72:   PetscMalloc1(lwork, &work);
 73:   MatDenseGetArray(jac->A, &a);
 74:   MatDenseGetArray(jac->U, &u);
 75:   MatDenseGetArray(jac->Vt, &v);
 76:   VecGetArray(jac->diag, &d);
 77: #if !defined(PETSC_USE_COMPLEX)
 78:   {
 79:     PetscBLASInt lierr;
 80:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 81:     PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &nb, &nb, a, &nb, d, u, &nb, v, &nb, work, &lwork, &lierr));
 83:     PetscFPTrapPop();
 84:   }
 85: #else
 86:   {
 87:     PetscBLASInt lierr;
 88:     PetscReal   *rwork, *dd;
 89:     PetscMalloc1(5 * nb, &rwork);
 90:     PetscMalloc1(nb, &dd);
 91:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 92:     PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &nb, &nb, a, &nb, dd, u, &nb, v, &nb, work, &lwork, rwork, &lierr));
 94:     PetscFree(rwork);
 95:     for (i = 0; i < n; i++) d[i] = dd[i];
 96:     PetscFree(dd);
 97:     PetscFPTrapPop();
 98:   }
 99: #endif
100:   MatDenseRestoreArray(jac->A, &a);
101:   MatDenseRestoreArray(jac->U, &u);
102:   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:     PetscViewerASCIIAddTab(jac->monitor, ((PetscObject)pc)->tablevel);
108:     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:       PetscViewerASCIIPrintf(jac->monitor, "    SVD: singular values:\n");
111:       for (i = 0; i < n; i++) {
112:         if (i % 5 == 0) {
113:           if (i != 0) PetscViewerASCIIPrintf(jac->monitor, "\n");
114:           PetscViewerASCIIPrintf(jac->monitor, "        ");
115:         }
116:         PetscViewerASCIIPrintf(jac->monitor, " %14.12e", (double)PetscRealPart(d[i]));
117:       }
118:       PetscViewerASCIIPrintf(jac->monitor, "\n");
119:     } else { /* print 5 smallest and 5 largest */
120:       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:       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:     PetscViewerASCIISubtractTab(jac->monitor, ((PetscObject)pc)->tablevel);
124:   }
125:   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:   PetscInfo(pc, "Number of zero or nearly singular values %" PetscInt_FMT "\n", jac->nzero);
131:   VecRestoreArray(jac->diag, &d);
132:   PetscFree(work);
133:   return 0;
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:   MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
142:   *xred = NULL;
143:   switch (side) {
144:   case PC_LEFT:
145:     if (size == 1) *xred = x;
146:     else {
147:       if (!jac->left2red) VecScatterCreateToAll(x, &jac->left2red, &jac->leftred);
148:       if (amode & READ) {
149:         VecScatterBegin(jac->left2red, x, jac->leftred, INSERT_VALUES, SCATTER_FORWARD);
150:         VecScatterEnd(jac->left2red, x, jac->leftred, INSERT_VALUES, SCATTER_FORWARD);
151:       }
152:       *xred = jac->leftred;
153:     }
154:     break;
155:   case PC_RIGHT:
156:     if (size == 1) *xred = x;
157:     else {
158:       if (!jac->right2red) VecScatterCreateToAll(x, &jac->right2red, &jac->rightred);
159:       if (amode & READ) {
160:         VecScatterBegin(jac->right2red, x, jac->rightred, INSERT_VALUES, SCATTER_FORWARD);
161:         VecScatterEnd(jac->right2red, x, jac->rightred, INSERT_VALUES, SCATTER_FORWARD);
162:       }
163:       *xred = jac->rightred;
164:     }
165:     break;
166:   default:
167:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Side must be LEFT or RIGHT");
168:   }
169:   return 0;
170: }

172: static PetscErrorCode PCSVDRestoreVec(PC pc, PCSide side, AccessMode amode, Vec x, Vec *xred)
173: {
174:   PC_SVD     *jac = (PC_SVD *)pc->data;
175:   PetscMPIInt size;

177:   MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
178:   switch (side) {
179:   case PC_LEFT:
180:     if (size != 1 && amode & WRITE) {
181:       VecScatterBegin(jac->left2red, jac->leftred, x, INSERT_VALUES, SCATTER_REVERSE);
182:       VecScatterEnd(jac->left2red, jac->leftred, x, INSERT_VALUES, SCATTER_REVERSE);
183:     }
184:     break;
185:   case PC_RIGHT:
186:     if (size != 1 && amode & WRITE) {
187:       VecScatterBegin(jac->right2red, jac->rightred, x, INSERT_VALUES, SCATTER_REVERSE);
188:       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:   return 0;
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:   PCSVDGetVec(pc, PC_RIGHT, READ, x, &xred);
216:   PCSVDGetVec(pc, PC_LEFT, WRITE, y, &yred);
217: #if !defined(PETSC_USE_COMPLEX)
218:   MatMultTranspose(jac->U, xred, work);
219: #else
220:   MatMultHermitianTranspose(jac->U, xred, work);
221: #endif
222:   VecPointwiseMult(work, work, jac->diag);
223: #if !defined(PETSC_USE_COMPLEX)
224:   MatMultTranspose(jac->Vt, work, yred);
225: #else
226:   MatMultHermitianTranspose(jac->Vt, work, yred);
227: #endif
228:   PCSVDRestoreVec(pc, PC_RIGHT, READ, x, &xred);
229:   PCSVDRestoreVec(pc, PC_LEFT, WRITE, y, &yred);
230:   return 0;
231: }

233: static PetscErrorCode PCMatApply_SVD(PC pc, Mat X, Mat Y)
234: {
235:   PC_SVD *jac = (PC_SVD *)pc->data;
236:   Mat     W;

238:   MatTransposeMatMult(jac->U, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &W);
239:   MatDiagonalScale(W, jac->diag, NULL);
240:   MatTransposeMatMult(jac->Vt, W, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y);
241:   MatDestroy(&W);
242:   return 0;
243: }

245: static PetscErrorCode PCApplyTranspose_SVD(PC pc, Vec x, Vec y)
246: {
247:   PC_SVD *jac  = (PC_SVD *)pc->data;
248:   Vec     work = jac->work, xred, yred;

250:   PCSVDGetVec(pc, PC_LEFT, READ, x, &xred);
251:   PCSVDGetVec(pc, PC_RIGHT, WRITE, y, &yred);
252:   MatMult(jac->Vt, xred, work);
253:   VecPointwiseMult(work, work, jac->diag);
254:   MatMult(jac->U, work, yred);
255:   PCSVDRestoreVec(pc, PC_LEFT, READ, x, &xred);
256:   PCSVDRestoreVec(pc, PC_RIGHT, WRITE, y, &yred);
257:   return 0;
258: }

260: static PetscErrorCode PCReset_SVD(PC pc)
261: {
262:   PC_SVD *jac = (PC_SVD *)pc->data;

264:   MatDestroy(&jac->A);
265:   MatDestroy(&jac->U);
266:   MatDestroy(&jac->Vt);
267:   VecDestroy(&jac->diag);
268:   VecDestroy(&jac->work);
269:   VecScatterDestroy(&jac->right2red);
270:   VecScatterDestroy(&jac->left2red);
271:   VecDestroy(&jac->rightred);
272:   VecDestroy(&jac->leftred);
273:   return 0;
274: }

276: /*
277:    PCDestroy_SVD - Destroys the private context for the SVD preconditioner
278:    that was created with PCCreate_SVD().

280:    Input Parameter:
281: .  pc - the preconditioner context

283:    Application Interface Routine: PCDestroy()
284: */
285: static PetscErrorCode PCDestroy_SVD(PC pc)
286: {
287:   PC_SVD *jac = (PC_SVD *)pc->data;

289:   PCReset_SVD(pc);
290:   PetscViewerDestroy(&jac->monitor);
291:   PetscFree(pc->data);
292:   return 0;
293: }

295: static PetscErrorCode PCSetFromOptions_SVD(PC pc, PetscOptionItems *PetscOptionsObject)
296: {
297:   PC_SVD   *jac = (PC_SVD *)pc->data;
298:   PetscBool flg;

300:   PetscOptionsHeadBegin(PetscOptionsObject, "SVD options");
301:   PetscOptionsReal("-pc_svd_zero_sing", "Singular values smaller than this treated as zero", "None", jac->zerosing, &jac->zerosing, NULL);
302:   PetscOptionsInt("-pc_svd_ess_rank", "Essential rank of operator (0 to use entire operator)", "None", jac->essrank, &jac->essrank, NULL);
303:   PetscOptionsViewer("-pc_svd_monitor", "Monitor the conditioning, and extremal singular values", "None", &jac->monitor, &jac->monitorformat, &flg);
304:   PetscOptionsHeadEnd();
305:   return 0;
306: }

308: static PetscErrorCode PCView_SVD(PC pc, PetscViewer viewer)
309: {
310:   PC_SVD   *svd = (PC_SVD *)pc->data;
311:   PetscBool iascii;

313:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
314:   if (iascii) {
315:     PetscViewerASCIIPrintf(viewer, "  All singular values smaller than %g treated as zero\n", (double)svd->zerosing);
316:     PetscViewerASCIIPrintf(viewer, "  Provided essential rank of the matrix %" PetscInt_FMT " (all other eigenvalues are zeroed)\n", svd->essrank);
317:   }
318:   return 0;
319: }

321: /*
322:    PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
323:    and sets this as the private data within the generic preconditioning
324:    context, PC, that was created within PCCreate().

326:    Input Parameter:
327: .  pc - the preconditioner context

329:    Application Interface Routine: PCCreate()
330: */

332: /*MC
333:      PCSVD - Use pseudo inverse defined by SVD of operator

335:    Level: advanced

337:   Options Database Keys:
338: +  -pc_svd_zero_sing <rtol> - Singular values smaller than this are treated as zero
339: -  -pc_svd_monitor - Print information on the extreme singular values of the operator

341:   Developer Note:
342:   This implementation automatically creates a redundant copy of the
343:    matrix on each process and uses a sequential SVD solve. Why does it do this instead
344:    of using the composable `PCREDUNDANT` object?

346: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCREDUNDANT`
347: M*/

349: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
350: {
351:   PC_SVD     *jac;
352:   PetscMPIInt size = 0;

354:   /*
355:      Creates the private data structure for this preconditioner and
356:      attach it to the PC object.
357:   */
358:   PetscNew(&jac);
359:   jac->zerosing = 1.e-12;
360:   pc->data      = (void *)jac;

362:   /*
363:       Set the pointers for the functions that are provided above.
364:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
365:       are called, they will automatically call these functions.  Note we
366:       choose not to provide a couple of these functions since they are
367:       not needed.
368:   */

370: #if defined(PETSC_HAVE_COMPLEX)
371:   MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
372: #endif
373:   if (size == 1) pc->ops->matapply = PCMatApply_SVD;
374:   pc->ops->apply           = PCApply_SVD;
375:   pc->ops->applytranspose  = PCApplyTranspose_SVD;
376:   pc->ops->setup           = PCSetUp_SVD;
377:   pc->ops->reset           = PCReset_SVD;
378:   pc->ops->destroy         = PCDestroy_SVD;
379:   pc->ops->setfromoptions  = PCSetFromOptions_SVD;
380:   pc->ops->view            = PCView_SVD;
381:   pc->ops->applyrichardson = NULL;
382:   return 0;
383: }