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 src/mat/matimpl.h

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

     15: PetscErrorCode PETSCMAT_DLLEXPORT 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 PETSCMAT_DLLEXPORT MatFDColoringView(MatFDColoring c,PetscViewer viewer)
     99: {
    100:   PetscErrorCode    ierr;
    101:   PetscInt          i,j;
    102:   PetscTruth        isdraw,iascii;
    103:   PetscViewerFormat format;

    107:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(c->comm);

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

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

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

    148:    Collective on MatFDColoring

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

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

    163:    Level: advanced

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

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

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

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

    185:    Collective on MatFDColoring

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

    191:    Options Database Keys:
    192: .  -mat_fd_coloring_freq <freq>  - Sets coloring frequency

    194:    Level: advanced

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

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

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

    211:   matfd->freq = freq;
    212:   return(0);
    213: }

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

    221:    Not Collective

    223:    Input Parameters:
    224: .  coloring - the coloring context

    226:    Output Parameters:
    227: .  freq - frequency (default is 1)

    229:    Options Database Keys:
    230: .  -mat_fd_coloring_freq <freq> - Sets coloring frequency

    232:    Level: advanced

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

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

    242: .seealso: MatFDColoringSetFrequency()
    243: @*/
    244: PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringGetFrequency(MatFDColoring matfd,PetscInt *freq)
    245: {
    248:   *freq = matfd->freq;
    249:   return(0);
    250: }

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

    257:    Collective on MatFDColoring

    259:    Input Parameters:
    260: +  coloring - the coloring context
    261: .  f - the function
    262: -  fctx - the optional user-defined function context

    264:    Level: intermediate

    266:    Notes:
    267:     In Fortran you must call MatFDColoringSetFunctionSNES() for a coloring object to
    268:   be used with the SNES solvers and MatFDColoringSetFunctionTS() if it is to be used
    269:   with the TS solvers.

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

    284: /*@
    285:    MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
    286:    the options database.

    288:    Collective on MatFDColoring

    290:    The Jacobian, F'(u), is estimated with the differencing approximation
    291: .vb
    292:        F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
    293:        h = error_rel*u[i]                 if  abs(u[i]) > umin
    294:          = +/- error_rel*umin             otherwise, with +/- determined by the sign of u[i]
    295:        dx_{i} = (0, ... 1, .... 0)
    296: .ve

    298:    Input Parameter:
    299: .  coloring - the coloring context

    301:    Options Database Keys:
    302: +  -mat_fd_coloring_err <err> - Sets <err> (square root
    303:            of relative error in the function)
    304: .  -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
    305: .  -mat_fd_coloring_freq <freq> - Sets frequency of computing a new Jacobian
    306: .  -mat_fd_coloring_view - Activates basic viewing
    307: .  -mat_fd_coloring_view_info - Activates viewing info
    308: -  -mat_fd_coloring_view_draw - Activates drawing

    310:     Level: intermediate

    312: .keywords: Mat, finite differences, parameters

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

    316: @*/
    317: PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringSetFromOptions(MatFDColoring matfd)
    318: {

    324:   PetscOptionsBegin(matfd->comm,matfd->prefix,"Jacobian computation via finite differences option","MatFD");
    325:     PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
    326:     PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
    327:     PetscOptionsInt("-mat_fd_coloring_freq","How often Jacobian is recomputed","MatFDColoringSetFrequency",matfd->freq,&matfd->freq,0);
    328:     /* not used here; but so they are presented in the GUI */
    329:     PetscOptionsName("-mat_fd_coloring_view","Print entire datastructure used for Jacobian","None",0);
    330:     PetscOptionsName("-mat_fd_coloring_view_info","Print number of colors etc for Jacobian","None",0);
    331:     PetscOptionsName("-mat_fd_coloring_view_draw","Plot nonzero structure ofJacobian","None",0);
    332:   PetscOptionsEnd();
    333:   return(0);
    334: }

    338: PetscErrorCode MatFDColoringView_Private(MatFDColoring fd)
    339: {
    341:   PetscTruth     flg;

    344:   PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_view",&flg);
    345:   if (flg) {
    346:     MatFDColoringView(fd,PETSC_VIEWER_STDOUT_(fd->comm));
    347:   }
    348:   PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_info",&flg);
    349:   if (flg) {
    350:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(fd->comm),PETSC_VIEWER_ASCII_INFO);
    351:     MatFDColoringView(fd,PETSC_VIEWER_STDOUT_(fd->comm));
    352:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(fd->comm));
    353:   }
    354:   PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_draw",&flg);
    355:   if (flg) {
    356:     MatFDColoringView(fd,PETSC_VIEWER_DRAW_(fd->comm));
    357:     PetscViewerFlush(PETSC_VIEWER_DRAW_(fd->comm));
    358:   }
    359:   return(0);
    360: }

    364: /*@C
    365:    MatFDColoringCreate - Creates a matrix coloring context for finite difference
    366:    computation of Jacobians.

    368:    Collective on Mat

    370:    Input Parameters:
    371: +  mat - the matrix containing the nonzero structure of the Jacobian
    372: -  iscoloring - the coloring of the matrix

    374:     Output Parameter:
    375: .   color - the new coloring context
    376:
    377:     Level: intermediate

    379: .seealso: MatFDColoringDestroy(),SNESDefaultComputeJacobianColor(), ISColoringCreate(),
    380:           MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
    381:           MatFDColoringSetFrequency(), MatFDColoringSetRecompute(), MatFDColoringView(),
    382:           MatFDColoringSetParameters()
    383: @*/
    384: PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
    385: {
    386:   MatFDColoring  c;
    387:   MPI_Comm       comm;
    389:   PetscInt       M,N;

    392:   PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
    393:   MatGetSize(mat,&M,&N);
    394:   if (M != N) SETERRQ(PETSC_ERR_SUP,"Only for square matrices");

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

    399:   if (mat->ops->fdcoloringcreate) {
    400:     (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
    401:   } else {
    402:     SETERRQ(PETSC_ERR_SUP,"Code not yet written for this matrix type");
    403:   }

    405:   c->error_rel         = PETSC_SQRT_MACHINE_EPSILON;
    406:   c->umin              = 100.0*PETSC_SQRT_MACHINE_EPSILON;
    407:   c->freq              = 1;
    408:   c->usersetsrecompute = PETSC_FALSE;
    409:   c->recompute         = PETSC_FALSE;
    410:   c->currentcolor      = -1;

    412:   *color = c;
    413:   PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
    414:   return(0);
    415: }

    419: /*@C
    420:     MatFDColoringDestroy - Destroys a matrix coloring context that was created
    421:     via MatFDColoringCreate().

    423:     Collective on MatFDColoring

    425:     Input Parameter:
    426: .   c - coloring context

    428:     Level: intermediate

    430: .seealso: MatFDColoringCreate()
    431: @*/
    432: PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringDestroy(MatFDColoring c)
    433: {
    435:   PetscInt       i;

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

    440:   for (i=0; i<c->ncolors; i++) {
    441:     if (c->columns[i])         {PetscFree(c->columns[i]);}
    442:     if (c->rows[i])            {PetscFree(c->rows[i]);}
    443:     if (c->columnsforrow[i])   {PetscFree(c->columnsforrow[i]);}
    444:     if (c->vscaleforrow && c->vscaleforrow[i]) {PetscFree(c->vscaleforrow[i]);}
    445:   }
    446:   PetscFree(c->ncolumns);
    447:   PetscFree(c->columns);
    448:   PetscFree(c->nrows);
    449:   PetscFree(c->rows);
    450:   PetscFree(c->columnsforrow);
    451:   if (c->vscaleforrow) {PetscFree(c->vscaleforrow);}
    452:   if (c->vscale)       {VecDestroy(c->vscale);}
    453:   if (c->w1) {
    454:     VecDestroy(c->w1);
    455:     VecDestroy(c->w2);
    456:     VecDestroy(c->w3);
    457:   }
    458:   PetscHeaderDestroy(c);
    459:   return(0);
    460: }

    464: /*@C
    465:     MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
    466:       that are currently being perturbed.

    468:     Not Collective

    470:     Input Parameters:
    471: .   coloring - coloring context created with MatFDColoringCreate()

    473:     Output Parameters:
    474: +   n - the number of local columns being perturbed
    475: -   cols - the column indices, in global numbering

    477:    Level: intermediate

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

    481: .keywords: coloring, Jacobian, finite differences
    482: @*/
    483: PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,PetscInt *cols[])
    484: {
    486:   if (coloring->currentcolor >= 0) {
    487:     *n    = coloring->ncolumns[coloring->currentcolor];
    488:     *cols = coloring->columns[coloring->currentcolor];
    489:   } else {
    490:     *n = 0;
    491:   }
    492:   return(0);
    493: }

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

    501:     Collective on MatFDColoring

    503:     Input Parameters:
    504: +   mat - location to store Jacobian
    505: .   coloring - coloring context created with MatFDColoringCreate()
    506: .   x1 - location at which Jacobian is to be computed
    507: -   sctx - optional context required by function (actually a SNES context)

    509:     Options Database Keys:
    510: +    -mat_fd_coloring_freq <freq> - Sets coloring frequency
    511: .    -mat_fd_coloring_view - Activates basic viewing or coloring
    512: .    -mat_fd_coloring_view_draw - Activates drawing of coloring
    513: -    -mat_fd_coloring_view_info - Activates viewing of coloring info

    515:     Level: intermediate

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

    519: .keywords: coloring, Jacobian, finite differences
    520: @*/
    521: PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
    522: {
    523:   PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f;
    525:   PetscInt       k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2;
    526:   PetscScalar    dx,mone = -1.0,*y,*xx,*w3_array;
    527:   PetscScalar    *vscale_array;
    528:   PetscReal      epsilon = coloring->error_rel,umin = coloring->umin;
    529:   Vec            w1,w2,w3;
    530:   void           *fctx = coloring->fctx;
    531:   PetscTruth     flg;

    539:   if (coloring->usersetsrecompute) {
    540:     if (!coloring->recompute) {
    541:       *flag = SAME_PRECONDITIONER;
    542:       PetscLogInfo((J,"MatFDColoringApply:Skipping Jacobian, since user called MatFDColorSetRecompute()\n"));
    543:       return(0);
    544:     } else {
    545:       coloring->recompute = PETSC_FALSE;
    546:     }
    547:   }

    549:   PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
    550:   if (J->ops->fdcoloringapply) {
    551:     (*J->ops->fdcoloringapply)(J,coloring,x1,flag,sctx);
    552:   } else {

    554:     if (!coloring->w1) {
    555:       VecDuplicate(x1,&coloring->w1);
    556:       PetscLogObjectParent(coloring,coloring->w1);
    557:       VecDuplicate(x1,&coloring->w2);
    558:       PetscLogObjectParent(coloring,coloring->w2);
    559:       VecDuplicate(x1,&coloring->w3);
    560:       PetscLogObjectParent(coloring,coloring->w3);
    561:     }
    562:     w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3;

    564:     MatSetUnfactored(J);
    565:     PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg);
    566:     if (flg) {
    567:       PetscLogInfo((coloring,"MatFDColoringApply: Not calling MatZeroEntries()\n"));
    568:     } else {
    569:       PetscTruth assembled;
    570:       MatAssembled(J,&assembled);
    571:       if (assembled) {
    572:         MatZeroEntries(J);
    573:       }
    574:     }

    576:     VecGetOwnershipRange(x1,&start,&end);
    577:     VecGetSize(x1,&N);
    578:
    579:     /*
    580:       This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets
    581:       coloring->F for the coarser grids from the finest
    582:     */
    583:     if (coloring->F) {
    584:       VecGetLocalSize(coloring->F,&m1);
    585:       VecGetLocalSize(w1,&m2);
    586:       if (m1 != m2) {
    587:         coloring->F = 0;
    588:       }
    589:     }

    591:     if (coloring->F) {
    592:       w1          = coloring->F; /* use already computed value of function */
    593:       coloring->F = 0;
    594:     } else {
    595:       PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
    596:       (*f)(sctx,x1,w1,fctx);
    597:       PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
    598:     }

    600:     /*
    601:        Compute all the scale factors and share with other processors
    602:     */
    603:     VecGetArray(x1,&xx);xx = xx - start;
    604:     VecGetArray(coloring->vscale,&vscale_array);vscale_array = vscale_array - start;
    605:     for (k=0; k<coloring->ncolors; k++) {
    606:       /*
    607:         Loop over each column associated with color adding the
    608:         perturbation to the vector w3.
    609:       */
    610:       for (l=0; l<coloring->ncolumns[k]; l++) {
    611:         col = coloring->columns[k][l];    /* column of the matrix we are probing for */
    612:         dx  = xx[col];
    613:         if (dx == 0.0) dx = 1.0;
    614: #if !defined(PETSC_USE_COMPLEX)
    615:         if (dx < umin && dx >= 0.0)      dx = umin;
    616:         else if (dx < 0.0 && dx > -umin) dx = -umin;
    617: #else
    618:         if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
    619:         else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
    620: #endif
    621:         dx                *= epsilon;
    622:         vscale_array[col] = 1.0/dx;
    623:       }
    624:     }
    625:     vscale_array = vscale_array + start;VecRestoreArray(coloring->vscale,&vscale_array);
    626:     VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
    627:     VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);

    629:     /*  VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD);
    630:         VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/

    632:     if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow;
    633:     else                        vscaleforrow = coloring->columnsforrow;

    635:     VecGetArray(coloring->vscale,&vscale_array);
    636:     /*
    637:       Loop over each color
    638:     */
    639:     for (k=0; k<coloring->ncolors; k++) {
    640:       coloring->currentcolor = k;
    641:       VecCopy(x1,w3);
    642:       VecGetArray(w3,&w3_array);w3_array = w3_array - start;
    643:       /*
    644:         Loop over each column associated with color adding the
    645:         perturbation to the vector w3.
    646:       */
    647:       for (l=0; l<coloring->ncolumns[k]; l++) {
    648:         col = coloring->columns[k][l];    /* column of the matrix we are probing for */
    649:         dx  = xx[col];
    650:         if (dx == 0.0) dx = 1.0;
    651: #if !defined(PETSC_USE_COMPLEX)
    652:         if (dx < umin && dx >= 0.0)      dx = umin;
    653:         else if (dx < 0.0 && dx > -umin) dx = -umin;
    654: #else
    655:         if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
    656:         else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
    657: #endif
    658:         dx            *= epsilon;
    659:         if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter");
    660:         w3_array[col] += dx;
    661:       }
    662:       w3_array = w3_array + start; VecRestoreArray(w3,&w3_array);

    664:       /*
    665:         Evaluate function at x1 + dx (here dx is a vector of perturbations)
    666:       */

    668:       PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
    669:       (*f)(sctx,w3,w2,fctx);
    670:       PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
    671:       VecAXPY(w2,mone,w1);

    673:       /*
    674:         Loop over rows of vector, putting results into Jacobian matrix
    675:       */
    676:       VecGetArray(w2,&y);
    677:       for (l=0; l<coloring->nrows[k]; l++) {
    678:         row    = coloring->rows[k][l];
    679:         col    = coloring->columnsforrow[k][l];
    680:         y[row] *= vscale_array[vscaleforrow[k][l]];
    681:         srow   = row + start;
    682:         MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);
    683:       }
    684:       VecRestoreArray(w2,&y);
    685:     }
    686:     coloring->currentcolor = -1;
    687:     VecRestoreArray(coloring->vscale,&vscale_array);
    688:     xx = xx + start; VecRestoreArray(x1,&xx);
    689:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    690:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    691:   }
    692:   PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);

    694:   PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test",&flg);
    695:   if (flg) {
    696:     MatNullSpaceTest(J->nullsp,J);
    697:   }
    698:   MatFDColoringView_Private(coloring);

    700:   return(0);
    701: }

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

    709:    Collective on Mat, MatFDColoring, and Vec

    711:     Input Parameters:
    712: +   mat - location to store Jacobian
    713: .   coloring - coloring context created with MatFDColoringCreate()
    714: .   x1 - location at which Jacobian is to be computed
    715: -   sctx - optional context required by function (actually a SNES context)

    717:    Options Database Keys:
    718: .  -mat_fd_coloring_freq <freq> - Sets coloring frequency

    720:    Level: intermediate

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

    724: .keywords: coloring, Jacobian, finite differences
    725: @*/
    726: PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApplyTS(Mat J,MatFDColoring coloring,PetscReal t,Vec x1,MatStructure *flag,void *sctx)
    727: {
    728:   PetscErrorCode (*f)(void*,PetscReal,Vec,Vec,void*)=(PetscErrorCode (*)(void*,PetscReal,Vec,Vec,void *))coloring->f;
    730:   PetscInt       k,N,start,end,l,row,col,srow,**vscaleforrow;
    731:   PetscScalar    dx,mone = -1.0,*y,*xx,*w3_array;
    732:   PetscScalar    *vscale_array;
    733:   PetscReal      epsilon = coloring->error_rel,umin = coloring->umin;
    734:   Vec            w1,w2,w3;
    735:   void           *fctx = coloring->fctx;
    736:   PetscTruth     flg;

    743:   PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
    744:   if (!coloring->w1) {
    745:     VecDuplicate(x1,&coloring->w1);
    746:     PetscLogObjectParent(coloring,coloring->w1);
    747:     VecDuplicate(x1,&coloring->w2);
    748:     PetscLogObjectParent(coloring,coloring->w2);
    749:     VecDuplicate(x1,&coloring->w3);
    750:     PetscLogObjectParent(coloring,coloring->w3);
    751:   }
    752:   w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3;

    754:   MatSetUnfactored(J);
    755:   PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg);
    756:   if (flg) {
    757:     PetscLogInfo((coloring,"MatFDColoringApply: Not calling MatZeroEntries()\n"));
    758:   } else {
    759:     PetscTruth assembled;
    760:     MatAssembled(J,&assembled);
    761:     if (assembled) {
    762:       MatZeroEntries(J);
    763:     }
    764:   }

    766:   VecGetOwnershipRange(x1,&start,&end);
    767:   VecGetSize(x1,&N);
    768:   PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
    769:   (*f)(sctx,t,x1,w1,fctx);
    770:   PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);

    772:   /*
    773:       Compute all the scale factors and share with other processors
    774:   */
    775:   VecGetArray(x1,&xx);xx = xx - start;
    776:   VecGetArray(coloring->vscale,&vscale_array);vscale_array = vscale_array - start;
    777:   for (k=0; k<coloring->ncolors; k++) {
    778:     /*
    779:        Loop over each column associated with color adding the
    780:        perturbation to the vector w3.
    781:     */
    782:     for (l=0; l<coloring->ncolumns[k]; l++) {
    783:       col = coloring->columns[k][l];    /* column of the matrix we are probing for */
    784:       dx  = xx[col];
    785:       if (dx == 0.0) dx = 1.0;
    786: #if !defined(PETSC_USE_COMPLEX)
    787:       if (dx < umin && dx >= 0.0)      dx = umin;
    788:       else if (dx < 0.0 && dx > -umin) dx = -umin;
    789: #else
    790:       if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
    791:       else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
    792: #endif
    793:       dx                *= epsilon;
    794:       vscale_array[col] = 1.0/dx;
    795:     }
    796:   }
    797:   vscale_array = vscale_array - start;VecRestoreArray(coloring->vscale,&vscale_array);
    798:   VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
    799:   VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);

    801:   if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow;
    802:   else                        vscaleforrow = coloring->columnsforrow;

    804:   VecGetArray(coloring->vscale,&vscale_array);
    805:   /*
    806:       Loop over each color
    807:   */
    808:   for (k=0; k<coloring->ncolors; k++) {
    809:     VecCopy(x1,w3);
    810:     VecGetArray(w3,&w3_array);w3_array = w3_array - start;
    811:     /*
    812:        Loop over each column associated with color adding the
    813:        perturbation to the vector w3.
    814:     */
    815:     for (l=0; l<coloring->ncolumns[k]; l++) {
    816:       col = coloring->columns[k][l];    /* column of the matrix we are probing for */
    817:       dx  = xx[col];
    818:       if (dx == 0.0) dx = 1.0;
    819: #if !defined(PETSC_USE_COMPLEX)
    820:       if (dx < umin && dx >= 0.0)      dx = umin;
    821:       else if (dx < 0.0 && dx > -umin) dx = -umin;
    822: #else
    823:       if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
    824:       else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
    825: #endif
    826:       dx            *= epsilon;
    827:       w3_array[col] += dx;
    828:     }
    829:     w3_array = w3_array + start; VecRestoreArray(w3,&w3_array);

    831:     /*
    832:        Evaluate function at x1 + dx (here dx is a vector of perturbations)
    833:     */
    834:     PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
    835:     (*f)(sctx,t,w3,w2,fctx);
    836:     PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
    837:     VecAXPY(w2,mone,w1);

    839:     /*
    840:        Loop over rows of vector, putting results into Jacobian matrix
    841:     */
    842:     VecGetArray(w2,&y);
    843:     for (l=0; l<coloring->nrows[k]; l++) {
    844:       row    = coloring->rows[k][l];
    845:       col    = coloring->columnsforrow[k][l];
    846:       y[row] *= vscale_array[vscaleforrow[k][l]];
    847:       srow   = row + start;
    848:       MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);
    849:     }
    850:     VecRestoreArray(w2,&y);
    851:   }
    852:   VecRestoreArray(coloring->vscale,&vscale_array);
    853:   xx    = xx + start; VecRestoreArray(x1,&xx);
    854:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    855:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    856:   PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);
    857:   return(0);
    858: }

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

    869:    Collective on MatFDColoring

    871:    Input  Parameters:
    872: .  c - the coloring context

    874:    Level: intermediate

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

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

    880: .keywords: Mat, finite differences, coloring
    881: @*/
    882: PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringSetRecompute(MatFDColoring c)
    883: {
    886:   c->usersetsrecompute = PETSC_TRUE;
    887:   c->recompute         = PETSC_TRUE;
    888:   return(0);
    889: }