1 // Copyright (C) 2005, 2009 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Eclipse Public License. 4 // 5 // $Id$ 6 // 7 // Authors: Andreas Waechter IBM 2009-11-05 8 // (based on IpLowRankAugSystemSolver.hpp rev 1324) 9 10 #ifndef __IP_LOWRANKSSAUGSYSTEMSOLVER_HPP__ 11 #define __IP_LOWRANKSSAUGSYSTEMSOLVER_HPP__ 12 13 #include "IpAugSystemSolver.hpp" 14 #include "IpDiagMatrix.hpp" 15 #include "IpCompoundMatrix.hpp" 16 #include "IpCompoundVector.hpp" 17 #include "IpExpandedMultiVectorMatrix.hpp" 18 19 namespace Ipopt 20 { 21 22 /** Solver for the augmented system with LowRankUpdateSymMatrix 23 * Hessian matrices. This version works with only one backsolve 24 * (so it is better for iterative linear solvers), by augmenting 25 * the regular augmented system. 26 */ 27 class LowRankSSAugSystemSolver : public AugSystemSolver 28 { 29 public: 30 /**@name Constructors/Destructors */ 31 //@{ 32 /** Constructor using an existing augmented system solver. the 33 * max_rank argument is the maximal rank that can appear. */ 34 LowRankSSAugSystemSolver(AugSystemSolver& aug_system_solver, 35 Index max_rank); 36 37 /** Default destructor */ 38 virtual ~LowRankSSAugSystemSolver(); 39 //@} 40 41 /** overloaded from AlgorithmStrategyObject */ 42 bool InitializeImpl(const OptionsList& options, 43 const std::string& prefix); 44 45 /** Set up the augmented system and solve it for a given right hand 46 * side. 47 */ 48 virtual ESymSolverStatus Solve( 49 const SymMatrix* W, 50 double W_factor, 51 const Vector* D_x, 52 double delta_x, 53 const Vector* D_s, 54 double delta_s, 55 const Matrix* J_c, 56 const Vector* D_c, 57 double delta_c, 58 const Matrix* J_d, 59 const Vector* D_d, 60 double delta_d, 61 const Vector& rhs_x, 62 const Vector& rhs_s, 63 const Vector& rhs_c, 64 const Vector& rhs_d, 65 Vector& sol_x, 66 Vector& sol_s, 67 Vector& sol_c, 68 Vector& sol_d, 69 bool check_NegEVals, 70 Index numberOfNegEVals); 71 72 /** Number of negative eigenvalues detected during last 73 * solve. Returns the number of negative eigenvalues of 74 * the most recent factorized matrix. This must not be called if 75 * the linear solver does not compute this quantities (see 76 * ProvidesInertia). 77 */ 78 virtual Index NumberOfNegEVals() const; 79 80 /** Query whether inertia is computed by linear solver. 81 * Returns true, if linear solver provides inertia. 82 */ 83 virtual bool ProvidesInertia() const; 84 85 /** Request to increase quality of solution for next solve. Ask 86 * underlying linear solver to increase quality of solution for 87 * the next solve (e.g. increase pivot tolerance). Returns 88 * false, if this is not possible (e.g. maximal pivot tolerance 89 * already used.) 90 */ 91 virtual bool IncreaseQuality(); 92 93 private: 94 /**@name Default Compiler Generated Methods 95 * (Hidden to avoid implicit creation/calling). 96 * These methods are not implemented and 97 * we do not want the compiler to implement 98 * them for us, so we declare them private 99 * and do not define them. This ensures that 100 * they will not be implicitly created/called. */ 101 //@{ 102 /** Default constructor. */ 103 LowRankSSAugSystemSolver(); 104 /** Copy Constructor */ 105 LowRankSSAugSystemSolver(const LowRankSSAugSystemSolver&); 106 107 /** Overloaded Equals Operator */ 108 void operator=(const LowRankSSAugSystemSolver&); 109 //@} 110 111 /** The augmented system solver object that should be used for the 112 * factorization of the augmented system without the low-rank 113 * update. 114 */ 115 SmartPtr<AugSystemSolver> aug_system_solver_; 116 117 /** Maximal rank of low rank Hessian update */ 118 Index max_rank_; 119 120 /**@name Tags and values to track in order to decide whether the 121 matrix has to be updated compared to the most recent call of 122 the Set method. 123 */ 124 //@{ 125 /** Tag for W matrix. If W has been given to Set as NULL, then 126 * this tag is set to 0 127 */ 128 TaggedObject::Tag w_tag_; 129 /** Most recent value of W_factor */ 130 double w_factor_; 131 /** Tag for D_x vector, representing the diagonal matrix D_x. If 132 * D_x has been given to Set as NULL, then this tag is set to 0 133 */ 134 TaggedObject::Tag d_x_tag_; 135 /** Most recent value of delta_x from Set method */ 136 double delta_x_; 137 /** Tag for D_s vector, representing the diagonal matrix D_s. If 138 * D_s has been given to Set as NULL, then this tag is set to 0 139 */ 140 TaggedObject::Tag d_s_tag_; 141 /** Most recent value of delta_s from Set method */ 142 double delta_s_; 143 /** Tag for J_c matrix. If J_c has been given to Set as NULL, then 144 * this tag is set to 0 145 */ 146 TaggedObject::Tag j_c_tag_; 147 /** Tag for D_c vector, representing the diagonal matrix D_c. If 148 * D_c has been given to Set as NULL, then this tag is set to 0 149 */ 150 TaggedObject::Tag d_c_tag_; 151 /** Most recent value of delta_c from Set method */ 152 double delta_c_; 153 /** Tag for J_d matrix. If J_d has been given to Set as NULL, then 154 * this tag is set to 0 155 */ 156 TaggedObject::Tag j_d_tag_; 157 /** Tag for D_d vector, representing the diagonal matrix D_d. If 158 * D_d has been given to Set as NULL, then this tag is set to 0 159 */ 160 TaggedObject::Tag d_d_tag_; 161 /** Most recent value of delta_d from Set method */ 162 double delta_d_; 163 //@} 164 165 /** Flag indicating if this is the first call */ 166 bool first_call_; 167 168 /** @name Information to be stored in order to resolve for the 169 * same matrix with a different right hand side. */ 170 //@{ 171 /** Hessian Matrix passed to the augmented system solver solving 172 * the matrix without the low-rank update. */ 173 SmartPtr<DiagMatrix> Wdiag_; 174 /** Artifical rows for Jac_c part for low rank data */ 175 SmartPtr<ExpandedMultiVectorMatrix> expanded_vu_; 176 /** Extended Jac_c to include expanded_vu_ */ 177 SmartPtr<CompoundMatrix> J_c_ext_; 178 /** Extended D_c diagonal */ 179 SmartPtr<CompoundVector> D_c_ext_; 180 /** Extended vector space for y_c */ 181 SmartPtr<CompoundVectorSpace> y_c_ext_space_; 182 /** Number of components in V, so that it can be used to correct 183 * the inertia */ 184 Index negEvalsCorrection_; 185 //@} 186 187 /** Stores the number of negative eigenvalues detected during most 188 * recent factorization. This is what is returned by 189 * NumberOfNegEVals() of this class. It usually is the number of 190 * negative eigenvalues retured from the aug_system_solver solve, 191 * but if a Cholesky factorization could not be performed, the 192 * returned value is one more than this what the 193 * aug_system_solver returned. */ 194 Index num_neg_evals_; 195 196 /** @name Internal functions */ 197 //@{ 198 /** Method for updating the factorization, including J1_, J2_, 199 * Vtilde1_, Utilde2, Wdiag_, compound_sol_vecspace_ */ 200 ESymSolverStatus UpdateExtendedData( 201 const SymMatrix* W, 202 double W_factor, 203 const Vector* D_x, 204 double delta_x, 205 const Vector* D_s, 206 double delta_s, 207 const Matrix& J_c, 208 const Vector* D_c, 209 double delta_c, 210 const Matrix& J_d, 211 const Vector* D_d, 212 double delta_d, 213 const Vector& proto_rhs_x, 214 const Vector& proto_rhs_s, 215 const Vector& proto_rhs_c, 216 const Vector& proto_rhs_d); 217 218 /** Method that compares the tags of the data for the matrix with 219 * those from the previous call. Returns true, if there was a 220 * change and the factorization has to be updated. */ 221 bool AugmentedSystemRequiresChange( 222 const SymMatrix* W, 223 double W_factor, 224 const Vector* D_x, 225 double delta_x, 226 const Vector* D_s, 227 double delta_s, 228 const Matrix& J_c, 229 const Vector* D_c, 230 double delta_c, 231 const Matrix& J_d, 232 const Vector* D_d, 233 double delta_d); 234 //@} 235 236 }; 237 238 } // namespace Ipopt 239 240 #endif 241