111320018SBarry Smith 2*39e2f89bSBarry Smith 311320018SBarry Smith #ifndef lint 4*39e2f89bSBarry Smith static char vcid[] = "$Id: snesj.c,v 1.6 1995/05/03 13:21:44 bsmith Exp bsmith $"; 511320018SBarry Smith #endif 611320018SBarry Smith 711320018SBarry Smith #include "draw.h" 8*39e2f89bSBarry Smith #include "snes.h" 911320018SBarry Smith #include "options.h" 1011320018SBarry Smith 1111320018SBarry Smith /*@ 1223242f5aSBarry Smith SNESDefaultComputeJacobian - Computes Jacobian using finite 1323242f5aSBarry Smith differences. Slow and expensive. 1411320018SBarry Smith 1511320018SBarry Smith Input Parameters: 1623242f5aSBarry Smith . x - compute Jacobian at this point 1723242f5aSBarry Smith . ctx - applications Function context 1811320018SBarry Smith 1911320018SBarry Smith Output Parameters: 2023242f5aSBarry Smith . J - Jacobian 2123242f5aSBarry Smith . B - preconditioner, same as Jacobian 2211320018SBarry Smith 2323242f5aSBarry Smith .keywords: finite differences, Jacobian 2411320018SBarry Smith 2523242f5aSBarry Smith .seealso: SNESSetJacobian, SNESTestJacobian 2611320018SBarry Smith @*/ 2723242f5aSBarry Smith int SNESDefaultComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B,int *flag, 2823242f5aSBarry Smith void *ctx) 2911320018SBarry Smith { 3023242f5aSBarry Smith Vec j1,j2,x2; 3123242f5aSBarry Smith int i,ierr,N,start,end,j; 32*39e2f89bSBarry Smith Scalar dx, mone = -1.0,*y,scale,*xx; 33*39e2f89bSBarry Smith double epsilon = 1.e-8; /* assumes double precision */ 3423242f5aSBarry Smith 356469c4f9SBarry Smith ierr = VecDuplicate(x1,&j1); CHKERR(ierr); 366469c4f9SBarry Smith ierr = VecDuplicate(x1,&j2); CHKERR(ierr); 376469c4f9SBarry Smith ierr = VecDuplicate(x1,&x2); CHKERR(ierr); 3823242f5aSBarry Smith 3923242f5aSBarry Smith ierr = VecGetSize(x1,&N); CHKERR(ierr); 4023242f5aSBarry Smith ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 41*39e2f89bSBarry Smith VecGetArray(x1,&xx); 4223242f5aSBarry Smith ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 4323242f5aSBarry Smith ierr = VecCopy(x1,x2); CHKERR(ierr); 44*39e2f89bSBarry Smith for ( i=0; i<N; i++ ) { 4523242f5aSBarry Smith if ( i>= start && i<end) { 46*39e2f89bSBarry Smith dx = xx[i-start]; 47*39e2f89bSBarry Smith if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-16; 48*39e2f89bSBarry Smith else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-16; 49*39e2f89bSBarry Smith dx *= epsilon; 50*39e2f89bSBarry Smith scale = -1.0/dx; 519d00d63dSBarry Smith VecSetValues(x2,1,&i,&dx,ADDVALUES); 5211320018SBarry Smith } 5323242f5aSBarry Smith ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 54*39e2f89bSBarry Smith if ( i>= start && i<end) { 55*39e2f89bSBarry Smith dx = -dx; 56*39e2f89bSBarry Smith VecSetValues(x2,1,&i,&dx,ADDVALUES); 57*39e2f89bSBarry Smith } 5823242f5aSBarry Smith ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 5923242f5aSBarry Smith VecScale(&scale,j2); 6023242f5aSBarry Smith VecGetArray(j2,&y); 6123242f5aSBarry Smith for ( j=start; j<end; j++ ) { 6223242f5aSBarry Smith if (y[j-start]) { 63edae2e7dSBarry Smith ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERR(ierr); 6423242f5aSBarry Smith } 6523242f5aSBarry Smith } 6623242f5aSBarry Smith VecRestoreArray(j2,&y); 6723242f5aSBarry Smith } 68ee50ffe9SBarry Smith MatAssemblyBegin(*J,FINAL_ASSEMBLY); 6923242f5aSBarry Smith VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 70ee50ffe9SBarry Smith MatAssemblyEnd(*J,FINAL_ASSEMBLY); 7111320018SBarry Smith return 0; 7211320018SBarry Smith } 7311320018SBarry Smith 7411320018SBarry Smith /*@ 7523242f5aSBarry Smith SNESTestJacobian - Tests whether a hand computed Jacobian 7623242f5aSBarry Smith matches one compute via finite differences 7711320018SBarry Smith 7811320018SBarry Smith Input Parameters: 7911320018SBarry Smith 8023242f5aSBarry Smith Output Parameters: 8123242f5aSBarry Smith 8211320018SBarry Smith @*/ 8323242f5aSBarry Smith int SNESTestJacobian(SNES snes) 8411320018SBarry Smith { 8511320018SBarry Smith 8623242f5aSBarry Smith /* compute both versions of Jacobian */ 8711320018SBarry Smith 8823242f5aSBarry Smith /* compare */ 8911320018SBarry Smith 9011320018SBarry Smith return 0; 9111320018SBarry Smith } 92