2: #include <petsc-private/snesimpl.h>    /*I  "petscsnes.h"  I*/

  6: /*@C
  7:    SNESComputeJacobianDefault - Computes the Jacobian using finite differences.

  9:    Collective on SNES

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

 15:    Output Parameters:
 16: +  J - Jacobian matrix (not altered in this routine)
 17: -  B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)

 19:    Options Database Key:
 20: +  -snes_fd - Activates SNESComputeJacobianDefault()
 21: .  -snes_test_err - Square root of function error tolerance, default square root of machine
 22:                     epsilon (1.e-8 in double, 3.e-4 in single)
 23: -  -mat_fd_type - Either wp or ds (see MATMFFD_WP or MATMFFD_DS)

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

 32:    An alternative routine that uses coloring to exploit matrix sparsity is
 33:    SNESComputeJacobianDefaultColor().

 35:    Level: intermediate

 37: .keywords: SNES, finite differences, Jacobian

 39: .seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF()
 40: @*/
 41: PetscErrorCode  SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
 42: {
 43:   Vec            j1a,j2a,x2;
 45:   PetscInt       i,N,start,end,j,value,root;
 46:   PetscScalar    dx,*y,*xx,wscale;
 47:   PetscReal      amax,epsilon = PETSC_SQRT_MACHINE_EPSILON;
 48:   PetscReal      dx_min = 1.e-16,dx_par = 1.e-1,unorm;
 49:   MPI_Comm       comm;
 50:   PetscErrorCode (*eval_fct)(SNES,Vec,Vec)=0;
 51:   PetscBool      assembled,use_wp = PETSC_TRUE,flg;
 52:   const char     *list[2] = {"ds","wp"};
 53:   PetscMPIInt    size;
 54:   const PetscInt *ranges;

 57:   PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);
 58:   eval_fct = SNESComputeFunction;

 60:   PetscObjectGetComm((PetscObject)x1,&comm);
 61:   MPI_Comm_size(comm,&size);
 62:   MatAssembled(B,&assembled);
 63:   if (assembled) {
 64:     MatZeroEntries(B);
 65:   }
 66:   if (!snes->nvwork) {
 67:     snes->nvwork = 3;

 69:     VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);
 70:     PetscLogObjectParents(snes,snes->nvwork,snes->vwork);
 71:   }
 72:   j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2];

 74:   VecGetSize(x1,&N);
 75:   VecGetOwnershipRange(x1,&start,&end);
 76:   (*eval_fct)(snes,x1,j1a);

 78:   PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg);
 79:   if (flg && !value) use_wp = PETSC_FALSE;

 81:   if (use_wp) {
 82:     VecNorm(x1,NORM_2,&unorm);
 83:   }
 84:   /* Compute Jacobian approximation, 1 column at a time.
 85:       x1 = current iterate, j1a = F(x1)
 86:       x2 = perturbed iterate, j2a = F(x2)
 87:    */
 88:   for (i=0; i<N; i++) {
 89:     VecCopy(x1,x2);
 90:     if (i>= start && i<end) {
 91:       VecGetArray(x1,&xx);
 92:       if (use_wp) dx = 1.0 + unorm;
 93:       else        dx = xx[i-start];
 94:       VecRestoreArray(x1,&xx);
 95:       if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par;
 96:       dx    *= epsilon;
 97:       wscale = 1.0/dx;
 98:       VecSetValues(x2,1,&i,&dx,ADD_VALUES);
 99:     } else {
100:       wscale = 0.0;
101:     }
102:     VecAssemblyBegin(x2);
103:     VecAssemblyEnd(x2);
104:     (*eval_fct)(snes,x2,j2a);
105:     VecAXPY(j2a,-1.0,j1a);
106:     /* Communicate scale=1/dx_i to all processors */
107:     VecGetOwnershipRanges(x1,&ranges);
108:     root = size;
109:     for (j=size-1; j>-1; j--) {
110:       root--;
111:       if (i>=ranges[j]) break;
112:     }
113:     MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm);

115:     VecScale(j2a,wscale);
116:     VecNorm(j2a,NORM_INFINITY,&amax); amax *= 1.e-14;
117:     VecGetArray(j2a,&y);
118:     for (j=start; j<end; j++) {
119:       if (PetscAbsScalar(y[j-start]) > amax || j == i) {
120:         MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES);
121:       }
122:     }
123:     VecRestoreArray(j2a,&y);
124:   }
125:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
126:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
127:   if (B != J) {
128:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
129:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
130:   }
131:   return(0);
132: }