Actual source code: fdmatrix.c

  1: #define PETSCMAT_DLL

  3: /*
  4:    This is where the abstract matrix operations are defined that are
  5:   used for finite difference computations of Jacobians using coloring.
  6: */

 8:  #include include/private/matimpl.h

 10: /* Logging support */
 11: PetscCookie  MAT_FDCOLORING_COOKIE = 0;

 15: PetscErrorCode  MatFDColoringSetF(MatFDColoring fd,Vec F)
 16: {
 18:   fd->F = F;
 19:   return(0);
 20: }

 24: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
 25: {
 26:   MatFDColoring  fd = (MatFDColoring)Aa;
 28:   PetscInt       i,j;
 29:   PetscReal      x,y;


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

 46: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
 47: {
 49:   PetscTruth     isnull;
 50:   PetscDraw      draw;
 51:   PetscReal      xr,yr,xl,yl,h,w;

 54:   PetscViewerDrawGetDraw(viewer,0,&draw);
 55:   PetscDrawIsNull(draw,&isnull); if (isnull) return(0);

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

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

 69: /*@C
 70:    MatFDColoringView - Views a finite difference coloring context.

 72:    Collective on MatFDColoring

 74:    Input  Parameters:
 75: +  c - the coloring context
 76: -  viewer - visualization context

 78:    Level: intermediate

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

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

 94: .seealso: MatFDColoringCreate()

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

107:   if (!viewer) {
108:     PetscViewerASCIIGetStdout(c->comm,&viewer);
109:   }

113:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
114:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
115:   if (isdraw) {
116:     MatFDColoringView_Draw(c,viewer);
117:   } else if (iascii) {
118:     PetscViewerASCIIPrintf(viewer,"MatFDColoring Object:\n");
119:     PetscViewerASCIIPrintf(viewer,"  Error tolerance=%G\n",c->error_rel);
120:     PetscViewerASCIIPrintf(viewer,"  Umin=%G\n",c->umin);
121:     PetscViewerASCIIPrintf(viewer,"  Number of colors=%D\n",c->ncolors);

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

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

150:    Collective on MatFDColoring

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

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

165:    Level: advanced

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

169: .seealso: MatFDColoringCreate()
170: @*/
171: PetscErrorCode  MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
172: {

176:   if (error != PETSC_DEFAULT) matfd->error_rel = error;
177:   if (umin != PETSC_DEFAULT)  matfd->umin      = umin;
178:   return(0);
179: }

183: /*@
184:    MatFDColoringSetFrequency - Sets the frequency for computing new Jacobian
185:    matrices. 

187:    Collective on MatFDColoring

189:    Input Parameters:
190: +  coloring - the coloring context
191: -  freq - frequency (default is 1)

193:    Options Database Keys:
194: .  -mat_fd_coloring_freq <freq>  - Sets coloring frequency

196:    Level: advanced

198:    Notes:
199:    Using a modified Newton strategy, where the Jacobian remains fixed for several
200:    iterations, can be cost effective in terms of overall nonlinear solution 
201:    efficiency.  This parameter indicates that a new Jacobian will be computed every
202:    <freq> nonlinear iterations.  

204: .keywords: Mat, finite differences, coloring, set, frequency

206: .seealso: MatFDColoringCreate(), MatFDColoringGetFrequency(), MatFDColoringSetRecompute()
207: @*/
208: PetscErrorCode  MatFDColoringSetFrequency(MatFDColoring matfd,PetscInt freq)
209: {

213:   matfd->freq = freq;
214:   return(0);
215: }

219: /*@
220:    MatFDColoringGetFrequency - Gets the frequency for computing new Jacobian
221:    matrices. 

223:    Not Collective

225:    Input Parameters:
226: .  coloring - the coloring context

228:    Output Parameters:
229: .  freq - frequency (default is 1)

231:    Options Database Keys:
232: .  -mat_fd_coloring_freq <freq> - Sets coloring frequency

234:    Level: advanced

236:    Notes:
237:    Using a modified Newton strategy, where the Jacobian remains fixed for several
238:    iterations, can be cost effective in terms of overall nonlinear solution 
239:    efficiency.  This parameter indicates that a new Jacobian will be computed every
240:    <freq> nonlinear iterations.  

242: .keywords: Mat, finite differences, coloring, get, frequency

244: .seealso: MatFDColoringSetFrequency()
245: @*/
246: PetscErrorCode  MatFDColoringGetFrequency(MatFDColoring matfd,PetscInt *freq)
247: {
250:   *freq = matfd->freq;
251:   return(0);
252: }

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

259:    Collective on MatFDColoring

261:    Input Parameters:
262: .  coloring - the coloring context

264:    Output Parameters:
265: +  f - the function
266: -  fctx - the optional user-defined function context

268:    Level: intermediate

270: .keywords: Mat, Jacobian, finite differences, set, function
271: @*/
272: PetscErrorCode  MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
273: {
276:   if (f) *f = matfd->f;
277:   if (fctx) *fctx = matfd->fctx;
278:   return(0);
279: }

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

286:    Collective on MatFDColoring

288:    Input Parameters:
289: +  coloring - the coloring context
290: .  f - the function
291: -  fctx - the optional user-defined function context

293:    Level: intermediate

295:    Notes:
296:     In Fortran you must call MatFDColoringSetFunctionSNES() for a coloring object to 
297:   be used with the SNES solvers and MatFDColoringSetFunctionTS() if it is to be used
298:   with the TS solvers.

300: .keywords: Mat, Jacobian, finite differences, set, function
301: @*/
302: PetscErrorCode  MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
303: {
306:   matfd->f    = f;
307:   matfd->fctx = fctx;
308:   return(0);
309: }

313: /*@
314:    MatFDColoringSetFromOptions - Sets coloring finite difference parameters from 
315:    the options database.

317:    Collective on MatFDColoring

319:    The Jacobian, F'(u), is estimated with the differencing approximation
320: .vb
321:        F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
322:        h = error_rel*u[i]                 if  abs(u[i]) > umin
323:          = +/- error_rel*umin             otherwise, with +/- determined by the sign of u[i]
324:        dx_{i} = (0, ... 1, .... 0)
325: .ve

327:    Input Parameter:
328: .  coloring - the coloring context

330:    Options Database Keys:
331: +  -mat_fd_coloring_err <err> - Sets <err> (square root
332:            of relative error in the function)
333: .  -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
334: .  -mat_fd_coloring_freq <freq> - Sets frequency of computing a new Jacobian
335: .  -mat_fd_type - "wp" or "ds" (see MATSNESMF_WP or MATSNESMF_DS)
336: .  -mat_fd_coloring_view - Activates basic viewing
337: .  -mat_fd_coloring_view_info - Activates viewing info
338: -  -mat_fd_coloring_view_draw - Activates drawing

340:     Level: intermediate

342: .keywords: Mat, finite differences, parameters

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

346: @*/
347: PetscErrorCode  MatFDColoringSetFromOptions(MatFDColoring matfd)
348: {
350:   PetscTruth     flg;
351:   char           value[3];


356:   PetscOptionsBegin(matfd->comm,matfd->prefix,"Jacobian computation via finite differences option","MatFD");
357:     PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
358:     PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
359:     PetscOptionsInt("-mat_fd_coloring_freq","How often Jacobian is recomputed","MatFDColoringSetFrequency",matfd->freq,&matfd->freq,0);
360:     PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,2,&flg);
361:     if (flg) {
362:       if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
363:       else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
364:       else SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
365:     }
366:     /* not used here; but so they are presented in the GUI */
367:     PetscOptionsName("-mat_fd_coloring_view","Print entire datastructure used for Jacobian","None",0);
368:     PetscOptionsName("-mat_fd_coloring_view_info","Print number of colors etc for Jacobian","None",0);
369:     PetscOptionsName("-mat_fd_coloring_view_draw","Plot nonzero structure ofJacobian","None",0);
370:   PetscOptionsEnd();
371:   return(0);
372: }

376: PetscErrorCode MatFDColoringView_Private(MatFDColoring fd)
377: {
379:   PetscTruth     flg;
380:   PetscViewer    viewer;

383:   PetscViewerASCIIGetStdout(fd->comm,&viewer);
384:   PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_view",&flg);
385:   if (flg) {
386:     MatFDColoringView(fd,viewer);
387:   }
388:   PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_info",&flg);
389:   if (flg) {
390:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
391:     MatFDColoringView(fd,viewer);
392:     PetscViewerPopFormat(viewer);
393:   }
394:   PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_draw",&flg);
395:   if (flg) {
396:     MatFDColoringView(fd,PETSC_VIEWER_DRAW_(fd->comm));
397:     PetscViewerFlush(PETSC_VIEWER_DRAW_(fd->comm));
398:   }
399:   return(0);
400: }

404: /*@
405:    MatFDColoringCreate - Creates a matrix coloring context for finite difference 
406:    computation of Jacobians.

408:    Collective on Mat

410:    Input Parameters:
411: +  mat - the matrix containing the nonzero structure of the Jacobian
412: -  iscoloring - the coloring of the matrix

414:     Output Parameter:
415: .   color - the new coloring context
416:    
417:     Level: intermediate

419: .seealso: MatFDColoringDestroy(),SNESDefaultComputeJacobianColor(), ISColoringCreate(),
420:           MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
421:           MatFDColoringSetFrequency(), MatFDColoringSetRecompute(), MatFDColoringView(),
422:           MatFDColoringSetParameters()
423: @*/
424: PetscErrorCode  MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
425: {
426:   MatFDColoring  c;
427:   MPI_Comm       comm;
429:   PetscInt       M,N;
430:   PetscMPIInt    size;

434:   MatGetSize(mat,&M,&N);
435:   if (M != N) SETERRQ(PETSC_ERR_SUP,"Only for square matrices");

437:   PetscObjectGetComm((PetscObject)mat,&comm);
438:   PetscHeaderCreate(c,_p_MatFDColoring,int,MAT_FDCOLORING_COOKIE,0,"MatFDColoring",comm,MatFDColoringDestroy,MatFDColoringView);

440:   MPI_Comm_size(comm,&size);
441:   c->ctype = iscoloring->ctype;

443:   if (mat->ops->fdcoloringcreate) {
444:     (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
445:   } else {
446:     SETERRQ(PETSC_ERR_SUP,"Code not yet written for this matrix type");
447:   }

449:   MatGetVecs(mat,PETSC_NULL,&c->w1);
450:   PetscLogObjectParent(c,c->w1);
451:   VecDuplicate(c->w1,&c->w2);
452:   PetscLogObjectParent(c,c->w2);

454:   c->error_rel         = PETSC_SQRT_MACHINE_EPSILON;
455:   c->umin              = 100.0*PETSC_SQRT_MACHINE_EPSILON;
456:   c->freq              = 1;
457:   c->usersetsrecompute = PETSC_FALSE;
458:   c->recompute         = PETSC_FALSE;
459:   c->currentcolor      = -1;
460:   c->htype             = "wp";

462:   *color = c;
464:   return(0);
465: }

469: /*@
470:     MatFDColoringDestroy - Destroys a matrix coloring context that was created
471:     via MatFDColoringCreate().

473:     Collective on MatFDColoring

475:     Input Parameter:
476: .   c - coloring context

478:     Level: intermediate

480: .seealso: MatFDColoringCreate()
481: @*/
482: PetscErrorCode  MatFDColoringDestroy(MatFDColoring c)
483: {
485:   PetscInt       i;

488:   if (--c->refct > 0) return(0);

490:   for (i=0; i<c->ncolors; i++) {
491:     PetscFree(c->columns[i]);
492:     PetscFree(c->rows[i]);
493:     PetscFree(c->columnsforrow[i]);
494:     if (c->vscaleforrow) {PetscFree(c->vscaleforrow[i]);}
495:   }
496:   PetscFree(c->ncolumns);
497:   PetscFree(c->columns);
498:   PetscFree(c->nrows);
499:   PetscFree(c->rows);
500:   PetscFree(c->columnsforrow);
501:   PetscFree(c->vscaleforrow);
502:   if (c->vscale)       {VecDestroy(c->vscale);}
503:   if (c->w1) {
504:     VecDestroy(c->w1);
505:     VecDestroy(c->w2);
506:   }
507:   if (c->w3){
508:     VecDestroy(c->w3);
509:   }
510:   PetscHeaderDestroy(c);
511:   return(0);
512: }

516: /*@C
517:     MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
518:       that are currently being perturbed.

520:     Not Collective

522:     Input Parameters:
523: .   coloring - coloring context created with MatFDColoringCreate()

525:     Output Parameters:
526: +   n - the number of local columns being perturbed
527: -   cols - the column indices, in global numbering

529:    Level: intermediate

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

533: .keywords: coloring, Jacobian, finite differences
534: @*/
535: PetscErrorCode  MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,PetscInt *cols[])
536: {
538:   if (coloring->currentcolor >= 0) {
539:     *n    = coloring->ncolumns[coloring->currentcolor];
540:     *cols = coloring->columns[coloring->currentcolor];
541:   } else {
542:     *n = 0;
543:   }
544:   return(0);
545: }

547: #include "petscda.h"      /*I      "petscda.h"    I*/ 

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

555:     Collective on MatFDColoring

557:     Input Parameters:
558: +   mat - location to store Jacobian
559: .   coloring - coloring context created with MatFDColoringCreate()
560: .   x1 - location at which Jacobian is to be computed
561: -   sctx - optional context required by function (actually a SNES context)

563:     Options Database Keys:
564: +    -mat_fd_coloring_freq <freq> - Sets coloring frequency
565: .    -mat_fd_type - "wp" or "ds"  (see MATSNESMF_WP or MATSNESMF_DS)
566: .    -mat_fd_coloring_view - Activates basic viewing or coloring
567: .    -mat_fd_coloring_view_draw - Activates drawing of coloring
568: -    -mat_fd_coloring_view_info - Activates viewing of coloring info

570:     Level: intermediate

572: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView()

574: .keywords: coloring, Jacobian, finite differences
575: @*/
576: PetscErrorCode  MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
577: {
578:   PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f;
580:   PetscInt       k,start,end,l,row,col,srow,**vscaleforrow,m1,m2;
581:   PetscScalar    dx,*y,*xx,*w3_array;
582:   PetscScalar    *vscale_array;
583:   PetscReal      epsilon = coloring->error_rel,umin = coloring->umin,unorm;
584:   Vec            w1=coloring->w1,w2=coloring->w2,w3;
585:   void           *fctx = coloring->fctx;
586:   PetscTruth     flg;
587:   PetscInt       ctype=coloring->ctype,N,col_start=0,col_end=0;
588:   Vec            x1_tmp;

594:   if (!f) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");

596:   if (coloring->usersetsrecompute) {
597:     if (!coloring->recompute) {
598:       *flag = SAME_PRECONDITIONER;
599:       PetscInfo(J,"Skipping Jacobian, since user called MatFDColorSetRecompute()\n");
600:       return(0);
601:     } else {
602:       coloring->recompute = PETSC_FALSE;
603:     }
604:   }

607:   MatSetUnfactored(J);
608:   PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg);
609:   if (flg) {
610:     PetscInfo(coloring,"Not calling MatZeroEntries()\n");
611:   } else {
612:     PetscTruth assembled;
613:     MatAssembled(J,&assembled);
614:     if (assembled) {
615:       MatZeroEntries(J);
616:     }
617:   }

619:   x1_tmp = x1;
620:   if (!coloring->vscale){
621:     VecDuplicate(x1_tmp,&coloring->vscale);
622:   }
623: 
624:   /*
625:     This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets
626:     coloring->F for the coarser grids from the finest
627:   */
628:   if (coloring->F) {
629:     VecGetLocalSize(coloring->F,&m1);
630:     VecGetLocalSize(w1,&m2);
631:     if (m1 != m2) {
632:       coloring->F = 0;
633:       }
634:     }

636:   if (coloring->htype[0] == 'w') { /* tacky test; need to make systematic if we add other approaches to computing h*/
637:     VecNorm(x1_tmp,NORM_2,&unorm);
638:   }
639:   VecGetOwnershipRange(w1,&start,&end); /* OwnershipRange is used by ghosted x! */
640: 
641:   /* Set w1 = F(x1) */
642:   if (coloring->F) {
643:     w1          = coloring->F; /* use already computed value of function */
644:     coloring->F = 0;
645:   } else {
647:     (*f)(sctx,x1_tmp,w1,fctx);
649:   }
650: 
651:   if (!coloring->w3) {
652:     VecDuplicate(x1_tmp,&coloring->w3);
653:     PetscLogObjectParent(coloring,coloring->w3);
654:   }
655:   w3 = coloring->w3;

657:     /* Compute all the local scale factors, including ghost points */
658:   VecGetLocalSize(x1_tmp,&N);
659:   VecGetArray(x1_tmp,&xx);
660:   VecGetArray(coloring->vscale,&vscale_array);
661:   if (ctype == IS_COLORING_GHOSTED){
662:     col_start = 0; col_end = N;
663:   } else if (ctype == IS_COLORING_GLOBAL){
664:     xx = xx - start;
665:     vscale_array = vscale_array - start;
666:     col_start = start; col_end = N + start;
667:   }
668:   for (col=col_start; col<col_end; col++){
669:     /* Loop over each local column, vscale[col] = 1./(epsilon*dx[col]) */
670:     if (coloring->htype[0] == 'w') {
671:       dx = 1.0 + unorm;
672:     } else {
673:       dx  = xx[col];
674:     }
675:     if (dx == 0.0) dx = 1.0;
676: #if !defined(PETSC_USE_COMPLEX)
677:     if (dx < umin && dx >= 0.0)      dx = umin;
678:     else if (dx < 0.0 && dx > -umin) dx = -umin;
679: #else
680:     if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
681:     else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
682: #endif
683:     dx               *= epsilon;
684:     vscale_array[col] = 1.0/dx;
685:   }
686:   if (ctype == IS_COLORING_GLOBAL)  vscale_array = vscale_array + start;
687:   VecRestoreArray(coloring->vscale,&vscale_array);
688:   if (ctype == IS_COLORING_GLOBAL){
689:     VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
690:     VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
691:   }
692: 
693:   if (coloring->vscaleforrow) {
694:     vscaleforrow = coloring->vscaleforrow;
695:   } else {
696:     SETERRQ(PETSC_ERR_ARG_NULL,"Null Object: coloring->vscaleforrow");
697:   }

699:   /*
700:     Loop over each color
701:   */
702:   VecGetArray(coloring->vscale,&vscale_array);
703:   for (k=0; k<coloring->ncolors; k++) {
704:     coloring->currentcolor = k;
705:     VecCopy(x1_tmp,w3);
706:     VecGetArray(w3,&w3_array);
707:     if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array - start;
708:     /*
709:       Loop over each column associated with color 
710:       adding the perturbation to the vector w3.
711:     */
712:     for (l=0; l<coloring->ncolumns[k]; l++) {
713:       col = coloring->columns[k][l];    /* local column of the matrix we are probing for */
714:       if (coloring->htype[0] == 'w') {
715:         dx = 1.0 + unorm;
716:       } else {
717:         dx  = xx[col];
718:       }
719:       if (dx == 0.0) dx = 1.0;
720: #if !defined(PETSC_USE_COMPLEX)
721:       if (dx < umin && dx >= 0.0)      dx = umin;
722:       else if (dx < 0.0 && dx > -umin) dx = -umin;
723: #else
724:       if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
725:       else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
726: #endif
727:       dx            *= epsilon;
728:       if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter");
729:       w3_array[col] += dx;
730:     }
731:     if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array + start;
732:     VecRestoreArray(w3,&w3_array);

734:     /*
735:       Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
736:                            w2 = F(x1 + dx) - F(x1)
737:     */
739:     (*f)(sctx,w3,w2,fctx);
741:     VecAXPY(w2,-1.0,w1);
742: 
743:     /*
744:       Loop over rows of vector, putting results into Jacobian matrix
745:     */
746:     VecGetArray(w2,&y);
747:     for (l=0; l<coloring->nrows[k]; l++) {
748:       row    = coloring->rows[k][l];             /* local row index */
749:       col    = coloring->columnsforrow[k][l];    /* global column index */
750:       y[row] *= vscale_array[vscaleforrow[k][l]];
751:       srow   = row + start;
752:       MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);
753:     }
754:     VecRestoreArray(w2,&y);
755:   } /* endof for each color */
756:   if (ctype == IS_COLORING_GLOBAL) xx = xx + start;
757:   VecRestoreArray(coloring->vscale,&vscale_array);
758:   VecRestoreArray(x1_tmp,&xx);
759: 
760:   coloring->currentcolor = -1;
761:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
762:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

765:   PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test",&flg);
766:   if (flg) {
767:     MatNullSpaceTest(J->nullsp,J);
768:   }
769:   MatFDColoringView_Private(coloring);
770:   return(0);
771: }

775: /*@
776:     MatFDColoringApplyTS - Given a matrix for which a MatFDColoring context 
777:     has been created, computes the Jacobian for a function via finite differences.

779:    Collective on Mat, MatFDColoring, and Vec

781:     Input Parameters:
782: +   mat - location to store Jacobian
783: .   coloring - coloring context created with MatFDColoringCreate()
784: .   x1 - location at which Jacobian is to be computed
785: -   sctx - optional context required by function (actually a SNES context)

787:    Options Database Keys:
788: .  -mat_fd_coloring_freq <freq> - Sets coloring frequency

790:    Level: intermediate

792: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView()

794: .keywords: coloring, Jacobian, finite differences
795: @*/
796: PetscErrorCode  MatFDColoringApplyTS(Mat J,MatFDColoring coloring,PetscReal t,Vec x1,MatStructure *flag,void *sctx)
797: {
798:   PetscErrorCode (*f)(void*,PetscReal,Vec,Vec,void*)=(PetscErrorCode (*)(void*,PetscReal,Vec,Vec,void *))coloring->f;
800:   PetscInt       k,N,start,end,l,row,col,srow,**vscaleforrow;
801:   PetscScalar    dx,*y,*xx,*w3_array;
802:   PetscScalar    *vscale_array;
803:   PetscReal      epsilon = coloring->error_rel,umin = coloring->umin;
804:   Vec            w1=coloring->w1,w2=coloring->w2,w3;
805:   void           *fctx = coloring->fctx;
806:   PetscTruth     flg;


814:   if (!coloring->w3) {
815:     VecDuplicate(x1,&coloring->w3);
816:     PetscLogObjectParent(coloring,coloring->w3);
817:   }
818:   w3 = coloring->w3;

820:   MatSetUnfactored(J);
821:   PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg);
822:   if (flg) {
823:     PetscInfo(coloring,"Not calling MatZeroEntries()\n");
824:   } else {
825:     PetscTruth assembled;
826:     MatAssembled(J,&assembled);
827:     if (assembled) {
828:       MatZeroEntries(J);
829:     }
830:   }

832:   VecGetOwnershipRange(x1,&start,&end);
833:   VecGetSize(x1,&N);
835:   (*f)(sctx,t,x1,w1,fctx);

838:   /* 
839:       Compute all the scale factors and share with other processors
840:   */
841:   VecGetArray(x1,&xx);xx = xx - start;
842:   VecGetArray(coloring->vscale,&vscale_array);vscale_array = vscale_array - start;
843:   for (k=0; k<coloring->ncolors; k++) {
844:     /*
845:        Loop over each column associated with color adding the 
846:        perturbation to the vector w3.
847:     */
848:     for (l=0; l<coloring->ncolumns[k]; l++) {
849:       col = coloring->columns[k][l];    /* column of the matrix we are probing for */
850:       dx  = xx[col];
851:       if (dx == 0.0) dx = 1.0;
852: #if !defined(PETSC_USE_COMPLEX)
853:       if (dx < umin && dx >= 0.0)      dx = umin;
854:       else if (dx < 0.0 && dx > -umin) dx = -umin;
855: #else
856:       if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
857:       else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
858: #endif
859:       dx                *= epsilon;
860:       vscale_array[col] = 1.0/dx;
861:     }
862:   }
863:   vscale_array = vscale_array - start;VecRestoreArray(coloring->vscale,&vscale_array);
864:   VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
865:   VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);

867:   if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow;
868:   else                        vscaleforrow = coloring->columnsforrow;

870:   VecGetArray(coloring->vscale,&vscale_array);
871:   /*
872:       Loop over each color
873:   */
874:   for (k=0; k<coloring->ncolors; k++) {
875:     VecCopy(x1,w3);
876:     VecGetArray(w3,&w3_array);w3_array = w3_array - start;
877:     /*
878:        Loop over each column associated with color adding the 
879:        perturbation to the vector w3.
880:     */
881:     for (l=0; l<coloring->ncolumns[k]; l++) {
882:       col = coloring->columns[k][l];    /* column of the matrix we are probing for */
883:       dx  = xx[col];
884:       if (dx == 0.0) dx = 1.0;
885: #if !defined(PETSC_USE_COMPLEX)
886:       if (dx < umin && dx >= 0.0)      dx = umin;
887:       else if (dx < 0.0 && dx > -umin) dx = -umin;
888: #else
889:       if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
890:       else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
891: #endif
892:       dx            *= epsilon;
893:       w3_array[col] += dx;
894:     }
895:     w3_array = w3_array + start; VecRestoreArray(w3,&w3_array);

897:     /*
898:        Evaluate function at x1 + dx (here dx is a vector of perturbations)
899:     */
901:     (*f)(sctx,t,w3,w2,fctx);
903:     VecAXPY(w2,-1.0,w1);

905:     /*
906:        Loop over rows of vector, putting results into Jacobian matrix
907:     */
908:     VecGetArray(w2,&y);
909:     for (l=0; l<coloring->nrows[k]; l++) {
910:       row    = coloring->rows[k][l];
911:       col    = coloring->columnsforrow[k][l];
912:       y[row] *= vscale_array[vscaleforrow[k][l]];
913:       srow   = row + start;
914:       MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);
915:     }
916:     VecRestoreArray(w2,&y);
917:   }
918:   VecRestoreArray(coloring->vscale,&vscale_array);
919:   xx    = xx + start; VecRestoreArray(x1,&xx);
920:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
921:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
923:   return(0);
924: }


929: /*@C
930:    MatFDColoringSetRecompute - Indicates that the next time a Jacobian preconditioner
931:      is needed it sholuld be recomputed. Once this is called and the new Jacobian is computed
932:      no additional Jacobian's will be computed (the same one will be used) until this is
933:      called again.

935:    Collective on MatFDColoring

937:    Input  Parameters:
938: .  c - the coloring context

940:    Level: intermediate

942:    Notes: The MatFDColoringSetFrequency() is ignored once this is called

944: .seealso: MatFDColoringCreate(), MatFDColoringSetFrequency()

946: .keywords: Mat, finite differences, coloring
947: @*/
948: PetscErrorCode  MatFDColoringSetRecompute(MatFDColoring c)
949: {
952:   c->usersetsrecompute = PETSC_TRUE;
953:   c->recompute         = PETSC_TRUE;
954:   return(0);
955: }