111320018SBarry Smith 211320018SBarry Smith #ifndef lint 3*9d00d63dSBarry Smith static char vcid[] = "$Id: snesj.c,v 1.4 1995/05/02 23:39:40 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 /*@ 1123242f5aSBarry Smith SNESDefaultComputeJacobian - Computes Jacobian using finite 1223242f5aSBarry Smith differences. Slow and expensive. 1311320018SBarry Smith 1411320018SBarry Smith Input Parameters: 1523242f5aSBarry Smith . x - compute Jacobian at this point 1623242f5aSBarry Smith . ctx - applications Function context 1711320018SBarry Smith 1811320018SBarry Smith Output Parameters: 1923242f5aSBarry Smith . J - Jacobian 2023242f5aSBarry Smith . B - preconditioner, same as Jacobian 2111320018SBarry Smith 2223242f5aSBarry Smith .keywords: finite differences, Jacobian 2311320018SBarry Smith 2423242f5aSBarry Smith .seealso: SNESSetJacobian, SNESTestJacobian 2511320018SBarry Smith @*/ 2623242f5aSBarry Smith int SNESDefaultComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B,int *flag, 2723242f5aSBarry Smith void *ctx) 2811320018SBarry Smith { 2923242f5aSBarry Smith Vec j1,j2,x2; 3023242f5aSBarry Smith int i,ierr,N,start,end,j; 3123242f5aSBarry Smith Scalar dx, mone = -1.0,*y,scale; 3211320018SBarry Smith double norm; 3323242f5aSBarry Smith 3423242f5aSBarry Smith ierr = VecCreate(x1,&j1); CHKERR(ierr); 3523242f5aSBarry Smith ierr = VecCreate(x1,&j2); CHKERR(ierr); 3623242f5aSBarry Smith ierr = VecCreate(x1,&x2); CHKERR(ierr); 3723242f5aSBarry Smith ierr = VecNorm(x1,&norm); CHKERR(ierr); 3823242f5aSBarry Smith dx = 1.e-8*norm; /* assumes double precision */ 3923242f5aSBarry Smith scale = -1.0/dx; 4023242f5aSBarry Smith 4123242f5aSBarry Smith ierr = VecGetSize(x1,&N); CHKERR(ierr); 4223242f5aSBarry Smith ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 4323242f5aSBarry Smith ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 4423242f5aSBarry Smith for ( i=0; i<N; i++ ) { 4523242f5aSBarry Smith ierr = VecCopy(x1,x2); CHKERR(ierr); 4623242f5aSBarry Smith if ( i>= start && i<end) { 47*9d00d63dSBarry Smith VecSetValues(x2,1,&i,&dx,ADDVALUES); 4811320018SBarry Smith } 4923242f5aSBarry Smith ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 5023242f5aSBarry Smith ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 5123242f5aSBarry Smith VecScale(&scale,j2); 5223242f5aSBarry Smith VecGetArray(j2,&y); 5323242f5aSBarry Smith for ( j=start; j<end; j++ ) { 5423242f5aSBarry Smith if (y[j-start]) { 55edae2e7dSBarry Smith ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERR(ierr); 5623242f5aSBarry Smith } 5723242f5aSBarry Smith } 5823242f5aSBarry Smith VecRestoreArray(j2,&y); 5923242f5aSBarry Smith } 60ee50ffe9SBarry Smith MatAssemblyBegin(*J,FINAL_ASSEMBLY); 6123242f5aSBarry Smith VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 62ee50ffe9SBarry Smith MatAssemblyEnd(*J,FINAL_ASSEMBLY); 6311320018SBarry Smith return 0; 6411320018SBarry Smith } 6511320018SBarry Smith 6611320018SBarry Smith /*@ 6723242f5aSBarry Smith SNESTestJacobian - Tests whether a hand computed Jacobian 6823242f5aSBarry Smith matches one compute via finite differences 6911320018SBarry Smith 7011320018SBarry Smith Input Parameters: 7111320018SBarry Smith 7223242f5aSBarry Smith Output Parameters: 7323242f5aSBarry Smith 7411320018SBarry Smith @*/ 7523242f5aSBarry Smith int SNESTestJacobian(SNES snes) 7611320018SBarry Smith { 7711320018SBarry Smith 7823242f5aSBarry Smith /* compute both versions of Jacobian */ 7911320018SBarry Smith 8023242f5aSBarry Smith /* compare */ 8111320018SBarry Smith 8211320018SBarry Smith return 0; 8311320018SBarry Smith } 84