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