111320018SBarry Smith 206be10caSBarry Smith 311320018SBarry Smith #ifndef lint 4*da88328cSLois Curfman McInnes static char vcid[] = "$Id: snesj.c,v 1.10 1995/05/12 04:18:43 bsmith Exp curfman $"; 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: 16*da88328cSLois Curfman McInnes . x1 - compute Jacobian at this point 17*da88328cSLois Curfman McInnes . ctx - application's function context, as set with SNESSetFunction() 1811320018SBarry Smith 1911320018SBarry Smith Output Parameters: 2023242f5aSBarry Smith . J - Jacobian 2123242f5aSBarry Smith . B - preconditioner, same as Jacobian 22*da88328cSLois Curfman McInnes . flag - matrix flag 2311320018SBarry Smith 245f3c43d9SLois Curfman McInnes Notes: 255f3c43d9SLois Curfman McInnes This routine is slow and expensive, and is not currently optimized 265f3c43d9SLois Curfman McInnes to take advantage of sparsity in the problem. Although 275f3c43d9SLois Curfman McInnes SNESDefaultComputeJacobian() is not recommended for general use 285f3c43d9SLois Curfman McInnes in large-scale applications, It can be useful in checking the 295f3c43d9SLois Curfman McInnes correctness of a user-provided Jacobian. 3011320018SBarry Smith 315f3c43d9SLois Curfman McInnes .keywords: SNES, finite differences, Jacobian 325f3c43d9SLois Curfman McInnes 335f3c43d9SLois Curfman McInnes .seealso: SNESSetJacobian(), SNESTestJacobian() 3411320018SBarry Smith @*/ 3523242f5aSBarry Smith int SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,int *flag, 3623242f5aSBarry Smith void *ctx) 3711320018SBarry Smith { 3823242f5aSBarry Smith Vec j1,j2,x2; 3923242f5aSBarry Smith int i,ierr,N,start,end,j; 4039e2f89bSBarry Smith Scalar dx, mone = -1.0,*y,scale,*xx; 41336a5e98SBarry Smith double epsilon = 1.e-8,amax; /* assumes double precision */ 4223242f5aSBarry Smith 43336a5e98SBarry Smith MatZeroEntries(*J); 446469c4f9SBarry Smith ierr = VecDuplicate(x1,&j1); CHKERR(ierr); 456469c4f9SBarry Smith ierr = VecDuplicate(x1,&j2); CHKERR(ierr); 466469c4f9SBarry Smith ierr = VecDuplicate(x1,&x2); CHKERR(ierr); 4723242f5aSBarry Smith 4823242f5aSBarry Smith ierr = VecGetSize(x1,&N); CHKERR(ierr); 4923242f5aSBarry Smith ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 5039e2f89bSBarry Smith VecGetArray(x1,&xx); 5123242f5aSBarry Smith ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 5239e2f89bSBarry Smith for ( i=0; i<N; i++ ) { 53336a5e98SBarry Smith ierr = VecCopy(x1,x2); CHKERR(ierr); 5423242f5aSBarry Smith if ( i>= start && i<end) { 5539e2f89bSBarry Smith dx = xx[i-start]; 56336a5e98SBarry Smith if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-1; 57336a5e98SBarry Smith else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 5839e2f89bSBarry Smith dx *= epsilon; 5939e2f89bSBarry Smith scale = -1.0/dx; 609d00d63dSBarry Smith VecSetValues(x2,1,&i,&dx,ADDVALUES); 6111320018SBarry Smith } 6223242f5aSBarry Smith ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 6323242f5aSBarry Smith ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 6423242f5aSBarry Smith VecScale(&scale,j2); 6523242f5aSBarry Smith VecGetArray(j2,&y); 66336a5e98SBarry Smith VecAMax(j2,0,&amax); amax *= 1.e-14; 6723242f5aSBarry Smith for ( j=start; j<end; j++ ) { 68336a5e98SBarry Smith if (y[j-start] > amax || y[j-start] < -amax) { 69edae2e7dSBarry Smith ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERR(ierr); 7023242f5aSBarry Smith } 7123242f5aSBarry Smith } 7223242f5aSBarry Smith VecRestoreArray(j2,&y); 7323242f5aSBarry Smith } 74ee50ffe9SBarry Smith MatAssemblyBegin(*J,FINAL_ASSEMBLY); 7523242f5aSBarry Smith VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 76ee50ffe9SBarry Smith MatAssemblyEnd(*J,FINAL_ASSEMBLY); 7711320018SBarry Smith return 0; 7811320018SBarry Smith } 7911320018SBarry Smith 80