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(PETSC_SUCCESS); 17 } 18 19 /*@C 20 TaoDefaultComputeGradient - computes the gradient using finite differences. 21 22 Collective 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) g[i - low] = (f2 - f) / (2.0 * h); 79 } 80 PetscCall(VecRestoreArray(G, &g)); 81 PetscCall(VecDestroy(&X)); 82 PetscFunctionReturn(PETSC_SUCCESS); 83 } 84 85 /*@C 86 TaoDefaultComputeHessian - Computes the Hessian using finite differences. 87 88 Collective 89 90 Input Parameters: 91 + tao - the Tao context 92 . V - compute Hessian at this point 93 - dummy - not used 94 95 Output Parameters: 96 + H - Hessian matrix (not altered in this routine) 97 - B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 98 99 Options Database Key: 100 . -tao_fd_hessian - activates TaoDefaultComputeHessian() 101 102 Level: advanced 103 104 Notes: 105 This routine is slow and expensive, and is not optimized 106 to take advantage of sparsity in the problem. Although 107 it is not recommended for general use 108 in large-scale applications, It can be useful in checking the 109 correctness of a user-provided Hessian. 110 111 .seealso: `Tao`, `TaoSetHessian()`, `TaoDefaultComputeHessianColor()`, `SNESComputeJacobianDefault()`, `TaoSetGradient()`, `TaoDefaultComputeGradient()` 112 @*/ 113 PetscErrorCode TaoDefaultComputeHessian(Tao tao, Vec V, Mat H, Mat B, void *dummy) 114 { 115 SNES snes; 116 DM dm; 117 118 PetscFunctionBegin; 119 PetscCall(PetscInfo(tao, "TAO Using finite differences w/o coloring to compute Hessian matrix\n")); 120 PetscCall(SNESCreate(PetscObjectComm((PetscObject)H), &snes)); 121 PetscCall(SNESSetFunction(snes, NULL, Fsnes, tao)); 122 PetscCall(SNESGetDM(snes, &dm)); 123 PetscCall(DMShellSetGlobalVector(dm, V)); 124 PetscCall(SNESSetUp(snes)); 125 if (H) { 126 PetscInt n, N; 127 128 PetscCall(VecGetSize(V, &N)); 129 PetscCall(VecGetLocalSize(V, &n)); 130 PetscCall(MatSetSizes(H, n, n, N, N)); 131 PetscCall(MatSetUp(H)); 132 } 133 if (B && B != H) { 134 PetscInt n, N; 135 136 PetscCall(VecGetSize(V, &N)); 137 PetscCall(VecGetLocalSize(V, &n)); 138 PetscCall(MatSetSizes(B, n, n, N, N)); 139 PetscCall(MatSetUp(B)); 140 } 141 PetscCall(SNESComputeJacobianDefault(snes, V, H, B, NULL)); 142 PetscCall(SNESDestroy(&snes)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 /*@C 147 TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences. 148 149 Collective 150 151 Input Parameters: 152 + tao - the Tao context 153 . V - compute Hessian at this point 154 - ctx - the color object of type `MatFDColoring` 155 156 Output Parameters: 157 + H - Hessian matrix (not altered in this routine) 158 - B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 159 160 Level: advanced 161 162 .seealso: `Tao`, `MatColoring`, `TaoSetHessian()`, `TaoDefaultComputeHessian()`, `SNESComputeJacobianDefaultColor()`, `TaoSetGradient()` 163 @*/ 164 PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H, Mat B, void *ctx) 165 { 166 MatFDColoring coloring = (MatFDColoring)ctx; 167 168 PetscFunctionBegin; 169 PetscValidHeaderSpecific(coloring, MAT_FDCOLORING_CLASSID, 5); 170 PetscCall(PetscInfo(tao, "TAO computing matrix using finite differences Hessian and coloring\n")); 171 PetscCall(MatFDColoringApply(B, coloring, V, ctx)); 172 if (H != B) { 173 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 174 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 175 } 176 PetscFunctionReturn(PETSC_SUCCESS); 177 } 178 179 PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao, Vec X, Mat H, Mat B, void *ctx) 180 { 181 PetscInt n, N; 182 PetscBool assembled; 183 184 PetscFunctionBegin; 185 PetscCheck(!B || B == H, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "Preconditioning Hessian matrix"); 186 PetscCall(MatAssembled(H, &assembled)); 187 if (!assembled) { 188 PetscCall(VecGetSize(X, &N)); 189 PetscCall(VecGetLocalSize(X, &n)); 190 PetscCall(MatSetSizes(H, n, n, N, N)); 191 PetscCall(MatSetType(H, MATMFFD)); 192 PetscCall(MatSetUp(H)); 193 PetscCall(MatMFFDSetFunction(H, (PetscErrorCode(*)(void *, Vec, Vec))TaoComputeGradient, tao)); 194 } 195 PetscCall(MatMFFDSetBase(H, X, NULL)); 196 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 197 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 198 PetscFunctionReturn(PETSC_SUCCESS); 199 } 200