1 2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 3 4 /*@C 5 SNESComputeJacobianDefault - Computes the Jacobian using finite differences. 6 7 Collective on SNES 8 9 Input Parameters: 10 + x1 - compute Jacobian at this point 11 - ctx - application's function context, as set with SNESSetFunction() 12 13 Output Parameters: 14 + J - Jacobian matrix (not altered in this routine) 15 - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J) 16 17 Options Database Key: 18 + -snes_fd - Activates SNESComputeJacobianDefault() 19 . -snes_test_err - Square root of function error tolerance, default square root of machine 20 epsilon (1.e-8 in double, 3.e-4 in single) 21 - -mat_fd_type - Either wp or ds (see MATMFFD_WP or MATMFFD_DS) 22 23 Notes: 24 This routine is slow and expensive, and is not currently optimized 25 to take advantage of sparsity in the problem. Although 26 SNESComputeJacobianDefault() is not recommended for general use 27 in large-scale applications, It can be useful in checking the 28 correctness of a user-provided Jacobian. 29 30 An alternative routine that uses coloring to exploit matrix sparsity is 31 SNESComputeJacobianDefaultColor(). 32 33 Level: intermediate 34 35 .keywords: SNES, finite differences, Jacobian 36 37 .seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF() 38 @*/ 39 PetscErrorCode SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx) 40 { 41 Vec j1a,j2a,x2; 42 PetscErrorCode ierr; 43 PetscInt i,N,start,end,j,value,root; 44 PetscScalar dx,*y,wscale; 45 const PetscScalar *xx; 46 PetscReal amax,epsilon = PETSC_SQRT_MACHINE_EPSILON; 47 PetscReal dx_min = 1.e-16,dx_par = 1.e-1,unorm; 48 MPI_Comm comm; 49 PetscBool assembled,use_wp = PETSC_TRUE,flg; 50 const char *list[2] = {"ds","wp"}; 51 PetscMPIInt size; 52 const PetscInt *ranges; 53 54 PetscFunctionBegin; 55 /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */ 56 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRQ(ierr); 57 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr); 58 59 ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr); 60 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 61 ierr = MatAssembled(B,&assembled);CHKERRQ(ierr); 62 if (assembled) { 63 ierr = MatZeroEntries(B);CHKERRQ(ierr); 64 } 65 if (!snes->nvwork) { 66 snes->nvwork = 3; 67 68 ierr = VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);CHKERRQ(ierr); 69 ierr = PetscLogObjectParents(snes,snes->nvwork,snes->vwork);CHKERRQ(ierr); 70 } 71 j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; 72 73 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 74 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 75 ierr = SNESComputeFunction(snes,x1,j1a);CHKERRQ(ierr); 76 77 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES");CHKERRQ(ierr); 78 ierr = PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg);CHKERRQ(ierr); 79 ierr = PetscOptionsEnd();CHKERRQ(ierr); 80 if (flg && !value) 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 = VecGetArrayRead(x1,&xx);CHKERRQ(ierr); 93 if (use_wp) dx = PetscSqrtReal(1.0 + unorm); 94 else dx = xx[i-start]; 95 ierr = VecRestoreArrayRead(x1,&xx);CHKERRQ(ierr); 96 if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par; 97 dx *= epsilon; 98 wscale = 1.0/dx; 99 ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr); 100 } else { 101 wscale = 0.0; 102 } 103 ierr = VecAssemblyBegin(x2);CHKERRQ(ierr); 104 ierr = VecAssemblyEnd(x2);CHKERRQ(ierr); 105 ierr = SNESComputeFunction(snes,x2,j2a);CHKERRQ(ierr); 106 ierr = VecAXPY(j2a,-1.0,j1a);CHKERRQ(ierr); 107 /* Communicate scale=1/dx_i to all processors */ 108 ierr = VecGetOwnershipRanges(x1,&ranges);CHKERRQ(ierr); 109 root = size; 110 for (j=size-1; j>-1; j--) { 111 root--; 112 if (i>=ranges[j]) break; 113 } 114 ierr = MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm);CHKERRQ(ierr); 115 116 ierr = VecScale(j2a,wscale);CHKERRQ(ierr); 117 ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr); amax *= 1.e-14; 118 ierr = VecGetArray(j2a,&y);CHKERRQ(ierr); 119 for (j=start; j<end; j++) { 120 if (PetscAbsScalar(y[j-start]) > amax || j == i) { 121 ierr = MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr); 122 } 123 } 124 ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr); 125 } 126 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 127 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 128 if (B != J) { 129 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 130 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 131 } 132 PetscFunctionReturn(0); 133 } 134 135 136