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 static PetscErrorCode Fsnes(SNES snes ,Vec X,Vec G,void*ctx) 10 { 11 PetscErrorCode ierr; 12 Tao tao = (Tao)ctx; 13 14 PetscFunctionBegin; 15 PetscValidHeaderSpecific(ctx,TAO_CLASSID,4); 16 ierr=TaoComputeGradient(tao,X,G);CHKERRQ(ierr); 17 PetscFunctionReturn(0); 18 } 19 20 /*@C 21 TaoDefaultComputeGradient - computes the gradient using finite differences. 22 23 Collective on Tao 24 25 Input Parameters: 26 + tao - the Tao context 27 . X - compute gradient at this point 28 - dummy - not used 29 30 Output Parameters: 31 . G - Gradient Vector 32 33 Options Database Key: 34 + -tao_fd_gradient - Activates TaoDefaultComputeGradient() 35 - -tao_fd_delta <delta> - change in x used to calculate finite differences 36 37 Level: advanced 38 39 Note: 40 This routine is slow and expensive, and is not currently optimized 41 to take advantage of sparsity in the problem. Although 42 TaoAppDefaultComputeGradient is not recommended for general use 43 in large-scale applications, It can be useful in checking the 44 correctness of a user-provided gradient. Use the tao method TAOTEST 45 to get an indication of whether your gradient is correct. 46 47 48 Note: 49 This finite difference gradient evaluation can be set using the routine TaoSetGradientRoutine() or by using the command line option -tao_fd_gradient 50 51 .seealso: TaoSetGradientRoutine() 52 53 @*/ 54 PetscErrorCode TaoDefaultComputeGradient(Tao tao,Vec X,Vec G,void *dummy) 55 { 56 PetscScalar *x,*g; 57 PetscReal f, f2; 58 PetscErrorCode ierr; 59 PetscInt low,high,N,i; 60 PetscBool flg; 61 PetscReal h=.5*PETSC_SQRT_MACHINE_EPSILON; 62 63 PetscFunctionBegin; 64 ierr = PetscOptionsGetReal(((PetscObject)tao)->options,((PetscObject)tao)->prefix,"-tao_fd_delta",&h,&flg);CHKERRQ(ierr); 65 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 66 ierr = VecGetOwnershipRange(X,&low,&high);CHKERRQ(ierr); 67 ierr = VecGetArray(G,&g);CHKERRQ(ierr); 68 for (i=0;i<N;i++) { 69 if (i>=low && i<high) { 70 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 71 x[i-low] -= h; 72 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 73 } 74 75 ierr = TaoComputeObjective(tao, X,&f);CHKERRQ(ierr); 76 77 if (i>=low && i<high) { 78 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 79 x[i-low] += 2*h; 80 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 81 } 82 83 ierr = TaoComputeObjective(tao,X,&f2);CHKERRQ(ierr); 84 85 if (i>=low && i<high) { 86 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 87 x[i-low] -= h; 88 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 89 } 90 if (i>=low && i<high) { 91 g[i-low]=(f2-f)/(2.0*h); 92 } 93 } 94 ierr = VecRestoreArray(G,&g);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 98 /*@C 99 TaoDefaultComputeHessian - Computes the Hessian using finite differences. 100 101 Collective on Tao 102 103 Input Parameters: 104 + tao - the Tao context 105 . V - compute Hessian at this point 106 - dummy - not used 107 108 Output Parameters: 109 + H - Hessian matrix (not altered in this routine) 110 - B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 111 112 Options Database Key: 113 + -tao_fd - Activates TaoDefaultComputeHessian() 114 - -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD 115 116 Level: advanced 117 118 Notes: 119 This routine is slow and expensive, and is not currently optimized 120 to take advantage of sparsity in the problem. Although 121 TaoDefaultComputeHessian() is not recommended for general use 122 in large-scale applications, It can be useful in checking the 123 correctness of a user-provided Hessian. 124 125 .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessianColor(), SNESComputeJacobianDefault(), TaoSetGradientRoutine(), TaoDefaultComputeGradient() 126 127 @*/ 128 PetscErrorCode TaoDefaultComputeHessian(Tao tao,Vec V,Mat H,Mat B,void *dummy) 129 { 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,tao);CHKERRQ(ierr); 148 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 149 ierr = VecDestroy(&G);CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 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 167 Level: advanced 168 169 170 .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessian(),SNESComputeJacobianDefaultColor(), TaoSetGradientRoutine() 171 172 @*/ 173 PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H,Mat B,void *ctx) 174 { 175 PetscErrorCode ierr; 176 MatFDColoring coloring = (MatFDColoring)ctx; 177 178 PetscFunctionBegin; 179 PetscValidHeaderSpecific(ctx,MAT_FDCOLORING_CLASSID,6); 180 ierr=PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n");CHKERRQ(ierr); 181 ierr = MatFDColoringApply(B,coloring,V,ctx);CHKERRQ(ierr); 182 if (H != B) { 183 ierr = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 184 ierr = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 185 } 186 PetscFunctionReturn(0); 187 } 188 189 190