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