1 /*$Id: snesj.c,v 1.64 2000/04/09 04:38:32 bsmith Exp bsmith $*/ 2 3 #include "src/snes/snesimpl.h" /*I "snes.h" I*/ 4 5 #undef __FUNC__ 6 #define __FUNC__ /*<a name=""></a>*/"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 PetscReal amax,epsilon = 1.e-8; /* assumes PetscReal precision */ 46 PetscReal 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 && PetscRealPart(dx) >= 0.0) dx = dx_par; 83 else if (PetscRealPart(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__ /*<a name=""></a>*/"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,MatStructure *flag,void *ctx) 145 { 146 int ierr; 147 148 PetscFunctionBegin; 149 ierr = SNESDefaultComputeJacobian(snes,x1,J,B,flag,ctx);CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNC__ 154 #define __FUNC__ /*<a name=""></a>*/"SNESDefaultComputeHessianColor" 155 /*@C 156 SNESDefaultComputeHessianColor - Computes the Hessian using colored finite differences. 157 158 Collective on SNES 159 160 Input Parameters: 161 + x1 - compute Hessian at this point 162 - ctx - application's gradient context, as set with SNESSetGradient() 163 164 Output Parameters: 165 + J - Hessian matrix (not altered in this routine) 166 . B - newly computed Hessian matrix to use with preconditioner (generally the same as J) 167 - flag - flag indicating whether the matrix sparsity structure has changed 168 169 Options Database Keys: 170 . -mat_fd_coloring_freq <freq> - Activates SNESDefaultComputeJacobianColor() 171 172 Level: intermediate 173 174 .keywords: SNES, finite differences, Hessian, coloring, sparse 175 176 .seealso: SNESSetHessian() 177 @*/ 178 int SNESDefaultComputeHessianColor(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 179 { 180 int ierr; 181 182 PetscFunctionBegin; 183 ierr = SNESDefaultComputeJacobianColor(snes,x1,J,B,flag,ctx);CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 187