1 2 #ifndef lint 3 static char vcid[] = "$Id: snesj.c,v 1.33 1996/08/08 14:46:41 bsmith Exp bsmith $"; 4 #endif 5 6 #include "draw.h" /*I "draw.h" I*/ 7 #include "src/snes/snesimpl.h" /*I "snes.h" I*/ 8 9 /*@C 10 SNESDefaultComputeJacobian - Computes the Jacobian using finite 11 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(), SNESTestJacobian() 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 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,"SNESDefaultComputeJacobian: 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 < 1.e-16 && dx >= 0.0) dx = 1.e-1; 75 else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 76 #else 77 if (abs(dx) < 1.e-16 && real(dx) >= 0.0) dx = 1.e-1; 78 else if (real(dx) < 0.0 && abs(dx) < 1.e-16) dx = -1.e-1; 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 /*@C 112 SNESDefaultComputeHessian - Computes the Hessian using finite 113 differences. 114 115 Input Parameters: 116 . x1 - compute Hessian at this point 117 . ctx - application's gradient context, as set with SNESSetGradient() 118 119 Output Parameters: 120 . J - Hessian 121 . B - preconditioner, same as Hessian 122 . flag - matrix flag 123 124 Options Database Key: 125 $ -snes_fd 126 127 Notes: 128 This routine is slow and expensive, and is not currently optimized 129 to take advantage of sparsity in the problem. Although 130 SNESDefaultComputeHessian() is not recommended for general use 131 in large-scale applications, It can be useful in checking the 132 correctness of a user-provided Hessian. 133 134 .keywords: SNES, finite differences, Hessian 135 136 .seealso: SNESSetHessian(), SNESTestHessian() 137 @*/ 138 int SNESDefaultComputeHessian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 139 { 140 return SNESDefaultComputeJacobian(snes,x1,J,B,flag,ctx); 141 } 142