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