1 #define PETSCSNES_DLL 2 3 #include "src/snes/snesimpl.h" /*I "petscsnes.h" I*/ 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "SNESDefaultComputeJacobian" 7 /*@C 8 SNESDefaultComputeJacobian - Computes the Jacobian using finite differences. 9 10 Collective on SNES 11 12 Input Parameters: 13 + x1 - compute Jacobian at this point 14 - ctx - application's function context, as set with SNESSetFunction() 15 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 20 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) 25 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. 32 33 An alternative routine that uses coloring to explot matrix sparsity is 34 SNESDefaultComputeJacobianColor(). 35 36 Level: intermediate 37 38 .keywords: SNES, finite differences, Jacobian 39 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; 45 PetscErrorCode ierr; 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; 53 54 PetscFunctionBegin; 55 ierr = PetscOptionsGetReal(snes->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr); 56 eval_fct = SNESComputeFunction; 57 58 ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr); 59 ierr = MatAssembled(*B,&assembled);CHKERRQ(ierr); 60 if (assembled) { 61 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 62 } 63 if (!snes->nvwork) { 64 ierr = VecDuplicateVecs(x1,3,&snes->vwork);CHKERRQ(ierr); 65 snes->nvwork = 3; 66 ierr = PetscLogObjectParents(snes,3,snes->vwork);CHKERRQ(ierr); 67 } 68 j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; 69 70 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 71 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 72 ierr = (*eval_fct)(snes,x1,j1a);CHKERRQ(ierr); 73 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 ierr = VecCopy(x1,x2);CHKERRQ(ierr); 80 if (i>= start && i<end) { 81 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr); 82 dx = xx[i-start]; 83 ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 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 ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr); 94 } else { 95 wscale = 0.0; 96 } 97 ierr = (*eval_fct)(snes,x2,j2a);CHKERRQ(ierr); 98 ierr = VecAXPY(j2a,mone,j1a);CHKERRQ(ierr); 99 /* Communicate scale to all processors */ 100 ierr = MPI_Allreduce(&wscale,&scale,1,MPIU_SCALAR,PetscSum_Op,comm);CHKERRQ(ierr); 101 ierr = VecScale(j2a,scale);CHKERRQ(ierr); 102 ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr); amax *= 1.e-14; 103 ierr = VecGetArray(j2a,&y);CHKERRQ(ierr); 104 for (j=start; j<end; j++) { 105 if (PetscAbsScalar(y[j-start]) > amax) { 106 ierr = MatSetValues(*B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr); 107 } 108 } 109 ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr); 110 } 111 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 112 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 113 if (*B != *J) { 114 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 115 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 116 } 117 *flag = DIFFERENT_NONZERO_PATTERN; 118 PetscFunctionReturn(0); 119 } 120 121 122