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 + snes - the SNES context 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 This routine ignores the maximum number of function evaluations set with SNESSetTolerances() and the function 36 evaluations it performs are not counted in what is returned by of SNESGetNumberFunctionEvals(). 37 38 Level: intermediate 39 40 .seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF() 41 @*/ 42 PetscErrorCode SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx) 43 { 44 Vec j1a,j2a,x2; 45 PetscErrorCode ierr; 46 PetscInt i,N,start,end,j,value,root,max_funcs = snes->max_funcs; 47 PetscScalar dx,*y,wscale; 48 const PetscScalar *xx; 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 PetscBool assembled,use_wp = PETSC_TRUE,flg; 53 const char *list[2] = {"ds","wp"}; 54 PetscMPIInt size; 55 const PetscInt *ranges; 56 DM dm; 57 DMSNES dms; 58 59 PetscFunctionBegin; 60 snes->max_funcs = PETSC_MAX_INT; 61 /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */ 62 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRQ(ierr); 63 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,NULL);CHKERRQ(ierr); 64 65 ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr); 66 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 67 ierr = MatAssembled(B,&assembled);CHKERRQ(ierr); 68 if (assembled) { 69 ierr = MatZeroEntries(B);CHKERRQ(ierr); 70 } 71 if (!snes->nvwork) { 72 if (snes->dm) { 73 ierr = DMGetGlobalVector(snes->dm,&j1a);CHKERRQ(ierr); 74 ierr = DMGetGlobalVector(snes->dm,&j2a);CHKERRQ(ierr); 75 ierr = DMGetGlobalVector(snes->dm,&x2);CHKERRQ(ierr); 76 } else { 77 snes->nvwork = 3; 78 ierr = VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);CHKERRQ(ierr); 79 ierr = PetscLogObjectParents(snes,snes->nvwork,snes->vwork);CHKERRQ(ierr); 80 j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; 81 } 82 } 83 84 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 85 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 86 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 87 ierr = DMGetDMSNES(dm,&dms);CHKERRQ(ierr); 88 if (dms->ops->computemffunction) { 89 ierr = SNESComputeMFFunction(snes,x1,j1a);CHKERRQ(ierr); 90 } else { 91 ierr = SNESComputeFunction(snes,x1,j1a);CHKERRQ(ierr); 92 } 93 94 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES");CHKERRQ(ierr); 95 ierr = PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg);CHKERRQ(ierr); 96 ierr = PetscOptionsEnd();CHKERRQ(ierr); 97 if (flg && !value) use_wp = PETSC_FALSE; 98 99 if (use_wp) { 100 ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr); 101 } 102 /* Compute Jacobian approximation, 1 column at a time. 103 x1 = current iterate, j1a = F(x1) 104 x2 = perturbed iterate, j2a = F(x2) 105 */ 106 for (i=0; i<N; i++) { 107 ierr = VecCopy(x1,x2);CHKERRQ(ierr); 108 if (i>= start && i<end) { 109 ierr = VecGetArrayRead(x1,&xx);CHKERRQ(ierr); 110 if (use_wp) dx = PetscSqrtReal(1.0 + unorm); 111 else dx = xx[i-start]; 112 ierr = VecRestoreArrayRead(x1,&xx);CHKERRQ(ierr); 113 if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par; 114 dx *= epsilon; 115 wscale = 1.0/dx; 116 ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr); 117 } else { 118 wscale = 0.0; 119 } 120 ierr = VecAssemblyBegin(x2);CHKERRQ(ierr); 121 ierr = VecAssemblyEnd(x2);CHKERRQ(ierr); 122 if (dms->ops->computemffunction) { 123 ierr = SNESComputeMFFunction(snes,x2,j2a);CHKERRQ(ierr); 124 } else { 125 ierr = SNESComputeFunction(snes,x2,j2a);CHKERRQ(ierr); 126 } 127 ierr = VecAXPY(j2a,-1.0,j1a);CHKERRQ(ierr); 128 /* Communicate scale=1/dx_i to all processors */ 129 ierr = VecGetOwnershipRanges(x1,&ranges);CHKERRQ(ierr); 130 root = size; 131 for (j=size-1; j>-1; j--) { 132 root--; 133 if (i>=ranges[j]) break; 134 } 135 ierr = MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm);CHKERRMPI(ierr); 136 ierr = VecScale(j2a,wscale);CHKERRQ(ierr); 137 ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr); amax *= 1.e-14; 138 ierr = VecGetArray(j2a,&y);CHKERRQ(ierr); 139 for (j=start; j<end; j++) { 140 if (PetscAbsScalar(y[j-start]) > amax || j == i) { 141 ierr = MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr); 142 } 143 } 144 ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr); 145 } 146 if (snes->dm) { 147 ierr = DMRestoreGlobalVector(snes->dm,&j1a);CHKERRQ(ierr); 148 ierr = DMRestoreGlobalVector(snes->dm,&j2a);CHKERRQ(ierr); 149 ierr = DMRestoreGlobalVector(snes->dm,&x2);CHKERRQ(ierr); 150 } 151 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 152 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 153 if (B != J) { 154 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 155 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 156 } 157 snes->max_funcs = max_funcs; 158 snes->nfuncs -= N; 159 PetscFunctionReturn(0); 160 } 161 162