Actual source code: fdmatrix.c

petsc-3.7.2 2016-06-05
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*/
  8: #include <petsc/private/isimpl.h>

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

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

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

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

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

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

 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:   PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);
 69:   PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
 70:   PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);
 71:   PetscDrawSave(draw);
 72:   return(0);
 73: }

 77: /*@C
 78:    MatFDColoringView - Views a finite difference coloring context.

 80:    Collective on MatFDColoring

 82:    Input  Parameters:
 83: +  c - the coloring context
 84: -  viewer - visualization context

 86:    Level: intermediate

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

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

102: .seealso: MatFDColoringCreate()

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

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

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

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

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

160:    Logically Collective on MatFDColoring

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

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

175:    Level: advanced

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

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

181: @*/
182: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
183: {
188:   if (error != PETSC_DEFAULT) matfd->error_rel = error;
189:   if (umin != PETSC_DEFAULT)  matfd->umin      = umin;
190:   return(0);
191: }

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

198:    Logically Collective on MatFDColoring

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

205:    Level: intermediate

207: .keywords: Mat, coloring

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

211: @*/
212: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
213: {
218:   if (brows != PETSC_DEFAULT) matfd->brows = brows;
219:   if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
220:   return(0);
221: }

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

228:    Collective on Mat

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

235:    Level: beginner

237: .keywords: MatFDColoring, setup

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

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

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

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

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

265:    Not Collective

267:    Input Parameters:
268: .  coloring - the coloring context

270:    Output Parameters:
271: +  f - the function
272: -  fctx - the optional user-defined function context

274:    Level: intermediate

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

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

280: @*/
281: PetscErrorCode  MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
282: {
285:   if (f) *f = matfd->f;
286:   if (fctx) *fctx = matfd->fctx;
287:   return(0);
288: }

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

295:    Logically Collective on MatFDColoring

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

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

306:    Level: advanced

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

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

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

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

320: @*/
321: PetscErrorCode  MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
322: {
325:   matfd->f    = f;
326:   matfd->fctx = fctx;
327:   return(0);
328: }

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

336:    Collective on MatFDColoring

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

346:    Input Parameter:
347: .  coloring - the coloring context

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

357:     Level: intermediate

359: .keywords: Mat, finite differences, parameters

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

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


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

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

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

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

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

425:    Collective on Mat

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

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

434:     Level: intermediate

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

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

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

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

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

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

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

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

486:     Collective on MatFDColoring

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

491:     Level: intermediate

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

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

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

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

531:     Not Collective

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

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

540:    Level: intermediate

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

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

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

564:     Collective on MatFDColoring

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

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

578:     Level: intermediate

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

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

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

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

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