1 #include <petsctao.h> /*I "petsctao.h" I*/ 2 #include <petsc-private/taoimpl.h> 3 #include <petscsnes.h> 4 5 /* 6 For finited difference computations of the Hessian, we use PETSc's SNESComputeJacobianDefault 7 */ 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "Fsnes" 11 static PetscErrorCode Fsnes(SNES snes ,Vec X,Vec G,void*ctx) 12 { 13 PetscErrorCode ierr; 14 Tao tao = (Tao)ctx; 15 16 PetscFunctionBegin; 17 PetscValidHeaderSpecific(ctx,TAO_CLASSID,4); 18 ierr=TaoComputeGradient(tao,X,G);CHKERRQ(ierr); 19 PetscFunctionReturn(0); 20 } 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "TaoDefaultComputeGradient" 24 /*@C 25 TaoDefaultComputeGradient - computes the gradient using finite differences. 26 27 Collective on Tao 28 29 Input Parameters: 30 + tao - the Tao context 31 . X - compute gradient at this point 32 - dummy - not used 33 34 Output Parameters: 35 . G - Gradient Vector 36 37 Options Database Key: 38 + -tao_fd_gradient - Activates TaoDefaultComputeGradient() 39 - -tao_fd_delta <delta> - change in x used to calculate finite differences 40 41 Level: advanced 42 43 Note: 44 This routine is slow and expensive, and is not currently optimized 45 to take advantage of sparsity in the problem. Although 46 TaoAppDefaultComputeGradient is not recommended for general use 47 in large-scale applications, It can be useful in checking the 48 correctness of a user-provided gradient. Use the tao method TAOTEST 49 to get an indication of whether your gradient is correct. 50 51 52 Note: 53 This finite difference gradient evaluation can be set using the routine TaoSetGradientRoutine() or by using the command line option -tao_fd_gradient 54 55 .seealso: TaoSetGradientRoutine() 56 57 @*/ 58 PetscErrorCode TaoDefaultComputeGradient(Tao tao,Vec X,Vec G,void *dummy) 59 { 60 PetscReal *g; 61 PetscReal f, f2; 62 PetscErrorCode ierr; 63 PetscInt low,high,N,i; 64 PetscBool flg; 65 PetscReal h=PETSC_SQRT_MACHINE_EPSILON; 66 67 PetscFunctionBegin; 68 ierr = TaoComputeObjective(tao, X,&f);CHKERRQ(ierr); 69 ierr = PetscOptionsGetReal(NULL,"-tao_fd_delta",&h,&flg);CHKERRQ(ierr); 70 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 71 ierr = VecGetOwnershipRange(X,&low,&high);CHKERRQ(ierr); 72 ierr = VecGetArray(G,&g);CHKERRQ(ierr); 73 for (i=0;i<N;i++) { 74 if (i>=low && i<high) { 75 PetscScalar *xx; 76 ierr = VecGetArray(X,&xx);CHKERRQ(ierr); 77 xx[i-low] += h; 78 ierr = VecRestoreArray(X,&xx);CHKERRQ(ierr); 79 } 80 81 ierr = TaoComputeObjective(tao,X,&f2);CHKERRQ(ierr); 82 83 if (i>=low && i<high) { 84 PetscScalar *xx; 85 ierr = VecGetArray(X,&xx);CHKERRQ(ierr); 86 xx[i-low] -= h; 87 ierr = VecRestoreArray(X,&xx);CHKERRQ(ierr); 88 } 89 if (i>=low && i<high) { 90 g[i-low]=(f2-f)/h; 91 } 92 } 93 ierr = VecRestoreArray(G,&g);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "TaoDefaultComputeHessian" 99 /*@C 100 TaoDefaultComputeHessian - Computes the Hessian using finite differences. 101 102 Collective on Tao 103 104 Input Parameters: 105 + tao - the Tao context 106 . V - compute Hessian at this point 107 - dummy - not used 108 109 Output Parameters: 110 + H - Hessian matrix (not altered in this routine) 111 . B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 112 - flag - flag indicating whether the matrix sparsity structure has changed 113 114 Options Database Key: 115 + -tao_fd - Activates TaoDefaultComputeHessian() 116 - -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD 117 118 Level: advanced 119 120 Notes: 121 This routine is slow and expensive, and is not currently optimized 122 to take advantage of sparsity in the problem. Although 123 TaoDefaultComputeHessian() is not recommended for general use 124 in large-scale applications, It can be useful in checking the 125 correctness of a user-provided Hessian. 126 127 128 129 .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessianColor(), SNESComputeJacobianDefault(), TaoSetGradientRoutine(), TaoDefaultComputeGradient() 130 131 @*/ 132 PetscErrorCode TaoDefaultComputeHessian(Tao tao,Vec V,Mat H,Mat B,MatStructure *flag,void *dummy) 133 { 134 PetscErrorCode ierr; 135 MPI_Comm comm; 136 Vec G; 137 SNES snes; 138 139 PetscFunctionBegin; 140 PetscValidHeaderSpecific(V,VEC_CLASSID,2); 141 ierr = VecDuplicate(V,&G);CHKERRQ(ierr); 142 143 ierr = PetscInfo(tao,"TAO Using finite differences w/o coloring to compute Hessian matrix\n");CHKERRQ(ierr); 144 145 ierr = TaoComputeGradient(tao,V,G);CHKERRQ(ierr); 146 147 ierr = PetscObjectGetComm((PetscObject)H,&comm);CHKERRQ(ierr); 148 ierr = SNESCreate(comm,&snes);CHKERRQ(ierr); 149 150 ierr = SNESSetFunction(snes,G,Fsnes,tao);CHKERRQ(ierr); 151 ierr = SNESComputeJacobianDefault(snes,V,H,B,tao);CHKERRQ(ierr); 152 *flag = SAME_NONZERO_PATTERN; 153 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 154 ierr = VecDestroy(&G);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "TaoDefaultComputeHessianColor" 160 /*@C 161 TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences. 162 163 Collective on Tao 164 165 Input Parameters: 166 + tao - the Tao context 167 . V - compute Hessian at this point 168 - ctx - the PetscColoring object (must be of type MatFDColoring) 169 170 Output Parameters: 171 + H - Hessian matrix (not altered in this routine) 172 . B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 173 - flag - flag indicating whether the matrix sparsity structure has changed 174 175 Level: advanced 176 177 178 .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessian(),SNESComputeJacobianDefaultColor(), TaoSetGradientRoutine() 179 180 @*/ 181 PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H,Mat B,MatStructure *flag,void *ctx) 182 { 183 PetscErrorCode ierr; 184 MatFDColoring coloring = (MatFDColoring)ctx; 185 186 PetscFunctionBegin; 187 PetscValidHeaderSpecific(ctx,MAT_FDCOLORING_CLASSID,6); 188 *flag = SAME_NONZERO_PATTERN; 189 190 ierr=PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n");CHKERRQ(ierr); 191 ierr = MatFDColoringApply(B,coloring,V,ctx);CHKERRQ(ierr); 192 *flag = SAME_NONZERO_PATTERN; 193 194 if (H != B) { 195 ierr = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 196 ierr = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 197 } 198 PetscFunctionReturn(0); 199 } 200 201 202