1 // Copyright (C) 2004, 2006 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Common Public License. 4 // 5 // $Id: IpPDFullSpaceSolver.hpp 759 2006-07-07 03:07:08Z andreasw $ 6 // 7 // Authors: Carl Laird, Andreas Waechter IBM 2004-08-13 8 9 #ifndef __IPPDFULLSPACESOLVER_HPP__ 10 #define __IPPDFULLSPACESOLVER_HPP__ 11 12 #include "IpPDSystemSolver.hpp" 13 #include "IpAugSystemSolver.hpp" 14 #include "IpPDPerturbationHandler.hpp" 15 16 namespace SimTKIpopt 17 { 18 19 /** This is the implemetation of the Primal-Dual System, using the 20 * full space approach with a direct linear solver. 21 * 22 * A note on the iterative refinement: We perform at least 23 * min_refinement_steps number of iterative refinement steps. If after 24 * one iterative refinement the quality of the solution (defined in 25 * ResidualRatio) does not improve or the maximal number of 26 * iterative refinement steps is exceeded before the tolerance 27 * residual_ratio_max_ is satisfied, we first ask the linear solver 28 * to solve the system more accurately (e.g. by increasing the 29 * pivot tolerance). If that doesn't help or is not possible, we 30 * treat the system, as if it is singular (i.e. increase delta's). 31 */ 32 class PDFullSpaceSolver: public PDSystemSolver 33 { 34 public: 35 /** @name /Destructor */ 36 //@{ 37 /** Constructor that takes in the Augmented System solver that 38 * is to be used inside 39 */ 40 PDFullSpaceSolver(AugSystemSolver& augSysSolver, 41 PDPerturbationHandler& perturbHandler); 42 43 /** Default destructor */ 44 virtual ~PDFullSpaceSolver(); 45 //@} 46 47 /* overloaded from AlgorithmStrategyObject */ 48 bool InitializeImpl(const OptionsList& options, 49 const std::string& prefix) override; 50 51 /** Solve the primal dual system, given one right hand side. 52 */ 53 virtual bool Solve(Number alpha, 54 Number beta, 55 const IteratesVector& rhs, 56 IteratesVector& res, 57 bool allow_inexact=false, 58 bool improve_solution=false) override; 59 60 /** Methods for IpoptType */ 61 //@{ 62 static void RegisterOptions(SmartPtr<RegisteredOptions> roptions); 63 //@} 64 65 private: 66 /**@name Default Compiler Generated Methods 67 * (Hidden to avoid implicit creation/calling). 68 * These methods are not implemented and 69 * we do not want the compiler to implement 70 * them for us, so we declare them private 71 * and do not define them. This ensures that 72 * they will not be implicitly created/called. */ 73 //@{ 74 /** Default Constructor */ 75 PDFullSpaceSolver(); 76 /** Overloaded Equals Operator */ 77 PDFullSpaceSolver& operator=(const PDFullSpaceSolver&); 78 //@} 79 80 /** @name Strategy objects to hold on to. */ 81 //@{ 82 /** Pointer to the Solver for the augmented system */ 83 SmartPtr<AugSystemSolver> augSysSolver_; 84 /** Pointer to the Perturbation Handler. */ 85 SmartPtr<PDPerturbationHandler> perturbHandler_; 86 //@} 87 88 /**@name Data about the correction made to the system */ 89 //@{ 90 /** A dummy cache to figure out if the deltas are still up to date*/ 91 CachedResults<void*> dummy_cache_; 92 /** Flag indicating if for the current matrix the solution quality 93 * of the augmented system solver has already been increased. */ 94 bool augsys_improved_; 95 //@} 96 97 /** @name Parameters */ 98 //@{ 99 /** Minimal number of iterative refinement performed per backsolve */ 100 Index min_refinement_steps_; 101 /** Maximal number of iterative refinement performed per backsolve */ 102 Index max_refinement_steps_; 103 /** Maximal allowed ratio of the norm of the residual over the 104 * norm of the right hand side and solution. */ 105 Number residual_ratio_max_; 106 /** If the residual_ratio is larger than this value after trying 107 * to improve the solution, the linear system is assumed to be 108 * singular and modified. */ 109 Number residual_ratio_singular_; 110 /** Factor defining require improvement to consider iterative 111 * refinement successful. */ 112 Number residual_improvement_factor_; 113 //@} 114 115 /** Internal function for a single backsolve (which will be used 116 * for iterative refinement on the outside). This method returns 117 * false, if for some reason the linear system could not be 118 * solved (e.g. when the regularization parameter becomes too 119 * large.) 120 */ 121 bool SolveOnce(bool resolve_unmodified, 122 bool pretend_singular, 123 const SymMatrix& W, 124 const Matrix& J_c, 125 const Matrix& J_d, 126 const Matrix& Px_L, 127 const Matrix& Px_U, 128 const Matrix& Pd_L, 129 const Matrix& Pd_U, 130 const Vector& z_L, 131 const Vector& z_U, 132 const Vector& v_L, 133 const Vector& v_U, 134 const Vector& slack_x_L, 135 const Vector& slack_x_U, 136 const Vector& slack_s_L, 137 const Vector& slack_s_U, 138 const Vector& sigma_x, 139 const Vector& sigma_s, 140 Number alpha, 141 Number beta, 142 const IteratesVector& rhs, 143 IteratesVector& res); 144 145 /** Internal function for computing the residual (resid) given the 146 * right hand side (rhs) and the solution of the system (res). 147 */ 148 void ComputeResiduals(const SymMatrix& W, 149 const Matrix& J_c, 150 const Matrix& J_d, 151 const Matrix& Px_L, 152 const Matrix& Px_U, 153 const Matrix& Pd_L, 154 const Matrix& Pd_U, 155 const Vector& z_L, 156 const Vector& z_U, 157 const Vector& v_L, 158 const Vector& v_U, 159 const Vector& slack_x_L, 160 const Vector& slack_x_U, 161 const Vector& slack_s_L, 162 const Vector& slack_s_U, 163 const Vector& sigma_x, 164 const Vector& sigma_s, 165 Number alpha, 166 Number beta, 167 const IteratesVector& rhs, 168 const IteratesVector& res, 169 IteratesVector& resid); 170 171 /** Internal function for computing the ratio of the residual 172 * compared to the right hand side and solution. The smaller 173 * this value, the better the solution. */ 174 Number ComputeResidualRatio(const IteratesVector& rhs, 175 const IteratesVector& res, 176 const IteratesVector& resid); 177 178 /** @name Auxilliary functions */ 179 //@{ 180 /** Compute \f$ x = S^{-1}(r + \alpha Z P^T d)\f$ */ 181 void SinvBlrmZPTdBr(Number alpha, const Vector& S, 182 const Vector& R, const Vector& Z, 183 const Matrix& P, const Vector&g, Vector& X); 184 //@} 185 }; 186 187 } // namespace Ipopt 188 189 #endif 190