1 2 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "SNESComputeJacobianDefault" 6 /*@C 7 SNESComputeJacobianDefault - 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 19 Options Database Key: 20 + -snes_fd - Activates SNESComputeJacobianDefault() 21 . -snes_test_err - Square root of function error tolerance, default square root of machine 22 epsilon (1.e-8 in double, 3.e-4 in single) 23 - -mat_fd_type - Either wp or ds (see MATMFFD_WP or MATMFFD_DS) 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 SNESComputeJacobianDefault() 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 exploit matrix sparsity is 33 SNESComputeJacobianDefaultColor(). 34 35 Level: intermediate 36 37 .keywords: SNES, finite differences, Jacobian 38 39 .seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF() 40 @*/ 41 PetscErrorCode SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx) 42 { 43 Vec j1a,j2a,x2; 44 PetscErrorCode ierr; 45 PetscInt i,N,start,end,j,value,root; 46 PetscScalar dx,*y,*xx,wscale; 47 PetscReal amax,epsilon = PETSC_SQRT_MACHINE_EPSILON; 48 PetscReal dx_min = 1.e-16,dx_par = 1.e-1,unorm; 49 MPI_Comm comm; 50 PetscErrorCode (*eval_fct)(SNES,Vec,Vec)=0; 51 PetscBool assembled,use_wp = PETSC_TRUE,flg; 52 const char *list[2] = {"ds","wp"}; 53 PetscMPIInt size; 54 const PetscInt *ranges; 55 56 PetscFunctionBegin; 57 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr); 58 eval_fct = SNESComputeFunction; 59 60 ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr); 61 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 62 ierr = MatAssembled(B,&assembled);CHKERRQ(ierr); 63 if (assembled) { 64 ierr = MatZeroEntries(B);CHKERRQ(ierr); 65 } 66 if (!snes->nvwork) { 67 snes->nvwork = 3; 68 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 = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES"); 79 ierr = PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg);CHKERRQ(ierr); 80 ierr = PetscOptionsEnd(); 81 if (flg && !value) 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) dx = 1.0 + unorm; 95 else dx = xx[i-start]; 96 ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 97 if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par; 98 dx *= epsilon; 99 wscale = 1.0/dx; 100 ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr); 101 } else { 102 wscale = 0.0; 103 } 104 ierr = VecAssemblyBegin(x2);CHKERRQ(ierr); 105 ierr = VecAssemblyEnd(x2);CHKERRQ(ierr); 106 ierr = (*eval_fct)(snes,x2,j2a);CHKERRQ(ierr); 107 ierr = VecAXPY(j2a,-1.0,j1a);CHKERRQ(ierr); 108 /* Communicate scale=1/dx_i to all processors */ 109 ierr = VecGetOwnershipRanges(x1,&ranges);CHKERRQ(ierr); 110 root = size; 111 for (j=size-1; j>-1; j--) { 112 root--; 113 if (i>=ranges[j]) break; 114 } 115 ierr = MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm);CHKERRQ(ierr); 116 117 ierr = VecScale(j2a,wscale);CHKERRQ(ierr); 118 ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr); amax *= 1.e-14; 119 ierr = VecGetArray(j2a,&y);CHKERRQ(ierr); 120 for (j=start; j<end; j++) { 121 if (PetscAbsScalar(y[j-start]) > amax || j == i) { 122 ierr = MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr); 123 } 124 } 125 ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr); 126 } 127 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 128 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 129 if (B != J) { 130 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 131 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 132 } 133 PetscFunctionReturn(0); 134 } 135 136 137