// Copyright (C) 2004, 2006 International Business Machines and others.
// All Rights Reserved.
// This code is published under the Eclipse Public License.
//
// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13

#ifndef __IPPDSYSTEMSOLVER_HPP__
#define __IPPDSYSTEMSOLVER_HPP__

#include "IpUtils.hpp"
#include "IpSymMatrix.hpp"
#include "IpAlgStrategy.hpp"
#include "IpIteratesVector.hpp"

namespace Ipopt
{

/** Pure Primal Dual System Solver Base Class.
 *  This is the base class for all derived Primal-Dual System Solver Types.
 *
 *  Here, we understand the primal-dual system as the following linear
 *  system:
 *
 *  \f$
 *  \left[\begin{array}{cccccccc}
 *  W & 0 & J_c^T & J_d^T & -P^x_L & P^x_U & 0 & 0 \\
 *  0 & 0 & 0 & -I & 0 & 0 & -P_L^d & P_U^d \\
 *  J_c & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
 *  J_d & -I & 0 & 0 & 0 & 0 & 0 & 0\\
 *  Z_L(P_L^x)^T & 0 & 0 & 0 & Sl^x_L & 0 & 0 & 0\\
 *  -Z_U(P_U^x)^T & 0 & 0 & 0 & 0 & Sl^x_U & 0 & 0\\
 *  0 & V_L(P_L^d)^T & 0 & 0 & 0 & 0 & Sl^s_L & 0 \\
 *  0 & -V_U(P_U^d)^T & 0 & 0 & 0 & 0 & 0 & Sl^s_U \\
 *  \end{array}\right]
 *  \left(\begin{array}{c}
 *  sol_x\\ sol_s\\ sol_c\\ sol_d\\ sol^z_L\\ sol^z_U\\ sol^v_L\\
 *  sol^v_U
 *  \end{array}\right) =
 *  \left(\begin{array}{c}
 *  rhs_x\\ rhs_s\\ rhs_c\\ rhs_d\\ rhs^z_L\\ rhs^z_U\\ rhs^v_L\\
 *  rhs^v_U
 *  \end{array}\right)
 *  \f$
 *
 *  Here, \f$Sl^x_L = (P^x_L)^T x - x_L\f$,
 *  \f$Sl^x_U = x_U - (P^x_U)^T x\f$, \f$Sl^d_L = (P^d_L)^T d(x) - d_L\f$,
 *  \f$Sl^d_U = d_U - (P^d_U)^T d(x)\f$.  The results returned to the
 *  caller is \f$res = \alpha * sol + \beta * res\f$.
 *
 *  The solution of this linear system (in order to compute the search
 *  direction of the algorithm) usually requires a considerable amount of
 *  computation time.  Therefore, it is important to tailor the solution
 *  of this system to the characteristics of the problem.  The purpose of
 *  this base class is to provide a generic interface to the algorithm
 *  that it can use whenever it requires a solution of the above system.
 *  Particular implementation can then be written to provide the methods
 *  defined here.
 *
 *  It is implicitly assumed here, that the upper left 2 by 2 block
 *  is possibly modified (implicitly or explicitly) so that its
 *  projection onto the null space of the overall constraint
 *  Jacobian \f$\left[\begin{array}{cc}J_c & 0\\J_d &
 *  -I\end{array}\right]\f$ is positive definite.  This is necessary
 *  to guarantee certain descent properties of the resulting search
 *  direction.  For example, in the full space implementation, a
 *  multiple of the identity might be added to the upper left 2 by 2
 *  block.
 *
 *  Note that the Solve method might be called several times for different
 *  right hand sides, but with identical data.  Therefore, if possible,
 *  an implementation of PDSystem should check whether the incoming data has
 *  changed, and not redo factorization etc. unless necessary.
 */
class IPOPTLIB_EXPORT PDSystemSolver: public AlgorithmStrategyObject
{
public:
   /** @name Constructor/Destructor */
   ///@{
   /** Default Constructor */
   PDSystemSolver()
   { }

   /** Destructor */
   virtual ~PDSystemSolver()
   { }
   ///@}

   /** overloaded from AlgorithmStrategyObject */
   virtual bool InitializeImpl(
      const OptionsList& options,
      const std::string& prefix
   ) = 0;

   /** Solve the primal dual system, given one right hand side.
    *
    *  If the flag allow_inexact is set to true, it is not necessary to
    *  solve the system to best accuracy; for example, we don't want
    *  iterative refinement during the computation of the second
    *  order correction.  On the other hand, if improve_solution is
    *  true, the solution given in res should be improved (here beta
    *  has to be zero, and res is assume to be the solution for the
    *  system using rhs, without the factor alpha...).
    *
    *  @return false, if a solution could not be computed (for
    *  example, when the Hessian regularization parameter becomes too
    *  large)
    */
   virtual bool Solve(
      Number                alpha,
      Number                beta,
      const IteratesVector& rhs,
      IteratesVector&       res,
      bool                  allow_inexact = false,
      bool                  improve_solution = false
   ) = 0;

private:
   /**@name Default Compiler Generated Methods
    * (Hidden to avoid implicit creation/calling).
    *
    * These methods are not implemented and
    * we do not want the compiler to implement
    * them for us, so we declare them private
    * and do not define them. This ensures that
    * they will not be implicitly created/called.
    */
   ///@{
   /** Default Assignment Operator */
   PDSystemSolver& operator=(
      const PDSystemSolver&
   );
   ///@}
};

} // namespace Ipopt

#endif