1 2 3 #ifndef lint 4 static char vcid[] = "$Id: snesj.c,v 1.12 1995/05/12 18:17:04 curfman Exp curfman $"; 5 #endif 6 7 #include "draw.h" 8 #include "snes.h" 9 #include "options.h" 10 11 /*@ 12 SNESDefaultComputeJacobian - Computes the Jacobian using finite 13 differences. 14 15 Input Parameters: 16 . x1 - compute Jacobian at this point 17 . ctx - application's function context, as set with SNESSetFunction() 18 19 Output Parameters: 20 . J - Jacobian 21 . B - preconditioner, same as Jacobian 22 . flag - matrix flag 23 24 Options Database Key: 25 $ -snes_fd 26 27 Notes: 28 This routine is slow and expensive, and is not currently optimized 29 to take advantage of sparsity in the problem. Although 30 SNESDefaultComputeJacobian() is not recommended for general use 31 in large-scale applications, It can be useful in checking the 32 correctness of a user-provided Jacobian. 33 34 .keywords: SNES, finite differences, Jacobian 35 36 .seealso: SNESSetJacobian(), SNESTestJacobian() 37 @*/ 38 int SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B, 39 MatStructure *flag,void *ctx) 40 { 41 Vec j1,j2,x2; 42 int i,ierr,N,start,end,j; 43 Scalar dx, mone = -1.0,*y,scale,*xx; 44 double epsilon = 1.e-8,amax; /* assumes double precision */ 45 46 MatZeroEntries(*J); 47 ierr = VecDuplicate(x1,&j1); CHKERR(ierr); 48 ierr = VecDuplicate(x1,&j2); CHKERR(ierr); 49 ierr = VecDuplicate(x1,&x2); CHKERR(ierr); 50 51 ierr = VecGetSize(x1,&N); CHKERR(ierr); 52 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERR(ierr); 53 VecGetArray(x1,&xx); 54 ierr = SNESComputeFunction(snes,x1,j1); CHKERR(ierr); 55 for ( i=0; i<N; i++ ) { 56 ierr = VecCopy(x1,x2); CHKERR(ierr); 57 if ( i>= start && i<end) { 58 dx = xx[i-start]; 59 if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-1; 60 else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 61 dx *= epsilon; 62 scale = -1.0/dx; 63 VecSetValues(x2,1,&i,&dx,ADDVALUES); 64 } 65 ierr = SNESComputeFunction(snes,x2,j2); CHKERR(ierr); 66 ierr = VecAXPY(&mone,j1,j2); CHKERR(ierr); 67 VecScale(&scale,j2); 68 VecGetArray(j2,&y); 69 VecAMax(j2,0,&amax); amax *= 1.e-14; 70 for ( j=start; j<end; j++ ) { 71 if (y[j-start] > amax || y[j-start] < -amax) { 72 ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERTVALUES); CHKERR(ierr); 73 } 74 } 75 VecRestoreArray(j2,&y); 76 } 77 MatAssemblyBegin(*J,FINAL_ASSEMBLY); 78 VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 79 MatAssemblyEnd(*J,FINAL_ASSEMBLY); 80 return 0; 81 } 82 83