1 2 #ifndef lint 3 static char vcid[] = "$Id: snesj.c,v 1.26 1996/02/05 19:32:38 curfman Exp curfman $"; 4 #endif 5 6 #include "draw.h" /*I "draw.h" I*/ 7 #include "snesimpl.h" /*I "snes.h" I*/ 8 9 /*@C 10 SNESDefaultComputeJacobian - Computes the Jacobian using finite 11 differences. 12 13 Input Parameters: 14 . x1 - compute Jacobian at this point 15 . ctx - application's function context, as set with SNESSetFunction() 16 17 Output Parameters: 18 . J - Jacobian 19 . B - preconditioner, same as Jacobian 20 . flag - matrix flag 21 22 Options Database Key: 23 $ -snes_fd 24 25 Notes: 26 This routine is slow and expensive, and is not currently optimized 27 to take advantage of sparsity in the problem. Although 28 SNESDefaultComputeJacobian() is not recommended for general use 29 in large-scale applications, It can be useful in checking the 30 correctness of a user-provided Jacobian. 31 32 .keywords: SNES, finite differences, Jacobian 33 34 .seealso: SNESSetJacobian(), SNESTestJacobian() 35 @*/ 36 int SNESDefaultComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 37 { 38 Vec j1,j2,x2; 39 int i,ierr,N,start,end,j; 40 Scalar dx, mone = -1.0,*y,scale,*xx,wscale; 41 double amax, epsilon = 1.e-8; /* assumes double precision */ 42 MPI_Comm comm; 43 int (*eval_fct)(SNES,Vec,Vec); 44 45 if (snes->method_class == SNES_NONLINEAR_EQUATIONS) 46 eval_fct = SNESComputeFunction; 47 else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) 48 eval_fct = SNESComputeGradient; 49 else SETERRQ(1,"SNESDefaultComputeJacobian: Invalid method class"); 50 51 PetscObjectGetComm((PetscObject)x1,&comm); 52 MatZeroEntries(*J); 53 ierr = VecDuplicate(x1,&j1); CHKERRQ(ierr); 54 ierr = VecDuplicate(x1,&j2); CHKERRQ(ierr); 55 ierr = VecDuplicate(x1,&x2); CHKERRQ(ierr); 56 PLogObjectParent(snes,j1); PLogObjectParent(snes,j2); 57 PLogObjectParent(snes,x2); 58 59 ierr = VecGetSize(x1,&N); CHKERRQ(ierr); 60 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERRQ(ierr); 61 VecGetArray(x1,&xx); 62 ierr = eval_fct(snes,x1,j1); CHKERRQ(ierr); 63 for ( i=0; i<N; i++ ) { 64 ierr = VecCopy(x1,x2); CHKERRQ(ierr); 65 if ( i>= start && i<end) { 66 dx = xx[i-start]; 67 #if !defined(PETSC_COMPLEX) 68 if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-1; 69 else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 70 #else 71 if (abs(dx) < 1.e-16 && real(dx) >= 0.0) dx = 1.e-1; 72 else if (real(dx) < 0.0 && abs(dx) < 1.e-16) dx = -1.e-1; 73 #endif 74 dx *= epsilon; 75 wscale = -1.0/dx; 76 VecSetValues(x2,1,&i,&dx,ADD_VALUES); 77 } 78 else { 79 wscale = 0.0; 80 } 81 ierr = eval_fct(snes,x2,j2); CHKERRQ(ierr); 82 ierr = VecAXPY(&mone,j1,j2); CHKERRQ(ierr); 83 /* communicate scale to all processors */ 84 #if !defined(PETSC_COMPLEX) 85 MPI_Allreduce(&wscale,&scale,1,MPI_DOUBLE,MPI_SUM,comm); 86 #else 87 MPI_Allreduce(&wscale,&scale,2,MPI_DOUBLE,MPI_SUM,comm); 88 #endif 89 scale = -scale; 90 VecScale(&scale,j2); 91 VecGetArray(j2,&y); 92 VecNorm(j2,NORM_INFINITY,&amax); amax *= 1.e-14; 93 for ( j=start; j<end; j++ ) { 94 if (PetscAbsScalar(y[j-start]) > amax) { 95 ierr = MatSetValues(*J,1,&j,1,&i,y+j-start,INSERT_VALUES); CHKERRQ(ierr); 96 } 97 } 98 VecRestoreArray(j2,&y); 99 } 100 MatAssemblyBegin(*J,FINAL_ASSEMBLY); 101 VecDestroy(x2); VecDestroy(j1); VecDestroy(j2); 102 MatAssemblyEnd(*J,FINAL_ASSEMBLY); 103 return 0; 104 } 105 106 /*@C 107 SNESDefaultComputeHessian - Computes the Hessian using finite 108 differences. 109 110 Input Parameters: 111 . x1 - compute Hessian at this point 112 . ctx - application's gradient context, as set with SNESSetGradient() 113 114 Output Parameters: 115 . J - Hessian 116 . B - preconditioner, same as Hessian 117 . flag - matrix flag 118 119 Options Database Key: 120 $ -snes_fd 121 122 Notes: 123 This routine is slow and expensive, and is not currently optimized 124 to take advantage of sparsity in the problem. Although 125 SNESDefaultComputeHessian() is not recommended for general use 126 in large-scale applications, It can be useful in checking the 127 correctness of a user-provided Hessian. 128 129 .keywords: SNES, finite differences, Hessian 130 131 .seealso: SNESSetHessian(), SNESTestHessian() 132 @*/ 133 int SNESDefaultComputeHessian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 134 { 135 return SNESDefaultComputeJacobian(snes,x1,J,B,flag,ctx); 136 } 137