1 2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscdm.h> 4 5 /*@C 6 SNESComputeJacobianDefault - Computes the Jacobian using finite differences. 7 8 Collective on SNES 9 10 Input Parameters: 11 + x1 - compute Jacobian at this point 12 - ctx - application's function context, as set with SNESSetFunction() 13 14 Output Parameters: 15 + J - Jacobian matrix (not altered in this routine) 16 - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J) 17 18 Options Database Key: 19 + -snes_fd - Activates SNESComputeJacobianDefault() 20 . -snes_test_err - Square root of function error tolerance, default square root of machine 21 epsilon (1.e-8 in double, 3.e-4 in single) 22 - -mat_fd_type - Either wp or ds (see MATMFFD_WP or MATMFFD_DS) 23 24 Notes: 25 This routine is slow and expensive, and is not currently optimized 26 to take advantage of sparsity in the problem. Although 27 SNESComputeJacobianDefault() is not recommended for general use 28 in large-scale applications, It can be useful in checking the 29 correctness of a user-provided Jacobian. 30 31 An alternative routine that uses coloring to exploit matrix sparsity is 32 SNESComputeJacobianDefaultColor(). 33 34 This routine ignores the maximum number of function evaluations set with SNESSetTolerances() and the function 35 evaluations it performs are not counted in what is returned by of SNESGetNumberFunctionEvals(). 36 37 Level: intermediate 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,max_funcs = snes->max_funcs; 46 PetscScalar dx,*y,wscale; 47 const PetscScalar *xx; 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 PetscBool assembled,use_wp = PETSC_TRUE,flg; 52 const char *list[2] = {"ds","wp"}; 53 PetscMPIInt size; 54 const PetscInt *ranges; 55 DM dm; 56 DMSNES dms; 57 58 PetscFunctionBegin; 59 snes->max_funcs = PETSC_MAX_INT; 60 /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */ 61 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRQ(ierr); 62 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,NULL);CHKERRQ(ierr); 63 64 ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr); 65 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 66 ierr = MatAssembled(B,&assembled);CHKERRQ(ierr); 67 if (assembled) { 68 ierr = MatZeroEntries(B);CHKERRQ(ierr); 69 } 70 if (!snes->nvwork) { 71 if (snes->dm) { 72 ierr = DMGetGlobalVector(snes->dm,&j1a);CHKERRQ(ierr); 73 ierr = DMGetGlobalVector(snes->dm,&j2a);CHKERRQ(ierr); 74 ierr = DMGetGlobalVector(snes->dm,&x2);CHKERRQ(ierr); 75 } else { 76 snes->nvwork = 3; 77 ierr = VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);CHKERRQ(ierr); 78 ierr = PetscLogObjectParents(snes,snes->nvwork,snes->vwork);CHKERRQ(ierr); 79 j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; 80 } 81 } 82 83 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 84 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 85 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 86 ierr = DMGetDMSNES(dm,&dms);CHKERRQ(ierr); 87 if (dms->ops->computemffunction) { 88 ierr = SNESComputeMFFunction(snes,x1,j1a);CHKERRQ(ierr); 89 } else { 90 ierr = SNESComputeFunction(snes,x1,j1a);CHKERRQ(ierr); 91 } 92 93 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES");CHKERRQ(ierr); 94 ierr = PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg);CHKERRQ(ierr); 95 ierr = PetscOptionsEnd();CHKERRQ(ierr); 96 if (flg && !value) use_wp = PETSC_FALSE; 97 98 if (use_wp) { 99 ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr); 100 } 101 /* Compute Jacobian approximation, 1 column at a time. 102 x1 = current iterate, j1a = F(x1) 103 x2 = perturbed iterate, j2a = F(x2) 104 */ 105 for (i=0; i<N; i++) { 106 ierr = VecCopy(x1,x2);CHKERRQ(ierr); 107 if (i>= start && i<end) { 108 ierr = VecGetArrayRead(x1,&xx);CHKERRQ(ierr); 109 if (use_wp) dx = PetscSqrtReal(1.0 + unorm); 110 else dx = xx[i-start]; 111 ierr = VecRestoreArrayRead(x1,&xx);CHKERRQ(ierr); 112 if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par; 113 dx *= epsilon; 114 wscale = 1.0/dx; 115 ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr); 116 } else { 117 wscale = 0.0; 118 } 119 ierr = VecAssemblyBegin(x2);CHKERRQ(ierr); 120 ierr = VecAssemblyEnd(x2);CHKERRQ(ierr); 121 if (dms->ops->computemffunction) { 122 ierr = SNESComputeMFFunction(snes,x2,j2a);CHKERRQ(ierr); 123 } else { 124 ierr = SNESComputeFunction(snes,x2,j2a);CHKERRQ(ierr); 125 } 126 ierr = VecAXPY(j2a,-1.0,j1a);CHKERRQ(ierr); 127 /* Communicate scale=1/dx_i to all processors */ 128 ierr = VecGetOwnershipRanges(x1,&ranges);CHKERRQ(ierr); 129 root = size; 130 for (j=size-1; j>-1; j--) { 131 root--; 132 if (i>=ranges[j]) break; 133 } 134 ierr = MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm);CHKERRMPI(ierr); 135 ierr = VecScale(j2a,wscale);CHKERRQ(ierr); 136 ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr); amax *= 1.e-14; 137 ierr = VecGetArray(j2a,&y);CHKERRQ(ierr); 138 for (j=start; j<end; j++) { 139 if (PetscAbsScalar(y[j-start]) > amax || j == i) { 140 ierr = MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr); 141 } 142 } 143 ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr); 144 } 145 if (snes->dm) { 146 ierr = DMRestoreGlobalVector(snes->dm,&j1a);CHKERRQ(ierr); 147 ierr = DMRestoreGlobalVector(snes->dm,&j2a);CHKERRQ(ierr); 148 ierr = DMRestoreGlobalVector(snes->dm,&x2);CHKERRQ(ierr); 149 } 150 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 151 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 152 if (B != J) { 153 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 154 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 155 } 156 snes->max_funcs = max_funcs; 157 snes->nfuncs -= N; 158 PetscFunctionReturn(0); 159 } 160 161