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