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