1 /*$Id: snesj.c,v 1.75 2001/09/11 18:06:40 bsmith Exp $*/ 2 3 #include "src/snes/snesimpl.h" /*I "petscsnes.h" I*/ 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "SNESDefaultComputeJacobian" 7 /*@C 8 SNESDefaultComputeJacobian - Computes the Jacobian using finite differences. 9 10 Collective on SNES 11 12 Input Parameters: 13 + x1 - compute Jacobian at this point 14 - ctx - application's function context, as set with SNESSetFunction() 15 16 Output Parameters: 17 + J - Jacobian matrix (not altered in this routine) 18 . B - newly computed Jacobian matrix to use with preconditioner (generally the same as J) 19 - flag - flag indicating whether the matrix sparsity structure has changed 20 21 Options Database Key: 22 + -snes_fd - Activates SNESDefaultComputeJacobian() 23 - -snes_test_err - Square root of function error tolerance, default square root of machine 24 epsilon (1.e-8 in double, 3.e-4 in single) 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 explot matrix sparsity is 34 SNESDefaultComputeJacobianColor(). 35 36 Level: intermediate 37 38 .keywords: SNES, finite differences, Jacobian 39 40 .seealso: SNESSetJacobian(), SNESDefaultComputeJacobianColor() 41 @*/ 42 int SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 43 { 44 Vec j1a,j2a,x2; 45 int i,ierr,N,start,end,j; 46 PetscScalar dx,mone = -1.0,*y,scale,*xx,wscale; 47 PetscReal amax,epsilon = PETSC_SQRT_MACHINE_EPSILON; 48 PetscReal dx_min = 1.e-16,dx_par = 1.e-1; 49 MPI_Comm comm; 50 int (*eval_fct)(SNES,Vec,Vec)=0; 51 52 PetscFunctionBegin; 53 ierr = PetscOptionsGetReal(snes->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr); 54 if (snes->method_class == SNES_NONLINEAR_EQUATIONS) eval_fct = SNESComputeFunction; 55 else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) eval_fct = SNESComputeGradient; 56 else SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid method class"); 57 58 ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr); 59 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 60 if (!snes->nvwork) { 61 ierr = VecDuplicateVecs(x1,3,&snes->vwork);CHKERRQ(ierr); 62 snes->nvwork = 3; 63 PetscLogObjectParents(snes,3,snes->vwork); 64 } 65 j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; 66 67 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 68 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 69 ierr = (*eval_fct)(snes,x1,j1a);CHKERRQ(ierr); 70 71 /* Compute Jacobian approximation, 1 column at a time. 72 x1 = current iterate, j1a = F(x1) 73 x2 = perturbed iterate, j2a = F(x2) 74 */ 75 for (i=0; i<N; i++) { 76 ierr = VecCopy(x1,x2);CHKERRQ(ierr); 77 if (i>= start && i<end) { 78 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr); 79 dx = xx[i-start]; 80 ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 81 #if !defined(PETSC_USE_COMPLEX) 82 if (dx < dx_min && dx >= 0.0) dx = dx_par; 83 else if (dx < 0.0 && dx > -dx_min) dx = -dx_par; 84 #else 85 if (PetscAbsScalar(dx) < dx_min && PetscRealPart(dx) >= 0.0) dx = dx_par; 86 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < dx_min) dx = -dx_par; 87 #endif 88 dx *= epsilon; 89 wscale = 1.0/dx; 90 ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr); 91 } else { 92 wscale = 0.0; 93 } 94 ierr = (*eval_fct)(snes,x2,j2a);CHKERRQ(ierr); 95 ierr = VecAXPY(&mone,j1a,j2a);CHKERRQ(ierr); 96 /* Communicate scale to all processors */ 97 ierr = MPI_Allreduce(&wscale,&scale,1,MPIU_SCALAR,PetscSum_Op,comm);CHKERRQ(ierr); 98 ierr = VecScale(&scale,j2a);CHKERRQ(ierr); 99 ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr); amax *= 1.e-14; 100 ierr = VecGetArray(j2a,&y);CHKERRQ(ierr); 101 for (j=start; j<end; j++) { 102 if (PetscAbsScalar(y[j-start]) > amax) { 103 ierr = MatSetValues(*B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr); 104 } 105 } 106 ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr); 107 } 108 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 109 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 110 *flag = DIFFERENT_NONZERO_PATTERN; 111 PetscFunctionReturn(0); 112 } 113 114 #undef __FUNCT__ 115 #define __FUNCT__ "SNESDefaultComputeHessian" 116 /*@C 117 SNESDefaultComputeHessian - Computes the Hessian using finite differences. 118 119 Collective on SNES 120 121 Input Parameters: 122 + x1 - compute Hessian at this point 123 - ctx - application's gradient context, as set with SNESSetGradient() 124 125 Output Parameters: 126 + J - Hessian matrix (not altered in this routine) 127 . B - newly computed Hessian matrix to use with preconditioner (generally the same as J) 128 - flag - flag indicating whether the matrix sparsity structure has changed 129 130 Options Database Key: 131 $ -snes_fd - Activates SNESDefaultComputeHessian() 132 133 134 Level: intermediate 135 136 Notes: 137 This routine is slow and expensive, and is not currently optimized 138 to take advantage of sparsity in the problem. Although 139 SNESDefaultComputeHessian() is not recommended for general use 140 in large-scale applications, It can be useful in checking the 141 correctness of a user-provided Hessian. 142 143 .keywords: SNES, finite differences, Hessian 144 145 .seealso: SNESSetHessian() 146 @*/ 147 int SNESDefaultComputeHessian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 148 { 149 int ierr; 150 151 PetscFunctionBegin; 152 ierr = SNESDefaultComputeJacobian(snes,x1,J,B,flag,ctx);CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 #undef __FUNCT__ 157 #define __FUNCT__ "SNESDefaultComputeHessianColor" 158 /*@C 159 SNESDefaultComputeHessianColor - Computes the Hessian using colored finite differences. 160 161 Collective on SNES 162 163 Input Parameters: 164 + x1 - compute Hessian at this point 165 - ctx - application's gradient context, as set with SNESSetGradient() 166 167 Output Parameters: 168 + J - Hessian matrix (not altered in this routine) 169 . B - newly computed Hessian matrix to use with preconditioner (generally the same as J) 170 - flag - flag indicating whether the matrix sparsity structure has changed 171 172 Options Database Keys: 173 . -mat_fd_coloring_freq <freq> - Activates SNESDefaultComputeJacobianColor() 174 175 Level: intermediate 176 177 .keywords: SNES, finite differences, Hessian, coloring, sparse 178 179 .seealso: SNESSetHessian() 180 @*/ 181 int SNESDefaultComputeHessianColor(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 182 { 183 int ierr; 184 185 PetscFunctionBegin; 186 ierr = SNESDefaultComputeJacobianColor(snes,x1,J,B,flag,ctx);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190