111320018SBarry Smith 211320018SBarry Smith #ifndef lint 3*5f3c43d9SLois Curfman McInnes static char vcid[] = "$Id: snesj.c,v 1.8 1995/05/11 17:44:07 bsmith Exp curfman $"; 411320018SBarry Smith #endif 511320018SBarry Smith 611320018SBarry Smith #include "draw.h" 739e2f89bSBarry Smith #include "snes.h" 811320018SBarry Smith #include "options.h" 911320018SBarry Smith 1011320018SBarry Smith /*@ 11*5f3c43d9SLois Curfman McInnes SNESDefaultComputeJacobian - Computes the Jacobian using finite 12*5f3c43d9SLois Curfman McInnes differences. 1311320018SBarry Smith 1411320018SBarry Smith Input Parameters: 1523242f5aSBarry Smith . x - compute Jacobian at this point 16*5f3c43d9SLois Curfman McInnes . ctx - application's Function context 1711320018SBarry Smith 1811320018SBarry Smith Output Parameters: 1923242f5aSBarry Smith . J - Jacobian 2023242f5aSBarry Smith . B - preconditioner, same as Jacobian 2111320018SBarry Smith 22*5f3c43d9SLois Curfman McInnes Notes: 23*5f3c43d9SLois Curfman McInnes This routine is slow and expensive, and is not currently optimized 24*5f3c43d9SLois Curfman McInnes to take advantage of sparsity in the problem. Although 25*5f3c43d9SLois Curfman McInnes SNESDefaultComputeJacobian() is not recommended for general use 26*5f3c43d9SLois Curfman McInnes in large-scale applications, It can be useful in checking the 27*5f3c43d9SLois Curfman McInnes correctness of a user-provided Jacobian. 2811320018SBarry Smith 29*5f3c43d9SLois Curfman McInnes .keywords: SNES, finite differences, Jacobian 30*5f3c43d9SLois Curfman McInnes 31*5f3c43d9SLois Curfman McInnes .seealso: SNESSetJacobian(), SNESTestJacobian() 3211320018SBarry Smith @*/ 3323242f5aSBarry Smith int SNESDefaultComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B,int *flag, 3423242f5aSBarry Smith void *ctx) 3511320018SBarry Smith { 3623242f5aSBarry Smith Vec j1,j2,x2; 3723242f5aSBarry Smith int i,ierr,N,start,end,j; 3839e2f89bSBarry Smith Scalar dx, mone = -1.0,*y,scale,*xx; 39336a5e98SBarry Smith double epsilon = 1.e-8,amax; /* assumes double precision */ 4023242f5aSBarry Smith 41336a5e98SBarry Smith MatZeroEntries(*J); 426469c4f9SBarry Smith ierr = VecDuplicate(x1,&j1); CHKERR(ierr); 436469c4f9SBarry Smith ierr = VecDuplicate(x1,&j2); CHKERR(ierr); 446469c4f9SBarry Smith ierr = VecDuplicate(x1,&x2); CHKERR(ierr); 4523242f5aSBarry Smith 4623242f5aSBarry Smith ierr = VecGetSize(x1,&N); CHKERR(ierr); 4723242f5aSBarry Smith ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 4839e2f89bSBarry Smith VecGetArray(x1,&xx); 4923242f5aSBarry Smith ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 5039e2f89bSBarry Smith for ( i=0; i<N; i++ ) { 51336a5e98SBarry Smith ierr = VecCopy(x1,x2); CHKERR(ierr); 5223242f5aSBarry Smith if ( i>= start && i<end) { 5339e2f89bSBarry Smith dx = xx[i-start]; 54336a5e98SBarry Smith if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-1; 55336a5e98SBarry Smith else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 5639e2f89bSBarry Smith dx *= epsilon; 5739e2f89bSBarry Smith scale = -1.0/dx; 589d00d63dSBarry Smith VecSetValues(x2,1,&i,&dx,ADDVALUES); 5911320018SBarry Smith } 6023242f5aSBarry Smith ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 6123242f5aSBarry Smith ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 6223242f5aSBarry Smith VecScale(&scale,j2); 6323242f5aSBarry Smith VecGetArray(j2,&y); 64336a5e98SBarry Smith VecAMax(j2,0,&amax); amax *= 1.e-14; 6523242f5aSBarry Smith for ( j=start; j<end; j++ ) { 66336a5e98SBarry Smith if (y[j-start] > amax || y[j-start] < -amax) { 67edae2e7dSBarry Smith ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERR(ierr); 6823242f5aSBarry Smith } 6923242f5aSBarry Smith } 7023242f5aSBarry Smith VecRestoreArray(j2,&y); 7123242f5aSBarry Smith } 72ee50ffe9SBarry Smith MatAssemblyBegin(*J,FINAL_ASSEMBLY); 7323242f5aSBarry Smith VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 74ee50ffe9SBarry Smith MatAssemblyEnd(*J,FINAL_ASSEMBLY); 7511320018SBarry Smith return 0; 7611320018SBarry Smith } 7711320018SBarry Smith 78