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