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