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