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