111320018SBarry Smith 211320018SBarry Smith #ifndef lint 3*78b31e54SBarry Smith static char vcid[] = "$Id: snesj.c,v 1.14 1995/05/16 00:36:51 curfman Exp bsmith $"; 411320018SBarry Smith #endif 511320018SBarry Smith 611320018SBarry Smith #include "draw.h" 739e2f89bSBarry Smith #include "snes.h" 811320018SBarry Smith 911320018SBarry Smith /*@ 105f3c43d9SLois Curfman McInnes SNESDefaultComputeJacobian - Computes the Jacobian using finite 115f3c43d9SLois Curfman McInnes differences. 1211320018SBarry Smith 1311320018SBarry Smith Input Parameters: 14da88328cSLois Curfman McInnes . x1 - compute Jacobian at this point 15da88328cSLois Curfman McInnes . ctx - application's function context, as set with SNESSetFunction() 1611320018SBarry Smith 1711320018SBarry Smith Output Parameters: 1823242f5aSBarry Smith . J - Jacobian 1923242f5aSBarry Smith . B - preconditioner, same as Jacobian 20da88328cSLois Curfman McInnes . flag - matrix flag 2111320018SBarry Smith 22ad960d00SLois Curfman McInnes Options Database Key: 23ad960d00SLois Curfman McInnes $ -snes_fd 24ad960d00SLois Curfman McInnes 255f3c43d9SLois Curfman McInnes Notes: 265f3c43d9SLois Curfman McInnes This routine is slow and expensive, and is not currently optimized 275f3c43d9SLois Curfman McInnes to take advantage of sparsity in the problem. Although 285f3c43d9SLois Curfman McInnes SNESDefaultComputeJacobian() is not recommended for general use 295f3c43d9SLois Curfman McInnes in large-scale applications, It can be useful in checking the 305f3c43d9SLois Curfman McInnes correctness of a user-provided Jacobian. 3111320018SBarry Smith 325f3c43d9SLois Curfman McInnes .keywords: SNES, finite differences, Jacobian 335f3c43d9SLois Curfman McInnes 345f3c43d9SLois Curfman McInnes .seealso: SNESSetJacobian(), SNESTestJacobian() 3511320018SBarry Smith @*/ 36eccfb7ebSLois Curfman McInnes int SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B, 37eccfb7ebSLois Curfman McInnes MatStructure *flag,void *ctx) 3811320018SBarry Smith { 3923242f5aSBarry Smith Vec j1,j2,x2; 4023242f5aSBarry Smith int i,ierr,N,start,end,j; 4139e2f89bSBarry Smith Scalar dx, mone = -1.0,*y,scale,*xx; 42336a5e98SBarry Smith double epsilon = 1.e-8,amax; /* assumes double precision */ 4323242f5aSBarry Smith 44336a5e98SBarry Smith MatZeroEntries(*J); 45*78b31e54SBarry Smith ierr = VecDuplicate(x1,&j1); CHKERRQ(ierr); 46*78b31e54SBarry Smith ierr = VecDuplicate(x1,&j2); CHKERRQ(ierr); 47*78b31e54SBarry Smith ierr = VecDuplicate(x1,&x2); CHKERRQ(ierr); 4823242f5aSBarry Smith 49*78b31e54SBarry Smith ierr = VecGetSize(x1,&N); CHKERRQ(ierr); 50*78b31e54SBarry Smith ierr = VecGetOwnershipRange(x1,&start,&end); CHKERRQ(ierr); 5139e2f89bSBarry Smith VecGetArray(x1,&xx); 52*78b31e54SBarry Smith ierr = SNESComputeFunction(snes,x1,j1); CHKERRQ(ierr); 5339e2f89bSBarry Smith for ( i=0; i<N; i++ ) { 54*78b31e54SBarry Smith ierr = VecCopy(x1,x2); CHKERRQ(ierr); 5523242f5aSBarry Smith if ( i>= start && i<end) { 5639e2f89bSBarry Smith dx = xx[i-start]; 57336a5e98SBarry Smith if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-1; 58336a5e98SBarry Smith else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 5939e2f89bSBarry Smith dx *= epsilon; 6039e2f89bSBarry Smith scale = -1.0/dx; 619d00d63dSBarry Smith VecSetValues(x2,1,&i,&dx,ADDVALUES); 6211320018SBarry Smith } 63*78b31e54SBarry Smith ierr = SNESComputeFunction(snes,x2,j2); CHKERRQ(ierr); 64*78b31e54SBarry Smith ierr = VecAXPY(&mone,j1,j2); CHKERRQ(ierr); 6523242f5aSBarry Smith VecScale(&scale,j2); 6623242f5aSBarry Smith VecGetArray(j2,&y); 67336a5e98SBarry Smith VecAMax(j2,0,&amax); amax *= 1.e-14; 6823242f5aSBarry Smith for ( j=start; j<end; j++ ) { 69336a5e98SBarry Smith if (y[j-start] > amax || y[j-start] < -amax) { 70*78b31e54SBarry Smith ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERRQ(ierr); 7123242f5aSBarry Smith } 7223242f5aSBarry Smith } 7323242f5aSBarry Smith VecRestoreArray(j2,&y); 7423242f5aSBarry Smith } 75ee50ffe9SBarry Smith MatAssemblyBegin(*J,FINAL_ASSEMBLY); 7623242f5aSBarry Smith VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 77ee50ffe9SBarry Smith MatAssemblyEnd(*J,FINAL_ASSEMBLY); 7811320018SBarry Smith return 0; 7911320018SBarry Smith } 8011320018SBarry Smith 81