111320018SBarry Smith 211320018SBarry Smith #ifndef lint 3*23242f5aSBarry Smith static char vcid[] = "$Id: snesj.c,v 1.1 1995/05/02 02:14:58 bsmith Exp bsmith $"; 411320018SBarry Smith #endif 511320018SBarry Smith 611320018SBarry Smith #include "draw.h" 711320018SBarry Smith #include "snesimpl.h" 811320018SBarry Smith #include "options.h" 911320018SBarry Smith 1011320018SBarry Smith /*@ 11*23242f5aSBarry Smith SNESDefaultComputeJacobian - Computes Jacobian using finite 12*23242f5aSBarry Smith differences. Slow and expensive. 1311320018SBarry Smith 1411320018SBarry Smith Input Parameters: 15*23242f5aSBarry Smith . x - compute Jacobian at this point 16*23242f5aSBarry Smith . ctx - applications Function context 1711320018SBarry Smith 1811320018SBarry Smith Output Parameters: 19*23242f5aSBarry Smith . J - Jacobian 20*23242f5aSBarry Smith . B - preconditioner, same as Jacobian 2111320018SBarry Smith 22*23242f5aSBarry Smith .keywords: finite differences, Jacobian 2311320018SBarry Smith 24*23242f5aSBarry Smith .seealso: SNESSetJacobian, SNESTestJacobian 2511320018SBarry Smith @*/ 26*23242f5aSBarry Smith int SNESDefaultComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B,int *flag, 27*23242f5aSBarry Smith void *ctx) 2811320018SBarry Smith { 29*23242f5aSBarry Smith Vec j1,j2,x2; 30*23242f5aSBarry Smith int i,ierr,N,start,end,j; 31*23242f5aSBarry Smith Scalar dx, mone = -1.0,*y,scale; 3211320018SBarry Smith double norm; 33*23242f5aSBarry Smith 34*23242f5aSBarry Smith ierr = VecCreate(x1,&j1); CHKERR(ierr); 35*23242f5aSBarry Smith ierr = VecCreate(x1,&j2); CHKERR(ierr); 36*23242f5aSBarry Smith ierr = VecCreate(x1,&x2); CHKERR(ierr); 37*23242f5aSBarry Smith ierr = VecNorm(x1,&norm); CHKERR(ierr); 38*23242f5aSBarry Smith dx = 1.e-8*norm; /* assumes double precision */ 39*23242f5aSBarry Smith scale = -1.0/dx; 40*23242f5aSBarry Smith 41*23242f5aSBarry Smith ierr = VecGetSize(x1,&N); CHKERR(ierr); 42*23242f5aSBarry Smith ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 43*23242f5aSBarry Smith ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 44*23242f5aSBarry Smith for ( i=0; i<N; i++ ) { 45*23242f5aSBarry Smith ierr = VecCopy(x1,x2); CHKERR(ierr); 46*23242f5aSBarry Smith if ( i>= start && i<end) { 47*23242f5aSBarry Smith VecSetValues(x2,1,&i,&dx,AddValues); 4811320018SBarry Smith } 49*23242f5aSBarry Smith ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 50*23242f5aSBarry Smith ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 51*23242f5aSBarry Smith VecScale(&scale,j2); 52*23242f5aSBarry Smith VecGetArray(j2,&y); 53*23242f5aSBarry Smith for ( j=start; j<end; j++ ) { 54*23242f5aSBarry Smith if (y[j-start]) { 55*23242f5aSBarry Smith ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,InsertValues); CHKERR(ierr); 56*23242f5aSBarry Smith } 57*23242f5aSBarry Smith } 58*23242f5aSBarry Smith VecRestoreArray(j2,&y); 59*23242f5aSBarry Smith } 60*23242f5aSBarry Smith MatBeginAssembly(*J,FINAL_ASSEMBLY); 61*23242f5aSBarry Smith VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 62*23242f5aSBarry Smith MatEndAssembly(*J,FINAL_ASSEMBLY); 6311320018SBarry Smith return 0; 6411320018SBarry Smith } 6511320018SBarry Smith 6611320018SBarry Smith /*@ 67*23242f5aSBarry Smith SNESTestJacobian - Tests whether a hand computed Jacobian 68*23242f5aSBarry Smith matches one compute via finite differences 6911320018SBarry Smith 7011320018SBarry Smith Input Parameters: 7111320018SBarry Smith 72*23242f5aSBarry Smith Output Parameters: 73*23242f5aSBarry Smith 7411320018SBarry Smith @*/ 75*23242f5aSBarry Smith int SNESTestJacobian(SNES snes) 7611320018SBarry Smith { 7711320018SBarry Smith 78*23242f5aSBarry Smith /* compute both versions of Jacobian */ 7911320018SBarry Smith 80*23242f5aSBarry Smith /* compare */ 8111320018SBarry Smith 8211320018SBarry Smith return 0; 8311320018SBarry Smith } 84