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