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