Actual source code: svd.c

petsc-3.4.4 2014-03-13
  2: #include <petsc-private/pcimpl.h>   /*I "petscpc.h" I*/
  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:     PetscInt M;
 55:     MatGetSize(pc->pmat,&M,NULL);
 56:     MatGetRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,M,MAT_INITIAL_MATRIX,&redmat);
 57:     MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
 58:     MatDestroy(&redmat);
 59:   } else {
 60:     MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
 61:   }
 62:   if (!jac->diag) {    /* assume square matrices */
 63:     MatGetVecs(jac->A,&jac->diag,&jac->work);
 64:   }
 65:   if (!jac->U) {
 66:     MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
 67:     MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);
 68:   }
 69:   MatGetSize(pc->pmat,&n,NULL);
 70:   PetscBLASIntCast(n,&nb);
 71:   lwork = 5*nb;
 72:   PetscMalloc(lwork*sizeof(PetscScalar),&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:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
 82:     if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
 83:     PetscFPTrapPop();
 84:   }
 85: #else
 86:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex");
 87: #endif
 88:   MatDenseRestoreArray(jac->A,&a);
 89:   MatDenseRestoreArray(jac->U,&u);
 90:   MatDenseRestoreArray(jac->Vt,&v);
 91:   for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
 92:   jac->nzero = n-1-i;
 93:   if (jac->monitor) {
 94:     PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
 95:     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);
 96:     if (n >= 10) {              /* print 5 smallest and 5 largest */
 97:       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]));
 98:       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]));
 99:     } else {                    /* print all singular values */
100:       char     buf[256],*p;
101:       size_t   left = sizeof(buf),used;
102:       PetscInt thisline;
103:       for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
104:         PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));
105:         left -= used;
106:         p    += used;
107:         if (thisline > 4 || i==0) {
108:           PetscViewerASCIIPrintf(jac->monitor,"    SVD: singular values:%s\n",buf);
109:           p        = buf;
110:           thisline = 0;
111:         }
112:       }
113:     }
114:     PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
115:   }
116:   PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));
117:   for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
118:   for (; i<n; i++) d[i] = 0.0;
119:   if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
120:   PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
121:   VecRestoreArray(jac->diag,&d);
122: #if defined(foo)
123:   {
124:     PetscViewer viewer;
125:     PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);
126:     MatView(jac->A,viewer);
127:     MatView(jac->U,viewer);
128:     MatView(jac->Vt,viewer);
129:     VecView(jac->diag,viewer);
130:     PetscViewerDestroy(viewer);
131:   }
132: #endif
133:   PetscFree(work);
134:   return(0);
135: #endif
136: }

140: static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
141: {
142:   PC_SVD         *jac = (PC_SVD*)pc->data;
144:   PetscMPIInt    size;

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

179: static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
180: {
181:   PC_SVD         *jac = (PC_SVD*)pc->data;
183:   PetscMPIInt    size;

186:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
187:   switch (side) {
188:   case PC_LEFT:
189:     if (size != 1 && amode & WRITE) {
190:       VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
191:       VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);
192:     }
193:     break;
194:   case PC_RIGHT:
195:     if (size != 1 && amode & WRITE) {
196:       VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
197:       VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);
198:     }
199:     break;
200:   default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
201:   }
202:   *xred = NULL;
203:   return(0);
204: }

206: /* -------------------------------------------------------------------------- */
207: /*
208:    PCApply_SVD - Applies the SVD preconditioner to a vector.

210:    Input Parameters:
211: .  pc - the preconditioner context
212: .  x - input vector

214:    Output Parameter:
215: .  y - output vector

217:    Application Interface Routine: PCApply()
218:  */
221: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
222: {
223:   PC_SVD         *jac = (PC_SVD*)pc->data;
224:   Vec            work = jac->work,xred,yred;

228:   PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);
229:   PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);
230:   MatMultTranspose(jac->U,xred,work);
231:   VecPointwiseMult(work,work,jac->diag);
232:   MatMultTranspose(jac->Vt,work,yred);
233:   PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);
234:   PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);
235:   return(0);
236: }

240: static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
241: {
242:   PC_SVD         *jac = (PC_SVD*)pc->data;
243:   Vec            work = jac->work,xred,yred;

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

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

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

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

282:    Input Parameter:
283: .  pc - the preconditioner context

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

295:   PCReset_SVD(pc);
296:   PetscViewerDestroy(&jac->monitor);
297:   PetscFree(pc->data);
298:   return(0);
299: }

303: static PetscErrorCode PCSetFromOptions_SVD(PC pc)
304: {
306:   PC_SVD         *jac = (PC_SVD*)pc->data;
307:   PetscBool      flg,set;

310:   PetscOptionsHead("SVD options");
311:   PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);
312:   PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);
313:   PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
314:   if (set) {                    /* Should make PCSVDSetMonitor() */
315:     if (flg && !jac->monitor) {
316:       PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);
317:     } else if (!flg) {
318:       PetscViewerDestroy(&jac->monitor);
319:     }
320:   }
321:   PetscOptionsTail();
322:   return(0);
323: }

325: /* -------------------------------------------------------------------------- */
326: /*
327:    PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
328:    and sets this as the private data within the generic preconditioning
329:    context, PC, that was created within PCCreate().

331:    Input Parameter:
332: .  pc - the preconditioner context

334:    Application Interface Routine: PCCreate()
335: */

337: /*MC
338:      PCSVD - Use pseudo inverse defined by SVD of operator

340:    Level: advanced

342:   Concepts: SVD

344:   Options Database:
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: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC
349: M*/

353: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
354: {
355:   PC_SVD         *jac;

359:   /*
360:      Creates the private data structure for this preconditioner and
361:      attach it to the PC object.
362:   */
363:   PetscNewLog(pc,PC_SVD,&jac);
364:   jac->zerosing = 1.e-12;
365:   pc->data      = (void*)jac;

367:   /*
368:       Set the pointers for the functions that are provided above.
369:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
370:       are called, they will automatically call these functions.  Note we
371:       choose not to provide a couple of these functions since they are
372:       not needed.
373:   */
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            = 0;
381:   pc->ops->applyrichardson = 0;
382:   return(0);
383: }