1 #define PETSCSNES_DLL 2 3 #include "private/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 - -mat_fd_type - Either wp or ds (see MATMFFD_WP or MATMFFD_DS) 26 27 Notes: 28 This routine is slow and expensive, and is not currently optimized 29 to take advantage of sparsity in the problem. Although 30 SNESDefaultComputeJacobian() is not recommended for general use 31 in large-scale applications, It can be useful in checking the 32 correctness of a user-provided Jacobian. 33 34 An alternative routine that uses coloring to exploit matrix sparsity is 35 SNESDefaultComputeJacobianColor(). 36 37 Level: intermediate 38 39 .keywords: SNES, finite differences, Jacobian 40 41 .seealso: SNESSetJacobian(), SNESDefaultComputeJacobianColor(), MatCreateSNESMF() 42 @*/ 43 PetscErrorCode PETSCSNES_DLLEXPORT SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 44 { 45 Vec j1a,j2a,x2; 46 PetscErrorCode ierr; 47 PetscInt i,N,start,end,j,value,root; 48 PetscScalar dx,*y,*xx,wscale; 49 PetscReal amax,epsilon = PETSC_SQRT_MACHINE_EPSILON; 50 PetscReal dx_min = 1.e-16,dx_par = 1.e-1,unorm; 51 MPI_Comm comm; 52 PetscErrorCode (*eval_fct)(SNES,Vec,Vec)=0; 53 PetscTruth assembled,use_wp = PETSC_TRUE,flg; 54 const char *list[2] = {"ds","wp"}; 55 PetscMPIInt size; 56 const PetscInt *ranges; 57 58 PetscFunctionBegin; 59 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr); 60 eval_fct = SNESComputeFunction; 61 62 ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr); 63 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 64 ierr = MatAssembled(*B,&assembled);CHKERRQ(ierr); 65 if (assembled) { 66 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 67 } 68 if (!snes->nvwork) { 69 snes->nvwork = 3; 70 ierr = VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);CHKERRQ(ierr); 71 ierr = PetscLogObjectParents(snes,snes->nvwork,snes->vwork);CHKERRQ(ierr); 72 } 73 j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; 74 75 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 76 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 77 ierr = (*eval_fct)(snes,x1,j1a);CHKERRQ(ierr); 78 79 ierr = PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESDefaultComputeJacobian",list,2,"wp",&value,&flg);CHKERRQ(ierr); 80 if (flg && !value) { 81 use_wp = PETSC_FALSE; 82 } 83 if (use_wp) { 84 ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr); 85 } 86 /* Compute Jacobian approximation, 1 column at a time. 87 x1 = current iterate, j1a = F(x1) 88 x2 = perturbed iterate, j2a = F(x2) 89 */ 90 for (i=0; i<N; i++) { 91 ierr = VecCopy(x1,x2);CHKERRQ(ierr); 92 if (i>= start && i<end) { 93 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr); 94 if (use_wp) { 95 dx = 1.0 + unorm; 96 } else { 97 dx = xx[i-start]; 98 } 99 ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 100 #if !defined(PETSC_USE_COMPLEX) 101 if (dx < dx_min && dx >= 0.0) dx = dx_par; 102 else if (dx < 0.0 && dx > -dx_min) dx = -dx_par; 103 #else 104 if (PetscAbsScalar(dx) < dx_min && PetscRealPart(dx) >= 0.0) dx = dx_par; 105 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < dx_min) dx = -dx_par; 106 #endif 107 dx *= epsilon; 108 wscale = 1.0/dx; 109 ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr); 110 } else { 111 wscale = 0.0; 112 } 113 ierr = VecAssemblyBegin(x2);CHKERRQ(ierr); 114 ierr = VecAssemblyEnd(x2);CHKERRQ(ierr); 115 ierr = (*eval_fct)(snes,x2,j2a);CHKERRQ(ierr); 116 ierr = VecAXPY(j2a,-1.0,j1a);CHKERRQ(ierr); 117 /* Communicate scale=1/dx_i to all processors */ 118 ierr = VecGetOwnershipRanges(x1,&ranges);CHKERRQ(ierr); 119 root = size; 120 for (j=size-1; j>-1; j--){ 121 root--; 122 if (i>=ranges[j]) break; 123 } 124 ierr = MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm);CHKERRQ(ierr); 125 126 ierr = VecScale(j2a,wscale);CHKERRQ(ierr); 127 ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr); amax *= 1.e-14; 128 ierr = VecGetArray(j2a,&y);CHKERRQ(ierr); 129 for (j=start; j<end; j++) { 130 if (PetscAbsScalar(y[j-start]) > amax || j == i) { 131 ierr = MatSetValues(*B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr); 132 } 133 } 134 ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr); 135 } 136 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 137 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 138 if (*B != *J) { 139 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 140 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 141 } 142 *flag = DIFFERENT_NONZERO_PATTERN; 143 PetscFunctionReturn(0); 144 } 145 146 147