1 2 3 #ifndef lint 4 static char vcid[] = "$Id: snesj.c,v 1.6 1995/05/03 13:21:44 bsmith Exp bsmith $"; 5 #endif 6 7 #include "draw.h" 8 #include "snes.h" 9 #include "options.h" 10 11 /*@ 12 SNESDefaultComputeJacobian - Computes Jacobian using finite 13 differences. Slow and expensive. 14 15 Input Parameters: 16 . x - compute Jacobian at this point 17 . ctx - applications Function context 18 19 Output Parameters: 20 . J - Jacobian 21 . B - preconditioner, same as Jacobian 22 23 .keywords: finite differences, Jacobian 24 25 .seealso: SNESSetJacobian, SNESTestJacobian 26 @*/ 27 int SNESDefaultComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B,int *flag, 28 void *ctx) 29 { 30 Vec j1,j2,x2; 31 int i,ierr,N,start,end,j; 32 Scalar dx, mone = -1.0,*y,scale,*xx; 33 double epsilon = 1.e-8; /* assumes double precision */ 34 35 ierr = VecDuplicate(x1,&j1); CHKERR(ierr); 36 ierr = VecDuplicate(x1,&j2); CHKERR(ierr); 37 ierr = VecDuplicate(x1,&x2); CHKERR(ierr); 38 39 ierr = VecGetSize(x1,&N); CHKERR(ierr); 40 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 41 VecGetArray(x1,&xx); 42 ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 43 ierr = VecCopy(x1,x2); CHKERR(ierr); 44 for ( i=0; i<N; i++ ) { 45 if ( i>= start && i<end) { 46 dx = xx[i-start]; 47 if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-16; 48 else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-16; 49 dx *= epsilon; 50 scale = -1.0/dx; 51 VecSetValues(x2,1,&i,&dx,ADDVALUES); 52 } 53 ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 54 if ( i>= start && i<end) { 55 dx = -dx; 56 VecSetValues(x2,1,&i,&dx,ADDVALUES); 57 } 58 ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 59 VecScale(&scale,j2); 60 VecGetArray(j2,&y); 61 for ( j=start; j<end; j++ ) { 62 if (y[j-start]) { 63 ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERR(ierr); 64 } 65 } 66 VecRestoreArray(j2,&y); 67 } 68 MatAssemblyBegin(*J,FINAL_ASSEMBLY); 69 VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 70 MatAssemblyEnd(*J,FINAL_ASSEMBLY); 71 return 0; 72 } 73 74 /*@ 75 SNESTestJacobian - Tests whether a hand computed Jacobian 76 matches one compute via finite differences 77 78 Input Parameters: 79 80 Output Parameters: 81 82 @*/ 83 int SNESTestJacobian(SNES snes) 84 { 85 86 /* compute both versions of Jacobian */ 87 88 /* compare */ 89 90 return 0; 91 } 92