2: /*
  3:    This file contains utility routines for finite difference
  4:    approximation of Jacobian matrices.  This functionality for
  5:    the TS component will eventually be incorporated as part of
  6:    the base PETSc libraries.
  7: */
  8: #include <petsc-private/tsimpl.h>
  9: #include <petsc-private/snesimpl.h>
 10: #include <petsc-private/fortranimpl.h>

 12: PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
 13: PetscErrorCode RHSJacobianFD(TS,PetscReal,Vec,Mat*,Mat*,MatStructure *,void*);

 15: /* -------------------------------------------------------------------*/

 17: /* Temporary interface routine; this will be eliminated soon! */
 18: #ifdef PETSC_HAVE_FORTRAN_CAPS
 19: #define setcroutinefromfortran_ SETCROUTINEFROMFORTRAN
 20: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
 21: #define setcroutinefromfortran_ setcroutinefromfortran
 22: #endif

 24: EXTERN_C_BEGIN

 26: void PETSC_STDCALL setcroutinefromfortran_(TS *ts,Mat *A,Mat *B,PetscErrorCode *__ierr)
 27: {
 28:     *__TSSetRHSJacobian(*ts,*A,*B,RHSJacobianFD,PETSC_NULL);
 29: }

 31: EXTERN_C_END


 34: /* -------------------------------------------------------------------*/
 35: /*
 36:    RHSJacobianFD - Computes the Jacobian using finite differences.

 38:    Input Parameters:
 39: .  ts - TS context
 40: .  xx1 - compute Jacobian at this point
 41: .  ctx - application's function context, as set with SNESSetFunction()

 43:    Output Parameters:
 44: .  J - Jacobian
 45: .  B - preconditioner, same as Jacobian
 46: .  flag - matrix flag

 48:    Notes:
 49:    This routine is slow and expensive, and is not currently optimized
 50:    to take advantage of sparsity in the problem.

 52:    Sparse approximations using colorings are also available and
 53:    would be a much better alternative!
 54: */
 55: PetscErrorCode RHSJacobianFD(TS ts,PetscReal t,Vec xx1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
 56: {
 57:   Vec            jj1,jj2,xx2;
 58:   PetscInt       i,N,start,end,j;
 60:   PetscScalar    dx,*y,scale,*xx,wscale;
 61:   PetscReal      amax,epsilon = 1.e-8; /* assumes PetscReal precision */
 62:   PetscReal      dx_min = 1.e-16,dx_par = 1.e-1;
 63:   MPI_Comm       comm;
 64:   PetscBool      assembled;

 66:   VecDuplicate(xx1,&jj1);
 67:   VecDuplicate(xx1,&jj2);
 68:   VecDuplicate(xx1,&xx2);

 70:   PetscObjectGetComm((PetscObject)xx1,&comm);
 71:   MatAssembled(*J,&assembled);
 72:   if (assembled) {
 73:     MatZeroEntries(*J);
 74:   }

 76:   VecGetSize(xx1,&N);
 77:   VecGetOwnershipRange(xx1,&start,&end);
 78:   TSComputeRHSFunction(ts,t,xx1,jj1);

 80:   /* Compute Jacobian approximation, 1 column at a time.
 81:       xx1 = current iterate, jj1 = F(xx1)
 82:       xx2 = perturbed iterate, jj2 = F(xx2)
 83:    */
 84:   VecGetArray(xx1,&xx);
 85:   for (i=0; i<N; i++) {
 86:     VecCopy(xx1,xx2);
 87:     if (i>= start && i<end) {
 88:       dx = xx[i-start];
 89: #if !defined(PETSC_USE_COMPLEX)
 90:       if (dx < dx_min && dx >= 0.0) dx = dx_par;
 91:       else if (dx < 0.0 && dx > -dx_min) dx = -dx_par;
 92: #else
 93:       if (PetscAbsScalar(dx) < dx_min && PetscRealPart(dx) >= 0.0) dx = dx_par;
 94:       else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < dx_min) dx = -dx_par;
 95: #endif
 96:       dx *= epsilon;
 97:       wscale = 1.0/dx;
 98:       VecSetValues(xx2,1,&i,&dx,ADD_VALUES);
 99:     } else {
100:       wscale = 0.0;
101:     }
102:     TSComputeRHSFunction(ts,t,xx2,jj2);
103:     VecAXPY(jj2,-1.0,jj1);
104:     /* Communicate scale to all processors */
105:     MPI_Allreduce(&wscale,&scale,1,MPIU_SCALAR,MPIU_SUM,comm);
106:     VecScale(jj2,scale);
107:     VecGetArray(jj2,&y);
108:     VecNorm(jj2,NORM_INFINITY,&amax);
109:     amax *= 1.e-14;
110:     for (j=start; j<end; j++) {
111:       if (PetscAbsScalar(y[j-start]) > amax) {
112:         MatSetValues(*J,1,&j,1,&i,y+j-start,INSERT_VALUES);
113:       }
114:     }
115:     VecRestoreArray(jj2,&y);
116:   }

118:   VecRestoreArray(xx1,&xx);
119:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
120:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
121:   *flag =  DIFFERENT_NONZERO_PATTERN;

123:   VecDestroy(&jj1);
124:   VecDestroy(&jj2);
125:   VecDestroy(&xx2);

127:   return 0;
128: }