1 2 #ifndef lint 3 static char vcid[] = "$Id: snesj.c,v 1.8 1995/05/11 17:44:07 bsmith Exp curfman $"; 4 #endif 5 6 #include "draw.h" 7 #include "snes.h" 8 #include "options.h" 9 10 /*@ 11 SNESDefaultComputeJacobian - Computes the Jacobian using finite 12 differences. 13 14 Input Parameters: 15 . x - compute Jacobian at this point 16 . ctx - application's Function context 17 18 Output Parameters: 19 . J - Jacobian 20 . B - preconditioner, same as Jacobian 21 22 Notes: 23 This routine is slow and expensive, and is not currently optimized 24 to take advantage of sparsity in the problem. Although 25 SNESDefaultComputeJacobian() is not recommended for general use 26 in large-scale applications, It can be useful in checking the 27 correctness of a user-provided Jacobian. 28 29 .keywords: SNES, finite differences, Jacobian 30 31 .seealso: SNESSetJacobian(), SNESTestJacobian() 32 @*/ 33 int SNESDefaultComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B,int *flag, 34 void *ctx) 35 { 36 Vec j1,j2,x2; 37 int i,ierr,N,start,end,j; 38 Scalar dx, mone = -1.0,*y,scale,*xx; 39 double epsilon = 1.e-8,amax; /* assumes double precision */ 40 41 MatZeroEntries(*J); 42 ierr = VecDuplicate(x1,&j1); CHKERR(ierr); 43 ierr = VecDuplicate(x1,&j2); CHKERR(ierr); 44 ierr = VecDuplicate(x1,&x2); CHKERR(ierr); 45 46 ierr = VecGetSize(x1,&N); CHKERR(ierr); 47 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 48 VecGetArray(x1,&xx); 49 ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 50 for ( i=0; i<N; i++ ) { 51 ierr = VecCopy(x1,x2); CHKERR(ierr); 52 if ( i>= start && i<end) { 53 dx = xx[i-start]; 54 if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-1; 55 else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 56 dx *= epsilon; 57 scale = -1.0/dx; 58 VecSetValues(x2,1,&i,&dx,ADDVALUES); 59 } 60 ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 61 ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 62 VecScale(&scale,j2); 63 VecGetArray(j2,&y); 64 VecAMax(j2,0,&amax); amax *= 1.e-14; 65 for ( j=start; j<end; j++ ) { 66 if (y[j-start] > amax || y[j-start] < -amax) { 67 ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERR(ierr); 68 } 69 } 70 VecRestoreArray(j2,&y); 71 } 72 MatAssemblyBegin(*J,FINAL_ASSEMBLY); 73 VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 74 MatAssemblyEnd(*J,FINAL_ASSEMBLY); 75 return 0; 76 } 77 78