1 2 #ifndef lint 3 static char vcid[] = "$Id: snesj.c,v 1.41 1997/01/06 20:29:45 balay Exp bsmith $"; 4 #endif 5 6 #include "src/snes/snesimpl.h" /*I "snes.h" I*/ 7 8 #undef __FUNC__ 9 #define __FUNC__ "SNESDefaultComputeJacobian" 10 /*@C 11 SNESDefaultComputeJacobian - Computes the Jacobian using finite differences. 12 13 Input Parameters: 14 . x1 - compute Jacobian at this point 15 . ctx - application's function context, as set with SNESSetFunction() 16 17 Output Parameters: 18 . J - Jacobian 19 . B - preconditioner, same as Jacobian 20 . flag - matrix flag 21 22 Options Database Key: 23 $ -snes_fd 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 .keywords: SNES, finite differences, Jacobian 33 34 .seealso: SNESSetJacobian() 35 @*/ 36 int SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 37 { 38 Vec j1,j2,x2; 39 int i,ierr,N,start,end,j; 40 Scalar dx, mone = -1.0,*y,scale,*xx,wscale; 41 double amax, epsilon = 1.e-8; /* assumes double precision */ 42 double dx_min = 1.e-16, dx_par = 1.e-1; 43 MPI_Comm comm; 44 int (*eval_fct)(SNES,Vec,Vec); 45 46 if (snes->method_class == SNES_NONLINEAR_EQUATIONS) 47 eval_fct = SNESComputeFunction; 48 else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) 49 eval_fct = SNESComputeGradient; 50 else SETERRQ(1,0,"Invalid method class"); 51 52 PetscObjectGetComm((PetscObject)x1,&comm); 53 MatZeroEntries(*J); 54 if (!snes->nvwork) { 55 ierr = VecDuplicateVecs(x1,3,&snes->vwork); CHKERRQ(ierr); 56 snes->nvwork = 3; 57 PLogObjectParents(snes,3,snes->vwork); 58 } 59 j1 = snes->vwork[0]; j2 = snes->vwork[1]; x2 = snes->vwork[2]; 60 61 ierr = VecGetSize(x1,&N); CHKERRQ(ierr); 62 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERRQ(ierr); 63 VecGetArray(x1,&xx); 64 ierr = eval_fct(snes,x1,j1); CHKERRQ(ierr); 65 66 /* Compute Jacobian approximation, 1 column at a time. 67 x1 = current iterate, j1 = F(x1) 68 x2 = perturbed iterate, j2 = F(x2) 69 */ 70 for ( i=0; i<N; i++ ) { 71 ierr = VecCopy(x1,x2); CHKERRQ(ierr); 72 if ( i>= start && i<end) { 73 dx = xx[i-start]; 74 #if !defined(PETSC_COMPLEX) 75 if (dx < dx_min && dx >= 0.0) dx = dx_par; 76 else if (dx < 0.0 && dx > -dx_min) dx = -dx_par; 77 #else 78 if (abs(dx) < dx_min && real(dx) >= 0.0) dx = dx_par; 79 else if (real(dx) < 0.0 && abs(dx) < dx_min) dx = -dx_par; 80 #endif 81 dx *= epsilon; 82 wscale = 1.0/dx; 83 VecSetValues(x2,1,&i,&dx,ADD_VALUES); 84 } 85 else { 86 wscale = 0.0; 87 } 88 ierr = eval_fct(snes,x2,j2); CHKERRQ(ierr); 89 ierr = VecAXPY(&mone,j1,j2); CHKERRQ(ierr); 90 /* Communicate scale to all processors */ 91 #if !defined(PETSC_COMPLEX) 92 MPI_Allreduce(&wscale,&scale,1,MPI_DOUBLE,MPI_SUM,comm); 93 #else 94 MPI_Allreduce(&wscale,&scale,2,MPI_DOUBLE,MPI_SUM,comm); 95 #endif 96 VecScale(&scale,j2); 97 VecGetArray(j2,&y); 98 VecNorm(j2,NORM_INFINITY,&amax); amax *= 1.e-14; 99 for ( j=start; j<end; j++ ) { 100 if (PetscAbsScalar(y[j-start]) > amax) { 101 ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERT_VALUES); CHKERRQ(ierr); 102 } 103 } 104 VecRestoreArray(j2,&y); 105 } 106 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 107 ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 108 *flag = DIFFERENT_NONZERO_PATTERN; 109 return 0; 110 } 111 112 #undef __FUNC__ 113 #define __FUNC__ "SNESDefaultComputeHessian" 114 /*@C 115 SNESDefaultComputeHessian - Computes the Hessian using finite differences. 116 117 Input Parameters: 118 . x1 - compute Hessian at this point 119 . ctx - application's gradient context, as set with SNESSetGradient() 120 121 Output Parameters: 122 . J - Hessian 123 . B - preconditioner, same as Hessian 124 . flag - matrix flag 125 126 Options Database Key: 127 $ -snes_fd 128 129 Notes: 130 This routine is slow and expensive, and is not currently optimized 131 to take advantage of sparsity in the problem. Although 132 SNESDefaultComputeHessian() is not recommended for general use 133 in large-scale applications, It can be useful in checking the 134 correctness of a user-provided Hessian. 135 136 .keywords: SNES, finite differences, Hessian 137 138 .seealso: SNESSetHessian() 139 @*/ 140 int SNESDefaultComputeHessian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 141 { 142 return SNESDefaultComputeJacobian(snes,x1,J,B,flag,ctx); 143 } 144