1 2 #include "src/snes/snesimpl.h" /*I "petscsnes.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "SNESDefaultComputeJacobian" 6 /*@C 7 SNESDefaultComputeJacobian - Computes the Jacobian using finite differences. 8 9 Collective on SNES 10 11 Input Parameters: 12 + x1 - compute Jacobian at this point 13 - ctx - application's function context, as set with SNESSetFunction() 14 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) 18 - flag - flag indicating whether the matrix sparsity structure has changed 19 20 Options Database Key: 21 + -snes_fd - Activates SNESDefaultComputeJacobian() 22 - -snes_test_err - Square root of function error tolerance, default square root of machine 23 epsilon (1.e-8 in double, 3.e-4 in single) 24 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 SNESDefaultComputeJacobian() 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. 31 32 An alternative routine that uses coloring to explot matrix sparsity is 33 SNESDefaultComputeJacobianColor(). 34 35 Level: intermediate 36 37 .keywords: SNES, finite differences, Jacobian 38 39 .seealso: SNESSetJacobian(), SNESDefaultComputeJacobianColor() 40 @*/ 41 PetscErrorCode SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 42 { 43 Vec j1a,j2a,x2; 44 PetscErrorCode ierr; 45 int i,N,start,end,j; 46 PetscScalar dx,mone = -1.0,*y,scale,*xx,wscale; 47 PetscReal amax,epsilon = PETSC_SQRT_MACHINE_EPSILON; 48 PetscReal dx_min = 1.e-16,dx_par = 1.e-1; 49 MPI_Comm comm; 50 PetscErrorCode (*eval_fct)(SNES,Vec,Vec)=0; 51 52 PetscFunctionBegin; 53 ierr = PetscOptionsGetReal(snes->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr); 54 eval_fct = SNESComputeFunction; 55 56 ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr); 57 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 58 if (!snes->nvwork) { 59 ierr = VecDuplicateVecs(x1,3,&snes->vwork);CHKERRQ(ierr); 60 snes->nvwork = 3; 61 PetscLogObjectParents(snes,3,snes->vwork); 62 } 63 j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; 64 65 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 66 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 67 ierr = (*eval_fct)(snes,x1,j1a);CHKERRQ(ierr); 68 69 /* Compute Jacobian approximation, 1 column at a time. 70 x1 = current iterate, j1a = F(x1) 71 x2 = perturbed iterate, j2a = F(x2) 72 */ 73 for (i=0; i<N; i++) { 74 ierr = VecCopy(x1,x2);CHKERRQ(ierr); 75 if (i>= start && i<end) { 76 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr); 77 dx = xx[i-start]; 78 ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 79 #if !defined(PETSC_USE_COMPLEX) 80 if (dx < dx_min && dx >= 0.0) dx = dx_par; 81 else if (dx < 0.0 && dx > -dx_min) dx = -dx_par; 82 #else 83 if (PetscAbsScalar(dx) < dx_min && PetscRealPart(dx) >= 0.0) dx = dx_par; 84 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < dx_min) dx = -dx_par; 85 #endif 86 dx *= epsilon; 87 wscale = 1.0/dx; 88 ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr); 89 } else { 90 wscale = 0.0; 91 } 92 ierr = (*eval_fct)(snes,x2,j2a);CHKERRQ(ierr); 93 ierr = VecAXPY(&mone,j1a,j2a);CHKERRQ(ierr); 94 /* Communicate scale to all processors */ 95 ierr = MPI_Allreduce(&wscale,&scale,1,MPIU_SCALAR,PetscSum_Op,comm);CHKERRQ(ierr); 96 ierr = VecScale(&scale,j2a);CHKERRQ(ierr); 97 ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr); amax *= 1.e-14; 98 ierr = VecGetArray(j2a,&y);CHKERRQ(ierr); 99 for (j=start; j<end; j++) { 100 if (PetscAbsScalar(y[j-start]) > amax) { 101 ierr = MatSetValues(*B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr); 102 } 103 } 104 ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr); 105 } 106 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 107 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 108 *flag = DIFFERENT_NONZERO_PATTERN; 109 PetscFunctionReturn(0); 110 } 111 112 113