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 {READ=1, WRITE=2, READ_WRITE=3} AccessMode;

 22: /* -------------------------------------------------------------------------- */
 23: /*
 24:    PCSetUp_SVD - Prepares for the use of the SVD preconditioner
 25:                     by setting data structures and options.

 27:    Input Parameter:
 28: .  pc - the preconditioner context

 30:    Application Interface Routine: PCSetUp()

 32:    Notes:
 33:    The interface routine PCSetUp() is not usually called directly by
 34:    the user, but instead is called by PCApply() if necessary.
 35: */
 36: static PetscErrorCode PCSetUp_SVD(PC pc)
 37: {
 38:   PC_SVD         *jac = (PC_SVD*)pc->data;
 39:   PetscScalar    *a,*u,*v,*d,*work;
 40:   PetscBLASInt   nb,lwork;
 41:   PetscInt       i,n;
 42:   PetscMPIInt    size;

 44:   MatDestroy(&jac->A);
 45:   MPI_Comm_size(((PetscObject)pc->pmat)->comm,&size);
 46:   if (size > 1) {
 47:     Mat redmat;

 49:     MatCreateRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&redmat);
 50:     MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
 51:     MatDestroy(&redmat);
 52:   } else {
 53:     MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
 54:   }
 55:   if (!jac->diag) {    /* assume square matrices */
 56:     MatCreateVecs(jac->A,&jac->diag,&jac->work);
 57:   }
 58:   if (!jac->U) {
 59:     MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
 60:     MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);
 61:   }
 62:   MatGetSize(jac->A,&n,NULL);
 63:   if (!n) {
 64:     PetscInfo(pc,"Matrix has zero rows, skipping svd\n");
 65:     return 0;
 66:   }
 67:   PetscBLASIntCast(n,&nb);
 68:   lwork = 5*nb;
 69:   PetscMalloc1(lwork,&work);
 70:   MatDenseGetArray(jac->A,&a);
 71:   MatDenseGetArray(jac->U,&u);
 72:   MatDenseGetArray(jac->Vt,&v);
 73:   VecGetArray(jac->diag,&d);
 74: #if !defined(PETSC_USE_COMPLEX)
 75:   {
 76:     PetscBLASInt lierr;
 77:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 78:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
 80:     PetscFPTrapPop();
 81:   }
 82: #else
 83:   {
 84:     PetscBLASInt lierr;
 85:     PetscReal    *rwork,*dd;
 86:     PetscMalloc1(5*nb,&rwork);
 87:     PetscMalloc1(nb,&dd);
 88:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 89:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,dd,u,&nb,v,&nb,work,&lwork,rwork,&lierr));
 91:     PetscFree(rwork);
 92:     for (i=0; i<n; i++) d[i] = dd[i];
 93:     PetscFree(dd);
 94:     PetscFPTrapPop();
 95:   }
 96: #endif
 97:   MatDenseRestoreArray(jac->A,&a);
 98:   MatDenseRestoreArray(jac->U,&u);
 99:   MatDenseRestoreArray(jac->Vt,&v);
100:   for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
101:   jac->nzero = n-1-i;
102:   if (jac->monitor) {
103:     PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
104:     PetscViewerASCIIPrintf(jac->monitor,"    SVD: condition number %14.12e, %D of %D singular values are (nearly) zero\n",(double)PetscRealPart(d[0]/d[n-1]),jac->nzero,n);
105:     if (n < 10 || jac->monitorformat == PETSC_VIEWER_ALL) {
106:       PetscViewerASCIIPrintf(jac->monitor,"    SVD: singular values:\n");
107:       for (i=0; i<n; i++) {
108:         if (i%5 == 0) {
109:             if (i != 0) {
110:               PetscViewerASCIIPrintf(jac->monitor,"\n");
111:             }
112:             PetscViewerASCIIPrintf(jac->monitor,"        ");
113:           }
114:         PetscViewerASCIIPrintf(jac->monitor," %14.12e",(double)PetscRealPart(d[i]));
115:       }
116:       PetscViewerASCIIPrintf(jac->monitor,"\n");
117:     } else {              /* print 5 smallest and 5 largest */
118:       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:       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:     PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
122:   }
123:   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) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
127:   PetscInfo(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
128:   VecRestoreArray(jac->diag,&d);
129:   PetscFree(work);
130:   return 0;
131: }

133: static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
134: {
135:   PC_SVD         *jac = (PC_SVD*)pc->data;
136:   PetscMPIInt    size;

138:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
139:   *xred = NULL;
140:   switch (side) {
141:   case PC_LEFT:
142:     if (size == 1) *xred = x;
143:     else {
144:       if (!jac->left2red) VecScatterCreateToAll(x,&jac->left2red,&jac->leftred);
145:       if (amode & READ) {
146:         VecScatterBegin(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
147:         VecScatterEnd(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
148:       }
149:       *xred = jac->leftred;
150:     }
151:     break;
152:   case PC_RIGHT:
153:     if (size == 1) *xred = x;
154:     else {
155:       if (!jac->right2red) VecScatterCreateToAll(x,&jac->right2red,&jac->rightred);
156:       if (amode & READ) {
157:         VecScatterBegin(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
158:         VecScatterEnd(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
159:       }
160:       *xred = jac->rightred;
161:     }
162:     break;
163:   default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
164:   }
165:   return 0;
166: }

168: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
169: {
170:   PC_SVD         *jac = (PC_SVD*)pc->data;
171:   PetscMPIInt    size;

173:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
174:   switch (side) {
175:   case PC_LEFT:
176:     if (size != 1 && amode & WRITE) {
177:       VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
178:       VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
179:     }
180:     break;
181:   case PC_RIGHT:
182:     if (size != 1 && amode & WRITE) {
183:       VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
184:       VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
185:     }
186:     break;
187:   default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
188:   }
189:   *xred = NULL;
190:   return 0;
191: }

193: /* -------------------------------------------------------------------------- */
194: /*
195:    PCApply_SVD - Applies the SVD preconditioner to a vector.

197:    Input Parameters:
198: .  pc - the preconditioner context
199: .  x - input vector

201:    Output Parameter:
202: .  y - output vector

204:    Application Interface Routine: PCApply()
205:  */
206: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
207: {
208:   PC_SVD         *jac = (PC_SVD*)pc->data;
209:   Vec            work = jac->work,xred,yred;

211:   PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
212:   PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
213: #if !defined(PETSC_USE_COMPLEX)
214:   MatMultTranspose(jac->U,xred,work);
215: #else
216:   MatMultHermitianTranspose(jac->U,xred,work);
217: #endif
218:   VecPointwiseMult(work,work,jac->diag);
219: #if !defined(PETSC_USE_COMPLEX)
220:   MatMultTranspose(jac->Vt,work,yred);
221: #else
222:   MatMultHermitianTranspose(jac->Vt,work,yred);
223: #endif
224:   PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
225:   PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
226:   return 0;
227: }

229: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
230: {
231:   PC_SVD         *jac = (PC_SVD*)pc->data;
232:   Vec            work = jac->work,xred,yred;

234:   PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);
235:   PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);
236:   MatMult(jac->Vt,xred,work);
237:   VecPointwiseMult(work,work,jac->diag);
238:   MatMult(jac->U,work,yred);
239:   PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
240:   PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
241:   return 0;
242: }

244: static PetscErrorCode PCReset_SVD(PC pc)
245: {
246:   PC_SVD         *jac = (PC_SVD*)pc->data;

248:   MatDestroy(&jac->A);
249:   MatDestroy(&jac->U);
250:   MatDestroy(&jac->Vt);
251:   VecDestroy(&jac->diag);
252:   VecDestroy(&jac->work);
253:   VecScatterDestroy(&jac->right2red);
254:   VecScatterDestroy(&jac->left2red);
255:   VecDestroy(&jac->rightred);
256:   VecDestroy(&jac->leftred);
257:   return 0;
258: }

260: /* -------------------------------------------------------------------------- */
261: /*
262:    PCDestroy_SVD - Destroys the private context for the SVD preconditioner
263:    that was created with PCCreate_SVD().

265:    Input Parameter:
266: .  pc - the preconditioner context

268:    Application Interface Routine: PCDestroy()
269: */
270: static PetscErrorCode PCDestroy_SVD(PC pc)
271: {
272:   PC_SVD         *jac = (PC_SVD*)pc->data;

274:   PCReset_SVD(pc);
275:   PetscViewerDestroy(&jac->monitor);
276:   PetscFree(pc->data);
277:   return 0;
278: }

280: static PetscErrorCode PCSetFromOptions_SVD(PetscOptionItems *PetscOptionsObject,PC pc)
281: {
282:   PC_SVD         *jac = (PC_SVD*)pc->data;
283:   PetscBool      flg;

285:   PetscOptionsHead(PetscOptionsObject,"SVD options");
286:   PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
287:   PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
288:   PetscOptionsViewer("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",&jac->monitor,&jac->monitorformat,&flg);
289:   PetscOptionsTail();
290:   return 0;
291: }

293: static PetscErrorCode PCView_SVD(PC pc,PetscViewer viewer)
294: {
295:   PC_SVD         *svd = (PC_SVD*)pc->data;
296:   PetscBool      iascii;

298:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
299:   if (iascii) {
300:     PetscViewerASCIIPrintf(viewer,"  All singular values smaller than %g treated as zero\n",(double)svd->zerosing);
301:     PetscViewerASCIIPrintf(viewer,"  Provided essential rank of the matrix %D (all other eigenvalues are zeroed)\n",svd->essrank);
302:   }
303:   return 0;
304: }
305: /* -------------------------------------------------------------------------- */
306: /*
307:    PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
308:    and sets this as the private data within the generic preconditioning
309:    context, PC, that was created within PCCreate().

311:    Input Parameter:
312: .  pc - the preconditioner context

314:    Application Interface Routine: PCCreate()
315: */

317: /*MC
318:      PCSVD - Use pseudo inverse defined by SVD of operator

320:    Level: advanced

322:   Options Database:
323: +  -pc_svd_zero_sing <rtol> - Singular values smaller than this are treated as zero
324: -  -pc_svd_monitor - Print information on the extreme singular values of the operator

326:   Developer Note:
327:   This implementation automatically creates a redundant copy of the
328:    matrix on each process and uses a sequential SVD solve. Why does it do this instead
329:    of using the composable PCREDUNDANT object?

331: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC
332: M*/

334: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
335: {
336:   PC_SVD         *jac;

338:   /*
339:      Creates the private data structure for this preconditioner and
340:      attach it to the PC object.
341:   */
342:   PetscNewLog(pc,&jac);
343:   jac->zerosing = 1.e-12;
344:   pc->data      = (void*)jac;

346:   /*
347:       Set the pointers for the functions that are provided above.
348:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
349:       are called, they will automatically call these functions.  Note we
350:       choose not to provide a couple of these functions since they are
351:       not needed.
352:   */
353:   pc->ops->apply           = PCApply_SVD;
354:   pc->ops->applytranspose  = PCApplyTranspose_SVD;
355:   pc->ops->setup           = PCSetUp_SVD;
356:   pc->ops->reset           = PCReset_SVD;
357:   pc->ops->destroy         = PCDestroy_SVD;
358:   pc->ops->setfromoptions  = PCSetFromOptions_SVD;
359:   pc->ops->view            = PCView_SVD;
360:   pc->ops->applyrichardson = NULL;
361:   return 0;
362: }