Actual source code: fdmatrix.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  2: /*
  3:    This is where the abstract matrix operations are defined that are
  4:   used for finite difference computations of Jacobians using coloring.
  5: */

  7: #include <petsc-private/matimpl.h>        /*I "petscmat.h" I*/

 11: PetscErrorCode  MatFDColoringSetF(MatFDColoring fd,Vec F)
 12: {

 16:   if (F) {
 17:     VecCopy(F,fd->w1);
 18:     fd->fset = PETSC_TRUE;
 19:   } else {
 20:     fd->fset = PETSC_FALSE;
 21:   }
 22:   return(0);
 23: }

 25: #include <petscdraw.h>
 28: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
 29: {
 30:   MatFDColoring  fd = (MatFDColoring)Aa;
 32:   PetscInt       i,j,nz,row;
 33:   PetscReal      x,y;
 34:   MatEntry       *Jentry=fd->matentry;

 37:   /* loop over colors  */
 38:   nz = 0;
 39:   for (i=0; i<fd->ncolors; i++) {
 40:     for (j=0; j<fd->nrows[i]; j++) {
 41:       row = Jentry[nz].row;
 42:       y   = fd->M - row - fd->rstart;
 43:       x   = (PetscReal)Jentry[nz++].col;
 44:       PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);
 45:     }
 46:   }
 47:   return(0);
 48: }

 52: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
 53: {
 55:   PetscBool      isnull;
 56:   PetscDraw      draw;
 57:   PetscReal      xr,yr,xl,yl,h,w;

 60:   PetscViewerDrawGetDraw(viewer,0,&draw);
 61:   PetscDrawIsNull(draw,&isnull); if (isnull) return(0);

 63:   PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);

 65:   xr   = fd->N; yr  = fd->M; h = yr/10.0; w = xr/10.0;
 66:   xr  += w;     yr += h;    xl = -w;     yl = -h;
 67:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
 68:   PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
 69:   PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);
 70:   return(0);
 71: }

 75: /*@C
 76:    MatFDColoringView - Views a finite difference coloring context.

 78:    Collective on MatFDColoring

 80:    Input  Parameters:
 81: +  c - the coloring context
 82: -  viewer - visualization context

 84:    Level: intermediate

 86:    Notes:
 87:    The available visualization contexts include
 88: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 89: .     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 90:         output where only the first processor opens
 91:         the file.  All other processors send their
 92:         data to the first processor to print.
 93: -     PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure

 95:    Notes:
 96:      Since PETSc uses only a small number of basic colors (currently 33), if the coloring
 97:    involves more than 33 then some seemingly identical colors are displayed making it look
 98:    like an illegal coloring. This is just a graphical artifact.

100: .seealso: MatFDColoringCreate()

102: .keywords: Mat, finite differences, coloring, view
103: @*/
104: PetscErrorCode  MatFDColoringView(MatFDColoring c,PetscViewer viewer)
105: {
106:   PetscErrorCode    ierr;
107:   PetscInt          i,j;
108:   PetscBool         isdraw,iascii;
109:   PetscViewerFormat format;

113:   if (!viewer) {
114:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
115:   }

119:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
120:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
121:   if (isdraw) {
122:     MatFDColoringView_Draw(c,viewer);
123:   } else if (iascii) {
124:     PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer);
125:     PetscViewerASCIIPrintf(viewer,"  Error tolerance=%g\n",(double)c->error_rel);
126:     PetscViewerASCIIPrintf(viewer,"  Umin=%g\n",(double)c->umin);
127:     PetscViewerASCIIPrintf(viewer,"  Number of colors=%D\n",c->ncolors);

129:     PetscViewerGetFormat(viewer,&format);
130:     if (format != PETSC_VIEWER_ASCII_INFO) {
131:       PetscInt row,col,nz;
132:       nz = 0;
133:       for (i=0; i<c->ncolors; i++) {
134:         PetscViewerASCIIPrintf(viewer,"  Information for color %D\n",i);
135:         PetscViewerASCIIPrintf(viewer,"    Number of columns %D\n",c->ncolumns[i]);
136:         for (j=0; j<c->ncolumns[i]; j++) {
137:           PetscViewerASCIIPrintf(viewer,"      %D\n",c->columns[i][j]);
138:         }
139:         PetscViewerASCIIPrintf(viewer,"    Number of rows %D\n",c->nrows[i]);
140:         for (j=0; j<c->nrows[i]; j++) {
141:           row  = c->matentry[nz].row;
142:           col  = c->matentry[nz++].col;
143:           PetscViewerASCIIPrintf(viewer,"      %D %D \n",row,col);
144:         }
145:       }
146:     }
147:     PetscViewerFlush(viewer);
148:   }
149:   return(0);
150: }

154: /*@
155:    MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
156:    a Jacobian matrix using finite differences.

158:    Logically Collective on MatFDColoring

160:    The Jacobian is estimated with the differencing approximation
161: .vb
162:        F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
163:        h = error_rel*u[i]                 if  abs(u[i]) > umin
164:          = +/- error_rel*umin             otherwise, with +/- determined by the sign of u[i]
165:        dx_{i} = (0, ... 1, .... 0)
166: .ve

168:    Input Parameters:
169: +  coloring - the coloring context
170: .  error_rel - relative error
171: -  umin - minimum allowable u-value magnitude

173:    Level: advanced

175: .keywords: Mat, finite differences, coloring, set, parameters

177: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()

179: @*/
180: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
181: {
186:   if (error != PETSC_DEFAULT) matfd->error_rel = error;
187:   if (umin != PETSC_DEFAULT)  matfd->umin      = umin;
188:   return(0);
189: }

193: /*@
194:    MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.

196:    Logically Collective on MatFDColoring

198:    Input Parameters:
199: +  coloring - the coloring context
200: .  brows - number of rows in the block
201: -  bcols - number of columns in the block

203:    Level: intermediate

205: .keywords: Mat, coloring

207: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()

209: @*/
210: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
211: {
216:   if (brows != PETSC_DEFAULT) matfd->brows = brows;
217:   if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
218:   return(0);
219: }

223: /*@
224:    MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.

226:    Collective on Mat

228:    Input Parameters:
229: +  mat - the matrix containing the nonzero structure of the Jacobian
230: .  iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
231: -  color - the matrix coloring context

233:    Level: beginner

235: .keywords: MatFDColoring, setup

237: .seealso: MatFDColoringCreate(), MatFDColoringDestroy()
238: @*/
239: PetscErrorCode MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)
240: {

246:   if (color->setupcalled) return(0);

248:   PetscLogEventBegin(MAT_FDColoringSetUp,mat,0,0,0);
249:   if (mat->ops->fdcoloringsetup) {
250:     (*mat->ops->fdcoloringsetup)(mat,iscoloring,color);
251:   } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);

253:   color->setupcalled = PETSC_TRUE;
254:    PetscLogEventEnd(MAT_FDColoringSetUp,mat,0,0,0);
255:   return(0);
256: }

260: /*@C
261:    MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.

263:    Not Collective

265:    Input Parameters:
266: .  coloring - the coloring context

268:    Output Parameters:
269: +  f - the function
270: -  fctx - the optional user-defined function context

272:    Level: intermediate

274: .keywords: Mat, Jacobian, finite differences, set, function

276: .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()

278: @*/
279: PetscErrorCode  MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
280: {
283:   if (f) *f = matfd->f;
284:   if (fctx) *fctx = matfd->fctx;
285:   return(0);
286: }

290: /*@C
291:    MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.

293:    Logically Collective on MatFDColoring

295:    Input Parameters:
296: +  coloring - the coloring context
297: .  f - the function
298: -  fctx - the optional user-defined function context

300:    Calling sequence of (*f) function:
301:     For SNES:    PetscErrorCode (*f)(SNES,Vec,Vec,void*)
302:     If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored

304:    Level: advanced

306:    Notes: This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
307:      SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
308:      calling MatFDColoringApply()

310:    Fortran Notes:
311:     In Fortran you must call MatFDColoringSetFunction() for a coloring object to
312:   be used without SNES or within the SNES solvers.

314: .keywords: Mat, Jacobian, finite differences, set, function

316: .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()

318: @*/
319: PetscErrorCode  MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
320: {
323:   matfd->f    = f;
324:   matfd->fctx = fctx;
325:   return(0);
326: }

330: /*@
331:    MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
332:    the options database.

334:    Collective on MatFDColoring

336:    The Jacobian, F'(u), is estimated with the differencing approximation
337: .vb
338:        F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
339:        h = error_rel*u[i]                 if  abs(u[i]) > umin
340:          = +/- error_rel*umin             otherwise, with +/- determined by the sign of u[i]
341:        dx_{i} = (0, ... 1, .... 0)
342: .ve

344:    Input Parameter:
345: .  coloring - the coloring context

347:    Options Database Keys:
348: +  -mat_fd_coloring_err <err> - Sets <err> (square root
349:            of relative error in the function)
350: .  -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
351: .  -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
352: .  -mat_fd_coloring_view - Activates basic viewing
353: .  -mat_fd_coloring_view ::ascii_info - Activates viewing info
354: -  -mat_fd_coloring_view draw - Activates drawing

356:     Level: intermediate

358: .keywords: Mat, finite differences, parameters

360: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()

362: @*/
363: PetscErrorCode  MatFDColoringSetFromOptions(MatFDColoring matfd)
364: {
366:   PetscBool      flg;
367:   char           value[3];


372:   PetscObjectOptionsBegin((PetscObject)matfd);
373:   PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
374:   PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
375:   PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,3,&flg);
376:   if (flg) {
377:     if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
378:     else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
379:     else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
380:   }
381:   PetscOptionsInt("-mat_fd_coloring_brows","Number of block rows","MatFDColoringSetBlockSize",matfd->brows,&matfd->brows,NULL);
382:   PetscOptionsInt("-mat_fd_coloring_bcols","Number of block columns","MatFDColoringSetBlockSize",matfd->bcols,&matfd->bcols,&flg);
383:   if (flg && matfd->bcols > matfd->ncolors) {
384:     /* input bcols cannot be > matfd->ncolors, thus set it as ncolors */
385:     matfd->bcols = matfd->ncolors;
386:   }

388:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
389:   PetscObjectProcessOptionsHandlers((PetscObject)matfd);
390:   PetscOptionsEnd();
391:   return(0);
392: }

396: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
397: {
398:   PetscErrorCode    ierr;
399:   PetscBool         flg;
400:   PetscViewer       viewer;
401:   PetscViewerFormat format;

404:   if (prefix) {
405:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),prefix,optionname,&viewer,&format,&flg);
406:   } else {
407:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);
408:   }
409:   if (flg) {
410:     PetscViewerPushFormat(viewer,format);
411:     MatFDColoringView(fd,viewer);
412:     PetscViewerPopFormat(viewer);
413:     PetscViewerDestroy(&viewer);
414:   }
415:   return(0);
416: }

420: /*@
421:    MatFDColoringCreate - Creates a matrix coloring context for finite difference
422:    computation of Jacobians.

424:    Collective on Mat

426:    Input Parameters:
427: +  mat - the matrix containing the nonzero structure of the Jacobian
428: -  iscoloring - the coloring of the matrix; usually obtained with MatColoringCreate() or DMCreateColoring()

430:     Output Parameter:
431: .   color - the new coloring context

433:     Level: intermediate

435: .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
436:           MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
437:           MatFDColoringView(), MatFDColoringSetParameters(), MatColoringCreate(), DMCreateColoring()
438: @*/
439: PetscErrorCode  MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
440: {
441:   MatFDColoring  c;
442:   MPI_Comm       comm;
444:   PetscInt       M,N;

448:   if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();");
449:   PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
450:   MatGetSize(mat,&M,&N);
451:   if (M != N) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Only for square matrices");
452:   PetscObjectGetComm((PetscObject)mat,&comm);
453:   PetscHeaderCreate(c,_p_MatFDColoring,int,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);

455:   c->ctype = iscoloring->ctype;

457:   if (mat->ops->fdcoloringcreate) {
458:     (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
459:   } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);

461:   MatGetVecs(mat,NULL,&c->w1);
462:   PetscLogObjectParent((PetscObject)c,(PetscObject)c->w1);
463:   VecDuplicate(c->w1,&c->w2);
464:   PetscLogObjectParent((PetscObject)c,(PetscObject)c->w2);

466:   c->error_rel    = PETSC_SQRT_MACHINE_EPSILON;
467:   c->umin         = 100.0*PETSC_SQRT_MACHINE_EPSILON;
468:   c->currentcolor = -1;
469:   c->htype        = "wp";
470:   c->fset         = PETSC_FALSE;
471:   c->setupcalled  = PETSC_FALSE;

473:   *color = c;
474:   PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);
475:   PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
476:   return(0);
477: }

481: /*@
482:     MatFDColoringDestroy - Destroys a matrix coloring context that was created
483:     via MatFDColoringCreate().

485:     Collective on MatFDColoring

487:     Input Parameter:
488: .   c - coloring context

490:     Level: intermediate

492: .seealso: MatFDColoringCreate()
493: @*/
494: PetscErrorCode  MatFDColoringDestroy(MatFDColoring *c)
495: {
497:   PetscInt       i;
498:   MatFDColoring  color = *c;

501:   if (!*c) return(0);
502:   if (--((PetscObject)color)->refct > 0) {*c = 0; return(0);}

504:   for (i=0; i<color->ncolors; i++) {
505:     PetscFree(color->columns[i]);
506:   }
507:   PetscFree(color->ncolumns);
508:   PetscFree(color->columns);
509:   PetscFree(color->nrows);
510:   if (color->htype[0] == 'w') {
511:     PetscFree(color->matentry2);
512:   } else {
513:     PetscFree(color->matentry);
514:   }
515:   PetscFree(color->dy);
516:   if (color->vscale) {VecDestroy(&color->vscale);}
517:   VecDestroy(&color->w1);
518:   VecDestroy(&color->w2);
519:   VecDestroy(&color->w3);
520:   PetscHeaderDestroy(c);
521:   return(0);
522: }

526: /*@C
527:     MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
528:       that are currently being perturbed.

530:     Not Collective

532:     Input Parameters:
533: .   coloring - coloring context created with MatFDColoringCreate()

535:     Output Parameters:
536: +   n - the number of local columns being perturbed
537: -   cols - the column indices, in global numbering

539:    Level: intermediate

541: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()

543: .keywords: coloring, Jacobian, finite differences
544: @*/
545: PetscErrorCode  MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,PetscInt *cols[])
546: {
548:   if (coloring->currentcolor >= 0) {
549:     *n    = coloring->ncolumns[coloring->currentcolor];
550:     *cols = coloring->columns[coloring->currentcolor];
551:   } else {
552:     *n = 0;
553:   }
554:   return(0);
555: }

559: /*@
560:     MatFDColoringApply - Given a matrix for which a MatFDColoring context
561:     has been created, computes the Jacobian for a function via finite differences.

563:     Collective on MatFDColoring

565:     Input Parameters:
566: +   mat - location to store Jacobian
567: .   coloring - coloring context created with MatFDColoringCreate()
568: .   x1 - location at which Jacobian is to be computed
569: -   sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null

571:     Options Database Keys:
572: +    -mat_fd_type - "wp" or "ds"  (see MATMFFD_WP or MATMFFD_DS)
573: .    -mat_fd_coloring_view - Activates basic viewing or coloring
574: .    -mat_fd_coloring_view draw - Activates drawing of coloring
575: -    -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info

577:     Level: intermediate

579: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction()

581: .keywords: coloring, Jacobian, finite differences
582: @*/
583: PetscErrorCode  MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
584: {
586:   PetscBool      flg = PETSC_FALSE;

592:   if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
593:   if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
594:   if (!coloring->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetUp()");

596:   MatSetUnfactored(J);
597:   PetscOptionsGetBool(NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);
598:   if (flg) {
599:     PetscInfo(coloring,"Not calling MatZeroEntries()\n");
600:   } else {
601:     PetscBool assembled;
602:     MatAssembled(J,&assembled);
603:     if (assembled) {
604:       MatZeroEntries(J);
605:     }
606:   }

608:   PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
609:   (*J->ops->fdcoloringapply)(J,coloring,x1,sctx);
610:   PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);
611:   return(0);
612: }