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