1 /* 2 Context for bounded Newton-Krylov type optimization algorithms 3 */ 4 5 #if !defined(__TAO_BNK_H) 6 #define __TAO_BNK_H 7 #include <petsc/private/taoimpl.h> 8 #include <../src/tao/matrix/lmvmmat.h> 9 10 typedef struct { 11 Mat M; 12 Vec W; 13 14 Vec Xold; 15 Vec Gold; 16 Vec Diag; 17 18 /* Scalar values for the solution and step */ 19 PetscReal fold, f, gnorm, dnorm; 20 21 /* Parameters when updating the perturbation added to the Hessian matrix 22 according to the following scheme: 23 24 pert = sval; 25 26 do until convergence 27 shift Hessian by pert 28 solve Newton system 29 30 if (linear solver failed or did not compute a descent direction) 31 use steepest descent direction and increase perturbation 32 33 if (0 == pert) 34 initialize perturbation 35 pert = min(imax, max(imin, imfac * norm(G))) 36 else 37 increase perturbation 38 pert = min(pmax, max(pgfac * pert, pmgfac * norm(G))) 39 fi 40 else 41 use linear solver direction and decrease perturbation 42 43 pert = min(psfac * pert, pmsfac * norm(G)) 44 if (pert < pmin) 45 pert = 0 46 fi 47 fi 48 49 perform line search 50 function and gradient evaluation 51 check convergence 52 od 53 */ 54 PetscReal sval; /* Starting perturbation value, default zero */ 55 56 PetscReal imin; /* Minimum perturbation added during initialization */ 57 PetscReal imax; /* Maximum perturbation added during initialization */ 58 PetscReal imfac; /* Merit function factor during initialization */ 59 60 PetscReal pert; /* Current perturbation value */ 61 PetscReal pmin; /* Minimim perturbation value */ 62 PetscReal pmax; /* Maximum perturbation value */ 63 PetscReal pgfac; /* Perturbation growth factor */ 64 PetscReal psfac; /* Perturbation shrink factor */ 65 PetscReal pmgfac; /* Merit function growth factor */ 66 PetscReal pmsfac; /* Merit function shrink factor */ 67 68 /* Parameters when updating the trust-region radius based on steplength 69 if step < nu1 (very bad step) 70 radius = omega1 * min(norm(d), radius) 71 elif step < nu2 (bad step) 72 radius = omega2 * min(norm(d), radius) 73 elif step < nu3 (okay step) 74 radius = omega3 * radius; 75 elif step < nu4 (good step) 76 radius = max(omega4 * norm(d), radius) 77 else (very good step) 78 radius = max(omega5 * norm(d), radius) 79 fi 80 */ 81 PetscReal nu1; /* used to compute trust-region radius */ 82 PetscReal nu2; /* used to compute trust-region radius */ 83 PetscReal nu3; /* used to compute trust-region radius */ 84 PetscReal nu4; /* used to compute trust-region radius */ 85 86 PetscReal omega1; /* factor used for trust-region update */ 87 PetscReal omega2; /* factor used for trust-region update */ 88 PetscReal omega3; /* factor used for trust-region update */ 89 PetscReal omega4; /* factor used for trust-region update */ 90 PetscReal omega5; /* factor used for trust-region update */ 91 92 /* Parameters when updating the trust-region radius based on reduction 93 94 kappa = ared / pred 95 if kappa < eta1 (very bad step) 96 radius = alpha1 * min(norm(d), radius) 97 elif kappa < eta2 (bad step) 98 radius = alpha2 * min(norm(d), radius) 99 elif kappa < eta3 (okay step) 100 radius = alpha3 * radius; 101 elif kappa < eta4 (good step) 102 radius = max(alpha4 * norm(d), radius) 103 else (very good step) 104 radius = max(alpha5 * norm(d), radius) 105 fi 106 */ 107 PetscReal eta1; /* used to compute trust-region radius */ 108 PetscReal eta2; /* used to compute trust-region radius */ 109 PetscReal eta3; /* used to compute trust-region radius */ 110 PetscReal eta4; /* used to compute trust-region radius */ 111 112 PetscReal alpha1; /* factor used for trust-region update */ 113 PetscReal alpha2; /* factor used for trust-region update */ 114 PetscReal alpha3; /* factor used for trust-region update */ 115 PetscReal alpha4; /* factor used for trust-region update */ 116 PetscReal alpha5; /* factor used for trust-region update */ 117 118 /* Parameters when updating the trust-region radius based on interpolation 119 120 kappa = ared / pred 121 if kappa >= 1.0 - mu1 (very good step) 122 choose tau in [gamma3, gamma4] 123 radius = max(tau * norm(d), radius) 124 elif kappa >= 1.0 - mu2 (good step) 125 choose tau in [gamma2, gamma3] 126 if (tau >= 1.0) 127 radius = max(tau * norm(d), radius) 128 else 129 radius = tau * min(norm(d), radius) 130 fi 131 else (bad step) 132 choose tau in [gamma1, 1.0] 133 radius = tau * min(norm(d), radius) 134 fi 135 */ 136 PetscReal mu1; /* used for model agreement in interpolation */ 137 PetscReal mu2; /* used for model agreement in interpolation */ 138 139 PetscReal gamma1; /* factor used for interpolation */ 140 PetscReal gamma2; /* factor used for interpolation */ 141 PetscReal gamma3; /* factor used for interpolation */ 142 PetscReal gamma4; /* factor used for interpolation */ 143 144 PetscReal theta; /* factor used for interpolation */ 145 146 /* Parameters when initializing trust-region radius based on interpolation */ 147 PetscReal mu1_i; /* used for model agreement in interpolation */ 148 PetscReal mu2_i; /* used for model agreement in interpolation */ 149 150 PetscReal gamma1_i; /* factor used for interpolation */ 151 PetscReal gamma2_i; /* factor used for interpolation */ 152 PetscReal gamma3_i; /* factor used for interpolation */ 153 PetscReal gamma4_i; /* factor used for interpolation */ 154 155 PetscReal theta_i; /* factor used for interpolation */ 156 157 /* Other parameters */ 158 PetscReal min_radius; /* lower bound on initial radius value */ 159 PetscReal max_radius; /* upper bound on trust region radius */ 160 PetscReal epsilon; /* tolerance used when computing ared/pred */ 161 162 PetscInt newt; /* Newton directions attempted */ 163 PetscInt bfgs; /* BFGS directions attempted */ 164 PetscInt sgrad; /* Scaled gradient directions attempted */ 165 PetscInt grad; /* Gradient directions attempted */ 166 167 PetscInt pc_type; /* Preconditioner for the code */ 168 PetscInt bfgs_scale_type; /* Scaling matrix to used for the bfgs preconditioner */ 169 PetscInt init_type; /* Trust-region initialization method */ 170 PetscInt update_type; /* Trust-region update method */ 171 172 PetscInt ksp_atol; 173 PetscInt ksp_rtol; 174 PetscInt ksp_ctol; 175 PetscInt ksp_negc; 176 PetscInt ksp_dtol; 177 PetscInt ksp_iter; 178 PetscInt ksp_othr; 179 PetscBool is_nash, is_stcg, is_gltr; 180 } TAO_BNK; 181 182 #endif /* if !defined(__TAO_BNK_H) */ 183 184 #define BNK_NEWTON 0 185 #define BNK_BFGS 1 186 #define BNK_SCALED_GRADIENT 2 187 #define BNK_GRADIENT 3 188 189 #define BNK_PC_NONE 0 190 #define BNK_PC_AHESS 1 191 #define BNK_PC_BFGS 2 192 #define BNK_PC_PETSC 3 193 #define BNK_PC_TYPES 4 194 195 #define BFGS_SCALE_AHESS 0 196 #define BFGS_SCALE_PHESS 1 197 #define BFGS_SCALE_BFGS 2 198 #define BFGS_SCALE_TYPES 3 199 200 #define BNK_INIT_CONSTANT 0 201 #define BNK_INIT_DIRECTION 1 202 #define BNK_INIT_INTERPOLATION 2 203 #define BNK_INIT_TYPES 3 204 205 #define BNK_UPDATE_STEP 0 206 #define BNK_UPDATE_REDUCTION 1 207 #define BNK_UPDATE_INTERPOLATION 2 208 #define BNK_UPDATE_TYPES 3 209 210 static const char *BNK_PC[64] = {"none", "ahess", "bfgs", "petsc"}; 211 212 static const char *BFGS_SCALE[64] = {"ahess", "phess", "bfgs"}; 213 214 static const char *BNK_INIT[64] = {"constant", "direction", "interpolation"}; 215 216 static const char *BNK_UPDATE[64] = {"step", "reduction", "interpolation"}; 217 218 PETSC_INTERN PetscErrorCode TaoCreate_BNK(Tao); 219 220 PETSC_INTERN PetscErrorCode MatLMVMSolveShell(PC, Vec, Vec); 221 PETSC_INTERN PetscErrorCode TaoBNKInitialize(Tao); 222 PETSC_INTERN PetscErrorCode TaoBNKComputeStep(Tao, PetscInt*); 223 PETSC_INTERN PetscErrorCode TaoBNKUpdateTrustRadius(Tao, PetscReal, PetscReal, PetscInt, PetscBool*); 224