Actual source code: snesj.c

      1: #define PETSCSNES_DLL

     3:  #include src/snes/snesimpl.h

      7: /*@C
      8:    SNESDefaultComputeJacobian - Computes the Jacobian using finite differences.

     10:    Collective on SNES

     12:    Input Parameters:
     13: +  x1 - compute Jacobian at this point
     14: -  ctx - application's function context, as set with SNESSetFunction()

     16:    Output Parameters:
     17: +  J - Jacobian matrix (not altered in this routine)
     18: .  B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)
     19: -  flag - flag indicating whether the matrix sparsity structure has changed

     21:    Options Database Key:
     22: +  -snes_fd - Activates SNESDefaultComputeJacobian()
     23: -  -snes_test_err - Square root of function error tolerance, default square root of machine
     24:                     epsilon (1.e-8 in double, 3.e-4 in single)

     26:    Notes:
     27:    This routine is slow and expensive, and is not currently optimized
     28:    to take advantage of sparsity in the problem.  Although
     29:    SNESDefaultComputeJacobian() is not recommended for general use
     30:    in large-scale applications, It can be useful in checking the
     31:    correctness of a user-provided Jacobian.

     33:    An alternative routine that uses coloring to explot matrix sparsity is
     34:    SNESDefaultComputeJacobianColor().

     36:    Level: intermediate

     38: .keywords: SNES, finite differences, Jacobian

     40: .seealso: SNESSetJacobian(), SNESDefaultComputeJacobianColor()
     41: @*/
     42: PetscErrorCode PETSCSNES_DLLEXPORT SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
     43: {
     44:   Vec            j1a,j2a,x2;
     46:   PetscInt       i,N,start,end,j;
     47:   PetscScalar    dx,mone = -1.0,*y,scale,*xx,wscale;
     48:   PetscReal      amax,epsilon = PETSC_SQRT_MACHINE_EPSILON;
     49:   PetscReal      dx_min = 1.e-16,dx_par = 1.e-1;
     50:   MPI_Comm       comm;
     51:   PetscErrorCode (*eval_fct)(SNES,Vec,Vec)=0;
     52:   PetscTruth     assembled;

     55:   PetscOptionsGetReal(snes->prefix,"-snes_test_err",&epsilon,0);
     56:   eval_fct = SNESComputeFunction;

     58:   PetscObjectGetComm((PetscObject)x1,&comm);
     59:   MatAssembled(*B,&assembled);
     60:   if (assembled) {
     61:     MatZeroEntries(*B);
     62:   }
     63:   if (!snes->nvwork) {
     64:     VecDuplicateVecs(x1,3,&snes->vwork);
     65:     snes->nvwork = 3;
     66:     PetscLogObjectParents(snes,3,snes->vwork);
     67:   }
     68:   j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2];

     70:   VecGetSize(x1,&N);
     71:   VecGetOwnershipRange(x1,&start,&end);
     72:   (*eval_fct)(snes,x1,j1a);

     74:   /* Compute Jacobian approximation, 1 column at a time.
     75:       x1 = current iterate, j1a = F(x1)
     76:       x2 = perturbed iterate, j2a = F(x2)
     77:    */
     78:   for (i=0; i<N; i++) {
     79:     VecCopy(x1,x2);
     80:     if (i>= start && i<end) {
     81:       VecGetArray(x1,&xx);
     82:       dx = xx[i-start];
     83:       VecRestoreArray(x1,&xx);
     84: #if !defined(PETSC_USE_COMPLEX)
     85:       if (dx < dx_min && dx >= 0.0) dx = dx_par;
     86:       else if (dx < 0.0 && dx > -dx_min) dx = -dx_par;
     87: #else
     88:       if (PetscAbsScalar(dx) < dx_min && PetscRealPart(dx) >= 0.0) dx = dx_par;
     89:       else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < dx_min) dx = -dx_par;
     90: #endif
     91:       dx *= epsilon;
     92:       wscale = 1.0/dx;
     93:       VecSetValues(x2,1,&i,&dx,ADD_VALUES);
     94:     } else {
     95:       wscale = 0.0;
     96:     }
     97:     (*eval_fct)(snes,x2,j2a);
     98:     VecAXPY(j2a,mone,j1a);
     99:     /* Communicate scale to all processors */
    100:     MPI_Allreduce(&wscale,&scale,1,MPIU_SCALAR,PetscSum_Op,comm);
    101:     VecScale(j2a,scale);
    102:     VecNorm(j2a,NORM_INFINITY,&amax); amax *= 1.e-14;
    103:     VecGetArray(j2a,&y);
    104:     for (j=start; j<end; j++) {
    105:       if (PetscAbsScalar(y[j-start]) > amax) {
    106:         MatSetValues(*B,1,&j,1,&i,y+j-start,INSERT_VALUES);
    107:       }
    108:     }
    109:     VecRestoreArray(j2a,&y);
    110:   }
    111:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
    112:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
    113:   if (*B != *J) {
    114:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    115:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    116:   }
    117:   *flag =  DIFFERENT_NONZERO_PATTERN;
    118:   return(0);
    119: }