111320018SBarry Smith 206be10caSBarry Smith 311320018SBarry Smith #ifndef lint 4*ad960d00SLois Curfman McInnes static char vcid[] = "$Id: snesj.c,v 1.11 1995/05/12 18:16:21 curfman 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: 16da88328cSLois Curfman McInnes . x1 - compute Jacobian at this point 17da88328cSLois 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 22da88328cSLois Curfman McInnes . flag - matrix flag 2311320018SBarry Smith 24*ad960d00SLois Curfman McInnes Options Database Key: 25*ad960d00SLois Curfman McInnes $ -snes_fd 26*ad960d00SLois Curfman McInnes 275f3c43d9SLois Curfman McInnes Notes: 285f3c43d9SLois Curfman McInnes This routine is slow and expensive, and is not currently optimized 295f3c43d9SLois Curfman McInnes to take advantage of sparsity in the problem. Although 305f3c43d9SLois Curfman McInnes SNESDefaultComputeJacobian() is not recommended for general use 315f3c43d9SLois Curfman McInnes in large-scale applications, It can be useful in checking the 325f3c43d9SLois Curfman McInnes correctness of a user-provided Jacobian. 3311320018SBarry Smith 345f3c43d9SLois Curfman McInnes .keywords: SNES, finite differences, Jacobian 355f3c43d9SLois Curfman McInnes 365f3c43d9SLois Curfman McInnes .seealso: SNESSetJacobian(), SNESTestJacobian() 3711320018SBarry Smith @*/ 3823242f5aSBarry Smith int SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,int *flag, 3923242f5aSBarry Smith void *ctx) 4011320018SBarry Smith { 4123242f5aSBarry Smith Vec j1,j2,x2; 4223242f5aSBarry Smith int i,ierr,N,start,end,j; 4339e2f89bSBarry Smith Scalar dx, mone = -1.0,*y,scale,*xx; 44336a5e98SBarry Smith double epsilon = 1.e-8,amax; /* assumes double precision */ 4523242f5aSBarry Smith 46336a5e98SBarry Smith MatZeroEntries(*J); 476469c4f9SBarry Smith ierr = VecDuplicate(x1,&j1); CHKERR(ierr); 486469c4f9SBarry Smith ierr = VecDuplicate(x1,&j2); CHKERR(ierr); 496469c4f9SBarry Smith ierr = VecDuplicate(x1,&x2); CHKERR(ierr); 5023242f5aSBarry Smith 5123242f5aSBarry Smith ierr = VecGetSize(x1,&N); CHKERR(ierr); 5223242f5aSBarry Smith ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 5339e2f89bSBarry Smith VecGetArray(x1,&xx); 5423242f5aSBarry Smith ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 5539e2f89bSBarry Smith for ( i=0; i<N; i++ ) { 56336a5e98SBarry Smith ierr = VecCopy(x1,x2); CHKERR(ierr); 5723242f5aSBarry Smith if ( i>= start && i<end) { 5839e2f89bSBarry Smith dx = xx[i-start]; 59336a5e98SBarry Smith if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-1; 60336a5e98SBarry Smith else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 6139e2f89bSBarry Smith dx *= epsilon; 6239e2f89bSBarry Smith scale = -1.0/dx; 639d00d63dSBarry Smith VecSetValues(x2,1,&i,&dx,ADDVALUES); 6411320018SBarry Smith } 6523242f5aSBarry Smith ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 6623242f5aSBarry Smith ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 6723242f5aSBarry Smith VecScale(&scale,j2); 6823242f5aSBarry Smith VecGetArray(j2,&y); 69336a5e98SBarry Smith VecAMax(j2,0,&amax); amax *= 1.e-14; 7023242f5aSBarry Smith for ( j=start; j<end; j++ ) { 71336a5e98SBarry Smith if (y[j-start] > amax || y[j-start] < -amax) { 72edae2e7dSBarry Smith ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERR(ierr); 7323242f5aSBarry Smith } 7423242f5aSBarry Smith } 7523242f5aSBarry Smith VecRestoreArray(j2,&y); 7623242f5aSBarry Smith } 77ee50ffe9SBarry Smith MatAssemblyBegin(*J,FINAL_ASSEMBLY); 7823242f5aSBarry Smith VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 79ee50ffe9SBarry Smith MatAssemblyEnd(*J,FINAL_ASSEMBLY); 8011320018SBarry Smith return 0; 8111320018SBarry Smith } 8211320018SBarry Smith 83