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