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