Actual source code: fdmatrix.c

petsc-3.4.4 2014-03-13
  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;
 33:   PetscReal      x,y;

 36:   /* loop over colors  */
 37:   for (i=0; i<fd->ncolors; i++) {
 38:     for (j=0; j<fd->nrows[i]; j++) {
 39:       y    = fd->M - fd->rows[i][j] - fd->rstart;
 40:       x    = fd->columnsforrow[i][j];
 41:       PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);
 42:     }
 43:   }
 44:   return(0);
 45: }

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

 57:   PetscViewerDrawGetDraw(viewer,0,&draw);
 58:   PetscDrawIsNull(draw,&isnull); if (isnull) return(0);

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

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

 72: /*@C
 73:    MatFDColoringView - Views a finite difference coloring context.

 75:    Collective on MatFDColoring

 77:    Input  Parameters:
 78: +  c - the coloring context
 79: -  viewer - visualization context

 81:    Level: intermediate

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

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

 97: .seealso: MatFDColoringCreate()

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

110:   if (!viewer) {
111:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
112:   }

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

126:     PetscViewerGetFormat(viewer,&format);
127:     if (format != PETSC_VIEWER_ASCII_INFO) {
128:       for (i=0; i<c->ncolors; i++) {
129:         PetscViewerASCIIPrintf(viewer,"  Information for color %D\n",i);
130:         PetscViewerASCIIPrintf(viewer,"    Number of columns %D\n",c->ncolumns[i]);
131:         for (j=0; j<c->ncolumns[i]; j++) {
132:           PetscViewerASCIIPrintf(viewer,"      %D\n",c->columns[i][j]);
133:         }
134:         PetscViewerASCIIPrintf(viewer,"    Number of rows %D\n",c->nrows[i]);
135:         for (j=0; j<c->nrows[i]; j++) {
136:           PetscViewerASCIIPrintf(viewer,"      %D %D \n",c->rows[i][j],c->columnsforrow[i][j]);
137:         }
138:       }
139:     }
140:     PetscViewerFlush(viewer);
141:   }
142:   return(0);
143: }

147: /*@
148:    MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
149:    a Jacobian matrix using finite differences.

151:    Logically Collective on MatFDColoring

153:    The Jacobian is estimated with the differencing approximation
154: .vb
155:        F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
156:        h = error_rel*u[i]                 if  abs(u[i]) > umin
157:          = +/- error_rel*umin             otherwise, with +/- determined by the sign of u[i]
158:        dx_{i} = (0, ... 1, .... 0)
159: .ve

161:    Input Parameters:
162: +  coloring - the coloring context
163: .  error_rel - relative error
164: -  umin - minimum allowable u-value magnitude

166:    Level: advanced

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

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

172: @*/
173: PetscErrorCode  MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
174: {
179:   if (error != PETSC_DEFAULT) matfd->error_rel = error;
180:   if (umin != PETSC_DEFAULT)  matfd->umin      = umin;
181:   return(0);
182: }



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

191:    Not Collective

193:    Input Parameters:
194: .  coloring - the coloring context

196:    Output Parameters:
197: +  f - the function
198: -  fctx - the optional user-defined function context

200:    Level: intermediate

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

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

206: @*/
207: PetscErrorCode  MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
208: {
211:   if (f) *f = matfd->f;
212:   if (fctx) *fctx = matfd->fctx;
213:   return(0);
214: }

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

221:    Logically Collective on MatFDColoring

223:    Input Parameters:
224: +  coloring - the coloring context
225: .  f - the function
226: -  fctx - the optional user-defined function context

228:    Calling sequence of (*f) function:
229:     For SNES:    PetscErrorCode (*f)(SNES,Vec,Vec,void*)
230:     If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored

232:    Level: advanced

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

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

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

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

246: @*/
247: PetscErrorCode  MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
248: {
251:   matfd->f    = f;
252:   matfd->fctx = fctx;
253:   return(0);
254: }

258: /*@
259:    MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
260:    the options database.

262:    Collective on MatFDColoring

264:    The Jacobian, F'(u), is estimated with the differencing approximation
265: .vb
266:        F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
267:        h = error_rel*u[i]                 if  abs(u[i]) > umin
268:          = +/- error_rel*umin             otherwise, with +/- determined by the sign of u[i]
269:        dx_{i} = (0, ... 1, .... 0)
270: .ve

272:    Input Parameter:
273: .  coloring - the coloring context

275:    Options Database Keys:
276: +  -mat_fd_coloring_err <err> - Sets <err> (square root
277:            of relative error in the function)
278: .  -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
279: .  -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
280: .  -mat_fd_coloring_view - Activates basic viewing
281: .  -mat_fd_coloring_view ::ascii_info - Activates viewing info
282: -  -mat_fd_coloring_view draw - Activates drawing

284:     Level: intermediate

286: .keywords: Mat, finite differences, parameters

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

290: @*/
291: PetscErrorCode  MatFDColoringSetFromOptions(MatFDColoring matfd)
292: {
294:   PetscBool      flg;
295:   char           value[3];


300:   PetscObjectOptionsBegin((PetscObject)matfd);
301:   PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
302:   PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
303:   PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,3,&flg);
304:   if (flg) {
305:     if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
306:     else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
307:     else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
308:   }
309:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
310:   PetscObjectProcessOptionsHandlers((PetscObject)matfd);
311:   PetscOptionsEnd();
312:   return(0);
313: }

317: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
318: {
319:   PetscErrorCode    ierr;
320:   PetscBool         flg;
321:   PetscViewer       viewer;
322:   PetscViewerFormat format;

325:   if (prefix) {
326:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),prefix,optionname,&viewer,&format,&flg);
327:   } else {
328:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);
329:   }
330:   if (flg) {
331:     PetscViewerPushFormat(viewer,format);
332:     MatFDColoringView(fd,viewer);
333:     PetscViewerPopFormat(viewer);
334:     PetscViewerDestroy(&viewer);
335:   }
336:   return(0);
337: }

341: /*@
342:    MatFDColoringCreate - Creates a matrix coloring context for finite difference
343:    computation of Jacobians.

345:    Collective on Mat

347:    Input Parameters:
348: +  mat - the matrix containing the nonzero structure of the Jacobian
349: -  iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()

351:     Output Parameter:
352: .   color - the new coloring context

354:     Level: intermediate

356: .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
357:           MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
358:           MatFDColoringView(), MatFDColoringSetParameters(), MatGetColoring(), DMCreateColoring()
359: @*/
360: PetscErrorCode  MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
361: {
362:   MatFDColoring  c;
363:   MPI_Comm       comm;
365:   PetscInt       M,N;

368:   PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
369:   MatGetSize(mat,&M,&N);
370:   if (M != N) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Only for square matrices");
371:   PetscObjectGetComm((PetscObject)mat,&comm);
372:   PetscHeaderCreate(c,_p_MatFDColoring,int,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);

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

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

380:   MatGetVecs(mat,NULL,&c->w1);
381:   PetscLogObjectParent(c,c->w1);
382:   VecDuplicate(c->w1,&c->w2);
383:   PetscLogObjectParent(c,c->w2);

385:   c->error_rel    = PETSC_SQRT_MACHINE_EPSILON;
386:   c->umin         = 100.0*PETSC_SQRT_MACHINE_EPSILON;
387:   c->currentcolor = -1;
388:   c->htype        = "wp";
389:   c->fset         = PETSC_FALSE;

391:   *color = c;
392:   PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);
393:   PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
394:   return(0);
395: }

399: /*@
400:     MatFDColoringDestroy - Destroys a matrix coloring context that was created
401:     via MatFDColoringCreate().

403:     Collective on MatFDColoring

405:     Input Parameter:
406: .   c - coloring context

408:     Level: intermediate

410: .seealso: MatFDColoringCreate()
411: @*/
412: PetscErrorCode  MatFDColoringDestroy(MatFDColoring *c)
413: {
415:   PetscInt       i;

418:   if (!*c) return(0);
419:   if (--((PetscObject)(*c))->refct > 0) {*c = 0; return(0);}

421:   for (i=0; i<(*c)->ncolors; i++) {
422:     PetscFree((*c)->columns[i]);
423:     PetscFree((*c)->rows[i]);
424:     PetscFree((*c)->columnsforrow[i]);
425:     if ((*c)->vscaleforrow) {PetscFree((*c)->vscaleforrow[i]);}
426:   }
427:   PetscFree((*c)->ncolumns);
428:   PetscFree((*c)->columns);
429:   PetscFree((*c)->nrows);
430:   PetscFree((*c)->rows);
431:   PetscFree((*c)->columnsforrow);
432:   PetscFree((*c)->vscaleforrow);
433:   VecDestroy(&(*c)->vscale);
434:   VecDestroy(&(*c)->w1);
435:   VecDestroy(&(*c)->w2);
436:   VecDestroy(&(*c)->w3);
437:   PetscHeaderDestroy(c);
438:   return(0);
439: }

443: /*@C
444:     MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
445:       that are currently being perturbed.

447:     Not Collective

449:     Input Parameters:
450: .   coloring - coloring context created with MatFDColoringCreate()

452:     Output Parameters:
453: +   n - the number of local columns being perturbed
454: -   cols - the column indices, in global numbering

456:    Level: intermediate

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

460: .keywords: coloring, Jacobian, finite differences
461: @*/
462: PetscErrorCode  MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,PetscInt *cols[])
463: {
465:   if (coloring->currentcolor >= 0) {
466:     *n    = coloring->ncolumns[coloring->currentcolor];
467:     *cols = coloring->columns[coloring->currentcolor];
468:   } else {
469:     *n = 0;
470:   }
471:   return(0);
472: }

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

480:     Collective on MatFDColoring

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

488:     Options Database Keys:
489: +    -mat_fd_type - "wp" or "ds"  (see MATMFFD_WP or MATMFFD_DS)
490: .    -mat_fd_coloring_view - Activates basic viewing or coloring
491: .    -mat_fd_coloring_view draw - Activates drawing of coloring
492: -    -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info

494:     Level: intermediate

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

498: .keywords: coloring, Jacobian, finite differences
499: @*/
500: PetscErrorCode  MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
501: {

508:   if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
509:   if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
510:   (*J->ops->fdcoloringapply)(J,coloring,x1,flag,sctx);
511:   return(0);
512: }

516: PetscErrorCode  MatFDColoringApply_AIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
517: {
518:   PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f;
520:   PetscInt       k,start,end,l,row,col,srow,**vscaleforrow;
521:   PetscScalar    dx,*y,*xx,*w3_array;
522:   PetscScalar    *vscale_array;
523:   PetscReal      epsilon = coloring->error_rel,umin = coloring->umin,unorm;
524:   Vec            w1      = coloring->w1,w2=coloring->w2,w3;
525:   void           *fctx   = coloring->fctx;
526:   PetscBool      flg     = PETSC_FALSE;
527:   PetscInt       ctype   = coloring->ctype,N,col_start=0,col_end=0;
528:   Vec            x1_tmp;

534:   if (!f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");

536:   PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
537:   MatSetUnfactored(J);
538:   PetscOptionsGetBool(NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);
539:   if (flg) {
540:     PetscInfo(coloring,"Not calling MatZeroEntries()\n");
541:   } else {
542:     PetscBool assembled;
543:     MatAssembled(J,&assembled);
544:     if (assembled) {
545:       MatZeroEntries(J);
546:     }
547:   }

549:   x1_tmp = x1;
550:   if (!coloring->vscale) {
551:     VecDuplicate(x1_tmp,&coloring->vscale);
552:   }

554:   if (coloring->htype[0] == 'w') { /* tacky test; need to make systematic if we add other approaches to computing h*/
555:     VecNorm(x1_tmp,NORM_2,&unorm);
556:   }
557:   VecGetOwnershipRange(w1,&start,&end); /* OwnershipRange is used by ghosted x! */

559:   /* Set w1 = F(x1) */
560:   if (!coloring->fset) {
561:     PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
562:     (*f)(sctx,x1_tmp,w1,fctx);
563:     PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
564:   } else {
565:     coloring->fset = PETSC_FALSE;
566:   }

568:   if (!coloring->w3) {
569:     VecDuplicate(x1_tmp,&coloring->w3);
570:     PetscLogObjectParent(coloring,coloring->w3);
571:   }
572:   w3 = coloring->w3;

574:   /* Compute all the local scale factors, including ghost points */
575:   VecGetLocalSize(x1_tmp,&N);
576:   VecGetArray(x1_tmp,&xx);
577:   VecGetArray(coloring->vscale,&vscale_array);
578:   if (ctype == IS_COLORING_GHOSTED) {
579:     col_start = 0; col_end = N;
580:   } else if (ctype == IS_COLORING_GLOBAL) {
581:     xx           = xx - start;
582:     vscale_array = vscale_array - start;
583:     col_start    = start; col_end = N + start;
584:   }
585:   for (col=col_start; col<col_end; col++) {
586:     /* Loop over each local column, vscale[col] = 1./(epsilon*dx[col]) */
587:     if (coloring->htype[0] == 'w') {
588:       dx = 1.0 + unorm;
589:     } else {
590:       dx = xx[col];
591:     }
592:     if (dx == (PetscScalar)0.0) dx = 1.0;
593:     if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
594:     else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
595:     dx               *= epsilon;
596:     vscale_array[col] = (PetscScalar)1.0/dx;
597:   }
598:   if (ctype == IS_COLORING_GLOBAL) vscale_array = vscale_array + start;
599:   VecRestoreArray(coloring->vscale,&vscale_array);
600:   if (ctype == IS_COLORING_GLOBAL) {
601:     VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
602:     VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
603:   }

605:   if (coloring->vscaleforrow) {
606:     vscaleforrow = coloring->vscaleforrow;
607:   } else SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_NULL,"Null Object: coloring->vscaleforrow");

609:   /*
610:     Loop over each color
611:   */
612:   VecGetArray(coloring->vscale,&vscale_array);
613:   for (k=0; k<coloring->ncolors; k++) {
614:     coloring->currentcolor = k;

616:     VecCopy(x1_tmp,w3);
617:     VecGetArray(w3,&w3_array);
618:     if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array - start;
619:     /*
620:       Loop over each column associated with color
621:       adding the perturbation to the vector w3.
622:     */
623:     for (l=0; l<coloring->ncolumns[k]; l++) {
624:       col = coloring->columns[k][l];    /* local column of the matrix we are probing for */
625:       if (coloring->htype[0] == 'w') {
626:         dx = 1.0 + unorm;
627:       } else {
628:         dx = xx[col];
629:       }
630:       if (dx == (PetscScalar)0.0) dx = 1.0;
631:       if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
632:       else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
633:       dx            *= epsilon;
634:       if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter");
635:       w3_array[col] += dx;
636:     }
637:     if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array + start;
638:     VecRestoreArray(w3,&w3_array);

640:     /*
641:       Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
642:                            w2 = F(x1 + dx) - F(x1)
643:     */
644:     PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
645:     (*f)(sctx,w3,w2,fctx);
646:     PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
647:     VecAXPY(w2,-1.0,w1);

649:     /*
650:       Loop over rows of vector, putting results into Jacobian matrix
651:     */
652:     VecGetArray(w2,&y);
653:     for (l=0; l<coloring->nrows[k]; l++) {
654:       row     = coloring->rows[k][l];            /* local row index */
655:       col     = coloring->columnsforrow[k][l];   /* global column index */
656:       y[row] *= vscale_array[vscaleforrow[k][l]];
657:       srow    = row + start;
658:       MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);
659:     }
660:     VecRestoreArray(w2,&y);
661:   } /* endof for each color */
662:   if (ctype == IS_COLORING_GLOBAL) xx = xx + start;
663:   VecRestoreArray(coloring->vscale,&vscale_array);
664:   VecRestoreArray(x1_tmp,&xx);

666:   coloring->currentcolor = -1;

668:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
669:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
670:   PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);

672:   MatFDColoringViewFromOptions(coloring,NULL,"-mat_fd_coloring_view");
673:   return(0);
674: }