111320018SBarry Smith 2*06be10caSBarry Smith 311320018SBarry Smith #ifndef lint 4*06be10caSBarry Smith static char vcid[] = "$Id: snesj.c,v 1.9 1995/05/11 17:53:49 curfman Exp bsmith $"; 511320018SBarry Smith #endif 611320018SBarry Smith 711320018SBarry Smith #include "draw.h" 839e2f89bSBarry Smith #include "snes.h" 911320018SBarry Smith #include "options.h" 1011320018SBarry Smith 1111320018SBarry Smith /*@ 125f3c43d9SLois Curfman McInnes SNESDefaultComputeJacobian - Computes the Jacobian using finite 135f3c43d9SLois Curfman McInnes differences. 1411320018SBarry Smith 1511320018SBarry Smith Input Parameters: 1623242f5aSBarry Smith . x - compute Jacobian at this point 175f3c43d9SLois Curfman McInnes . ctx - application's Function context 1811320018SBarry Smith 1911320018SBarry Smith Output Parameters: 2023242f5aSBarry Smith . J - Jacobian 2123242f5aSBarry Smith . B - preconditioner, same as Jacobian 2211320018SBarry Smith 235f3c43d9SLois Curfman McInnes Notes: 245f3c43d9SLois Curfman McInnes This routine is slow and expensive, and is not currently optimized 255f3c43d9SLois Curfman McInnes to take advantage of sparsity in the problem. Although 265f3c43d9SLois Curfman McInnes SNESDefaultComputeJacobian() is not recommended for general use 275f3c43d9SLois Curfman McInnes in large-scale applications, It can be useful in checking the 285f3c43d9SLois Curfman McInnes correctness of a user-provided Jacobian. 2911320018SBarry Smith 305f3c43d9SLois Curfman McInnes .keywords: SNES, finite differences, Jacobian 315f3c43d9SLois Curfman McInnes 325f3c43d9SLois Curfman McInnes .seealso: SNESSetJacobian(), SNESTestJacobian() 3311320018SBarry Smith @*/ 3423242f5aSBarry Smith int SNESDefaultComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B,int *flag, 3523242f5aSBarry Smith void *ctx) 3611320018SBarry Smith { 3723242f5aSBarry Smith Vec j1,j2,x2; 3823242f5aSBarry Smith int i,ierr,N,start,end,j; 3939e2f89bSBarry Smith Scalar dx, mone = -1.0,*y,scale,*xx; 40336a5e98SBarry Smith double epsilon = 1.e-8,amax; /* assumes double precision */ 4123242f5aSBarry Smith 42336a5e98SBarry Smith MatZeroEntries(*J); 436469c4f9SBarry Smith ierr = VecDuplicate(x1,&j1); CHKERR(ierr); 446469c4f9SBarry Smith ierr = VecDuplicate(x1,&j2); CHKERR(ierr); 456469c4f9SBarry Smith ierr = VecDuplicate(x1,&x2); CHKERR(ierr); 4623242f5aSBarry Smith 4723242f5aSBarry Smith ierr = VecGetSize(x1,&N); CHKERR(ierr); 4823242f5aSBarry Smith ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 4939e2f89bSBarry Smith VecGetArray(x1,&xx); 5023242f5aSBarry Smith ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 5139e2f89bSBarry Smith for ( i=0; i<N; i++ ) { 52336a5e98SBarry Smith ierr = VecCopy(x1,x2); CHKERR(ierr); 5323242f5aSBarry Smith if ( i>= start && i<end) { 5439e2f89bSBarry Smith dx = xx[i-start]; 55336a5e98SBarry Smith if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-1; 56336a5e98SBarry Smith else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 5739e2f89bSBarry Smith dx *= epsilon; 5839e2f89bSBarry Smith scale = -1.0/dx; 599d00d63dSBarry Smith VecSetValues(x2,1,&i,&dx,ADDVALUES); 6011320018SBarry Smith } 6123242f5aSBarry Smith ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 6223242f5aSBarry Smith ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 6323242f5aSBarry Smith VecScale(&scale,j2); 6423242f5aSBarry Smith VecGetArray(j2,&y); 65336a5e98SBarry Smith VecAMax(j2,0,&amax); amax *= 1.e-14; 6623242f5aSBarry Smith for ( j=start; j<end; j++ ) { 67336a5e98SBarry Smith if (y[j-start] > amax || y[j-start] < -amax) { 68edae2e7dSBarry Smith ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERR(ierr); 6923242f5aSBarry Smith } 7023242f5aSBarry Smith } 7123242f5aSBarry Smith VecRestoreArray(j2,&y); 7223242f5aSBarry Smith } 73ee50ffe9SBarry Smith MatAssemblyBegin(*J,FINAL_ASSEMBLY); 7423242f5aSBarry Smith VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 75ee50ffe9SBarry Smith MatAssemblyEnd(*J,FINAL_ASSEMBLY); 7611320018SBarry Smith return 0; 7711320018SBarry Smith } 7811320018SBarry Smith 79