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