source: branches/dev/Algorithm/IpPDFullSpaceSolver.hpp @ 529

Last change on this file since 529 was 529, checked in by andreasw, 15 years ago
  • in Vector class, copy cached values in Copy and update them in Scal
  • perform iterative refinement only once for adaptive strategy
  • fix bugs in PDPerturbationHandler
  • avoid some overhead in CalculateQualityFunction?
  • minor changes in timing
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.0 KB
Line 
1// Copyright (C) 2004, 2005 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpPDFullSpaceSolver.hpp 529 2005-09-29 21:12:38Z 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
16namespace Ipopt
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);
50
51    /** Solve the primal dual system, given one right hand side.
52     */
53    virtual void Solve(Number alpha,
54                       Number beta,
55                       const IteratesVector& rhs,
56                       IteratesVector& res,
57                       bool allow_inexact=false,
58                       bool improve_solution=false);
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
Note: See TracBrowser for help on using the repository browser.