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