Actual source code: svd.c

petsc-master 2016-12-06
Report Typos and Errors

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

 19: typedef enum {READ=1, WRITE=2, READ_WRITE=3} AccessMode;

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

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

 29:    Application Interface Routine: PCSetUp()

 31:    Notes:
 32:    The interface routine PCSetUp() is not usually called directly by
 33:    the user, but instead is called by PCApply() if necessary.
 34: */
 37: static PetscErrorCode PCSetUp_SVD(PC pc)
 38: {
 39: #if defined(PETSC_MISSING_LAPACK_GESVD)
 40:   SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable\nNot able to provide singular value estimates.");
 41: #else
 42:   PC_SVD         *jac = (PC_SVD*)pc->data;
 44:   PetscScalar    *a,*u,*v,*d,*work;
 45:   PetscBLASInt   nb,lwork;
 46:   PetscInt       i,n;
 47:   PetscMPIInt    size;

 50:   MatDestroy(&jac->A);
 51:   MPI_Comm_size(((PetscObject)pc->pmat)->comm,&size);
 52:   if (size > 1) {
 53:     Mat          redmat;
 54:     MatCreateRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&redmat);
 55:     MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
 56:     MatDestroy(&redmat);
 57:   } else {
 58:     MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
 59:   }
 60:   if (!jac->diag) {    /* assume square matrices */
 61:     MatCreateVecs(jac->A,&jac->diag,&jac->work);
 62:   }
 63:   if (!jac->U) {
 64:     MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
 65:     MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);
 66:   }
 67:   MatGetSize(jac->A,&n,NULL);
 68:   if (!n) {
 69:     PetscInfo(pc,"Matrix has zero rows, skipping svd\n");
 70:     return(0);
 71:   }
 72:   PetscBLASIntCast(n,&nb);
 73:   lwork = 5*nb;
 74:   PetscMalloc1(lwork,&work);
 75:   MatDenseGetArray(jac->A,&a);
 76:   MatDenseGetArray(jac->U,&u);
 77:   MatDenseGetArray(jac->Vt,&v);
 78:   VecGetArray(jac->diag,&d);
 79: #if !defined(PETSC_USE_COMPLEX)
 80:   {
 81:     PetscBLASInt lierr;
 82:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 83:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
 84:     if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
 85:     PetscFPTrapPop();
 86:   }
 87: #else
 88:   {
 89:     PetscBLASInt lierr;
 90:     PetscReal    *rwork,*dd;
 91:     PetscMalloc1(5*nb,&rwork);
 92:     PetscMalloc1(nb,&dd);
 93:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 94:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,dd,u,&nb,v,&nb,work,&lwork,rwork,&lierr));
 95:     if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
 96:     PetscFree(rwork);
 97:     for (i=0; i<n; i++) d[i] = dd[i];
 98:     PetscFree(dd);
 99:     PetscFPTrapPop();
100:   }
101: #endif
102:   MatDenseRestoreArray(jac->A,&a);
103:   MatDenseRestoreArray(jac->U,&u);
104:   MatDenseRestoreArray(jac->Vt,&v);
105:   for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
106:   jac->nzero = n-1-i;
107:   if (jac->monitor) {
108:     PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
109:     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);
110:     if (n >= 10) {              /* print 5 smallest and 5 largest */
111:       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]));
112:       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]));
113:     } else {                    /* print all singular values */
114:       char     buf[256],*p;
115:       size_t   left = sizeof(buf),used;
116:       PetscInt thisline;
117:       for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
118:         PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));
119:         left -= used;
120:         p    += used;
121:         if (thisline > 4 || i==0) {
122:           PetscViewerASCIIPrintf(jac->monitor,"    SVD: singular values:%s\n",buf);
123:           p        = buf;
124:           thisline = 0;
125:         }
126:       }
127:     }
128:     PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
129:   }
130:   PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));
131:   for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
132:   for (; i<n; i++) d[i] = 0.0;
133:   if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
134:   PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
135:   VecRestoreArray(jac->diag,&d);
136: #if defined(foo)
137:   {
138:     PetscViewer viewer;
139:     PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);
140:     MatView(jac->A,viewer);
141:     MatView(jac->U,viewer);
142:     MatView(jac->Vt,viewer);
143:     VecView(jac->diag,viewer);
144:     PetscViewerDestroy(viewer);
145:   }
146: #endif
147:   PetscFree(work);
148:   return(0);
149: #endif
150: }

154: static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
155: {
156:   PC_SVD         *jac = (PC_SVD*)pc->data;
158:   PetscMPIInt    size;

161:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
162:   *xred = NULL;
163:   switch (side) {
164:   case PC_LEFT:
165:     if (size == 1) *xred = x;
166:     else {
167:       if (!jac->left2red) {VecScatterCreateToAll(x,&jac->left2red,&jac->leftred);}
168:       if (amode & READ) {
169:         VecScatterBegin(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
170:         VecScatterEnd(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);
171:       }
172:       *xred = jac->leftred;
173:     }
174:     break;
175:   case PC_RIGHT:
176:     if (size == 1) *xred = x;
177:     else {
178:       if (!jac->right2red) {VecScatterCreateToAll(x,&jac->right2red,&jac->rightred);}
179:       if (amode & READ) {
180:         VecScatterBegin(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
181:         VecScatterEnd(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);
182:       }
183:       *xred = jac->rightred;
184:     }
185:     break;
186:   default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
187:   }
188:   return(0);
189: }

193: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
194: {
195:   PC_SVD         *jac = (PC_SVD*)pc->data;
197:   PetscMPIInt    size;

200:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
201:   switch (side) {
202:   case PC_LEFT:
203:     if (size != 1 && amode & WRITE) {
204:       VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
205:       VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
206:     }
207:     break;
208:   case PC_RIGHT:
209:     if (size != 1 && amode & WRITE) {
210:       VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
211:       VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
212:     }
213:     break;
214:   default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
215:   }
216:   *xred = NULL;
217:   return(0);
218: }

220: /* -------------------------------------------------------------------------- */
221: /*
222:    PCApply_SVD - Applies the SVD preconditioner to a vector.

224:    Input Parameters:
225: .  pc - the preconditioner context
226: .  x - input vector

228:    Output Parameter:
229: .  y - output vector

231:    Application Interface Routine: PCApply()
232:  */
235: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
236: {
237:   PC_SVD         *jac = (PC_SVD*)pc->data;
238:   Vec            work = jac->work,xred,yred;

242:   PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
243:   PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
244: #if !defined(PETSC_USE_COMPLEX)
245:   MatMultTranspose(jac->U,xred,work);
246: #else
247:   MatMultHermitianTranspose(jac->U,xred,work);
248: #endif
249:   VecPointwiseMult(work,work,jac->diag);
250: #if !defined(PETSC_USE_COMPLEX)
251:   MatMultTranspose(jac->Vt,work,yred);
252: #else
253:   MatMultHermitianTranspose(jac->Vt,work,yred);
254: #endif
255:   PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
256:   PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
257:   return(0);
258: }

262: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
263: {
264:   PC_SVD         *jac = (PC_SVD*)pc->data;
265:   Vec            work = jac->work,xred,yred;

269:   PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);
270:   PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);
271: #if !defined(PETSC_USE_COMPLEX)
272:   MatMultTranspose(jac->Vt,work,yred);
273: #else
274:   MatMultHermitianTranspose(jac->Vt,work,yred);
275: #endif
276:   VecPointwiseMult(work,work,jac->diag);
277: #if !defined(PETSC_USE_COMPLEX)
278:   MatMultTranspose(jac->U,xred,work);
279: #else
280:   MatMultHermitianTranspose(jac->U,xred,work);
281: #endif
282:   PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);
283:   PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);
284:   return(0);
285: }

289: static PetscErrorCode PCReset_SVD(PC pc)
290: {
291:   PC_SVD         *jac = (PC_SVD*)pc->data;

295:   MatDestroy(&jac->A);
296:   MatDestroy(&jac->U);
297:   MatDestroy(&jac->Vt);
298:   VecDestroy(&jac->diag);
299:   VecDestroy(&jac->work);
300:   VecScatterDestroy(&jac->right2red);
301:   VecScatterDestroy(&jac->left2red);
302:   VecDestroy(&jac->rightred);
303:   VecDestroy(&jac->leftred);
304:   return(0);
305: }

307: /* -------------------------------------------------------------------------- */
308: /*
309:    PCDestroy_SVD - Destroys the private context for the SVD preconditioner
310:    that was created with PCCreate_SVD().

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

315:    Application Interface Routine: PCDestroy()
316: */
319: static PetscErrorCode PCDestroy_SVD(PC pc)
320: {
321:   PC_SVD         *jac = (PC_SVD*)pc->data;

325:   PCReset_SVD(pc);
326:   PetscViewerDestroy(&jac->monitor);
327:   PetscFree(pc->data);
328:   return(0);
329: }

333: static PetscErrorCode PCSetFromOptions_SVD(PetscOptionItems *PetscOptionsObject,PC pc)
334: {
336:   PC_SVD         *jac = (PC_SVD*)pc->data;
337:   PetscBool      flg,set;

340:   PetscOptionsHead(PetscOptionsObject,"SVD options");
341:   PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
342:   PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
343:   PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
344:   if (set) {                    /* Should make PCSVDSetMonitor() */
345:     if (flg && !jac->monitor) {
346:       PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);
347:     } else if (!flg) {
348:       PetscViewerDestroy(&jac->monitor);
349:     }
350:   }
351:   PetscOptionsTail();
352:   return(0);
353: }

357: static PetscErrorCode PCView_SVD(PC pc,PetscViewer viewer)
358: {
359:   PC_SVD         *svd = (PC_SVD*)pc->data;
361:   PetscBool      iascii;

364:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
365:   if (iascii) {
366:     PetscViewerASCIIPrintf(viewer,"  SVD: All singular values smaller than %g treated as zero\n",(double)svd->zerosing);
367:     PetscViewerASCIIPrintf(viewer,"  SVD: Provided essential rank of the matrix %D (all other eigenvalues are zeroed)\n",svd->essrank);
368:   }
369:   return(0);
370: }
371: /* -------------------------------------------------------------------------- */
372: /*
373:    PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
374:    and sets this as the private data within the generic preconditioning
375:    context, PC, that was created within PCCreate().

377:    Input Parameter:
378: .  pc - the preconditioner context

380:    Application Interface Routine: PCCreate()
381: */

383: /*MC
384:      PCSVD - Use pseudo inverse defined by SVD of operator

386:    Level: advanced

388:   Concepts: SVD

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

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

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

403: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
404: {
405:   PC_SVD         *jac;

409:   /*
410:      Creates the private data structure for this preconditioner and
411:      attach it to the PC object.
412:   */
413:   PetscNewLog(pc,&jac);
414:   jac->zerosing = 1.e-12;
415:   pc->data      = (void*)jac;

417:   /*
418:       Set the pointers for the functions that are provided above.
419:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
420:       are called, they will automatically call these functions.  Note we
421:       choose not to provide a couple of these functions since they are
422:       not needed.
423:   */
424:   pc->ops->apply           = PCApply_SVD;
425:   pc->ops->applytranspose  = PCApplyTranspose_SVD;
426:   pc->ops->setup           = PCSetUp_SVD;
427:   pc->ops->reset           = PCReset_SVD;
428:   pc->ops->destroy         = PCDestroy_SVD;
429:   pc->ops->setfromoptions  = PCSetFromOptions_SVD;
430:   pc->ops->view            = PCView_SVD;
431:   pc->ops->applyrichardson = 0;
432:   return(0);
433: }