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 "tao_fd_test" 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 ierr = VecSetValue(X,i,h,ADD_VALUES);CHKERRQ(ierr); 75 ierr = VecAssemblyBegin(X);CHKERRQ(ierr); 76 ierr = VecAssemblyEnd(X);CHKERRQ(ierr); 77 78 ierr = TaoComputeObjective(tao,X,&f2);CHKERRQ(ierr); 79 80 ierr = VecSetValue(X,i,-h,ADD_VALUES); 81 ierr = VecAssemblyBegin(X);CHKERRQ(ierr); 82 ierr = VecAssemblyEnd(X);CHKERRQ(ierr); 83 84 if (i>=low && i<high) { 85 g[i-low]=(f2-f)/h; 86 } 87 } 88 ierr = VecRestoreArray(G,&g);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "TaoDefaultComputeHessian" 94 /*@C 95 TaoDefaultComputeHessian - Computes the Hessian using finite differences. 96 97 Collective on Tao 98 99 Input Parameters: 100 + tao - the Tao context 101 . V - compute Hessian at this point 102 - dummy - not used 103 104 Output Parameters: 105 + H - Hessian matrix (not altered in this routine) 106 . B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 107 - flag - flag indicating whether the matrix sparsity structure has changed 108 109 Options Database Key: 110 + -tao_fd - Activates TaoDefaultComputeHessian() 111 - -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD 112 113 Level: advanced 114 115 Notes: 116 This routine is slow and expensive, and is not currently optimized 117 to take advantage of sparsity in the problem. Although 118 TaoDefaultComputeHessian() is not recommended for general use 119 in large-scale applications, It can be useful in checking the 120 correctness of a user-provided Hessian. 121 122 123 124 .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessianColor(), SNESComputeJacobianDefault(), TaoSetGradientRoutine(), TaoDefaultComputeGradient() 125 126 @*/ 127 PetscErrorCode TaoDefaultComputeHessian(Tao tao,Vec V,Mat *H,Mat *B,MatStructure *flag,void *dummy){ 128 PetscErrorCode ierr; 129 MPI_Comm comm; 130 Vec G; 131 SNES snes; 132 133 PetscFunctionBegin; 134 PetscValidHeaderSpecific(V,VEC_CLASSID,2); 135 ierr = VecDuplicate(V,&G);CHKERRQ(ierr); 136 137 ierr = PetscInfo(tao,"TAO Using finite differences w/o coloring to compute Hessian matrix\n");CHKERRQ(ierr); 138 139 ierr = TaoComputeGradient(tao,V,G);CHKERRQ(ierr); 140 141 ierr = PetscObjectGetComm((PetscObject)(*H),&comm);CHKERRQ(ierr); 142 ierr = SNESCreate(comm,&snes);CHKERRQ(ierr); 143 144 ierr = SNESSetFunction(snes,G,Fsnes,tao);CHKERRQ(ierr); 145 ierr = SNESComputeJacobianDefault(snes,V,H,B,flag,tao);CHKERRQ(ierr); 146 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 147 ierr = VecDestroy(&G);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "TaoDefaultComputeHessianColor" 153 /*@C 154 TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences. 155 156 Collective on Tao 157 158 Input Parameters: 159 + tao - the Tao context 160 . V - compute Hessian at this point 161 - ctx - the PetscColoring object (must be of type MatFDColoring) 162 163 Output Parameters: 164 + H - Hessian matrix (not altered in this routine) 165 . B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 166 - flag - flag indicating whether the matrix sparsity structure has changed 167 168 Level: advanced 169 170 171 .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessian(),SNESComputeJacobianDefaultColor(), TaoSetGradientRoutine() 172 173 @*/ 174 PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat *H,Mat *B,MatStructure *flag,void *ctx) 175 { 176 PetscErrorCode ierr; 177 MatFDColoring coloring = (MatFDColoring)ctx; 178 179 PetscFunctionBegin; 180 PetscValidHeaderSpecific(ctx,MAT_FDCOLORING_CLASSID,6); 181 *flag = SAME_NONZERO_PATTERN; 182 183 ierr=PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n");CHKERRQ(ierr); 184 ierr = MatFDColoringApply(*B,coloring,V,flag,ctx);CHKERRQ(ierr); 185 186 if (*H != *B) { 187 ierr = MatAssemblyBegin(*H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 188 ierr = MatAssemblyEnd(*H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 189 } 190 PetscFunctionReturn(0); 191 } 192 193 194