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