ConstrainedOptPack_QPSchur.hpp

// @HEADER
// ***********************************************************************
// 
// Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
//                  Copyright (2003) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#ifndef QPSCHUR_H
#define QPSCHUR_H

#include <ostream>
#include <map>
#include <vector>

#include "ConstrainedOptPack_Types.hpp"
#include "ConstrainedOptPack_MatrixSymAddDelUpdateableWithOpNonsingular.hpp"
#include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
#include "AbstractLinAlgPack_SpVectorClass.hpp"
#include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
#include "AbstractLinAlgPack_MatrixSymOp.hpp"
#include "AbstractLinAlgPack_MatrixOp.hpp"
#include "AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
#include "AbstractLinAlgPack_MatrixOpSerial.hpp"
#include "DenseLinAlgPack_DMatrixClass.hpp"
#include "Teuchos_StandardCompositionMacros.hpp"
#include "Teuchos_StandardMemberCompositionMacros.hpp"
#include "StopWatchPack_stopwatch.hpp"

namespace ConstrainedOptPack {

namespace QPSchurPack {

template < class T >
class vector_one_based_checked : public std::vector<T>
{
  typedef vector_one_based_checked this_t;
public:
  T& operator()( typename this_t::size_type i )
  {
#ifdef LINALGPACK_CHECK_RANGE
      return this->at(i-1);
#else
      return this->operator[](i-1);
#endif
  }
  T operator()( typename this_t::size_type i ) const
  {
#ifdef LINALGPACK_CHECK_RANGE
      return this->at(i-1);
#else
      return this->operator[](i-1);
#endif
  }
}; // end class vector_one_based

class Constraints;
class QP;

class QP {
public:

  // /////////////////
  // Public Types

  typedef vector_one_based_checked<EBounds>   x_init_t;
  typedef vector_one_based_checked<size_type>   l_x_X_map_t;
  typedef vector_one_based_checked<size_type>   i_x_X_map_t;

  typedef QPSchurPack::Constraints        Constraints;

  // /////////////////
  // Public Interface

  virtual ~QP()
  {}

  // ///////////////////////////////////////
  // Initial active set independent members 

  virtual size_type n() const = 0;
  virtual size_type m() const = 0;
  virtual const DVectorSlice g() const = 0;
  virtual const MatrixSymOp& G() const = 0;
  virtual const MatrixOp& A() const = 0;

  // /////////////////////////////////////
  // Initial active set specific members

  virtual size_type n_R() const = 0;

  virtual const x_init_t& x_init() const = 0;

  virtual const l_x_X_map_t& l_x_X_map() const = 0;

  virtual const i_x_X_map_t& i_x_X_map() const = 0;

  /* The bounds of the initially fixed variables.
   *
   * For 1 <= l <= n_X:
   *
   \verbatim
             / xl(i_x_X_map(l))                     : if x_init(i_x_X_map(l)) == LOWER
   b_X(l) =  | xu(i_x_X_map(l))                     : if x_init(i_x_X_map(l)) == UPPER
             \ xl(i_x_X_map(l)) = xu(i_x_X_map(l))  : if x_init(i_x_X_map(l)) == EQUALITY
   \endverbatim
   *
   */
  virtual const DVectorSlice b_X() const = 0;

  virtual const GenPermMatrixSlice& Q_R() const = 0;

  virtual const GenPermMatrixSlice& Q_X() const = 0;

  virtual const MatrixSymOpNonsing& Ko() const = 0;

  virtual const DVectorSlice fo() const = 0;

  // //////////////////////////////////////////////////////////
  // Additional constaints for cl_bar <= A_bar'*x <= cu_bar

  virtual Constraints& constraints() = 0;

  virtual const Constraints& constraints() const = 0;

  virtual void dump_qp( std::ostream& out );

};  // end class QP

class Constraints {
public:

  enum EPickPolicy { ANY_VIOLATED, MOST_VIOLATED };

  virtual ~Constraints() {}

  virtual size_type n() const = 0;
  
  virtual size_type m_breve() const = 0;

  virtual const MatrixOp& A_bar() const = 0;
  
  virtual void pick_violated_policy( EPickPolicy pick_policy ) = 0;
  virtual EPickPolicy pick_violated_policy() const = 0;

  virtual void pick_violated(
     const DVectorSlice& x, size_type* j_viol, value_type* constr_val
    ,value_type* viol_bnd_val, value_type* norm_2_constr, EBounds* bnd, bool* can_ignore
    ) const = 0;

  virtual void ignore( size_type j ) = 0;

  virtual value_type get_bnd( size_type j, EBounds bnd ) const = 0;

};  // end class Constraints

} // end namespace QPSchurPack 

class QPSchur {
public:


  typedef QPSchurPack::QP               QP;
  typedef MatrixSymAddDelUpdateable     MSADU;
  class TestFailed : public std::logic_error
  {public: TestFailed(const std::string& what_arg) : std::logic_error(what_arg) {}};
  class InconsistantConstraintsException : public std::logic_error
  {public: InconsistantConstraintsException(const std::string& what_arg) : std::logic_error(what_arg) {}};
  class NumericalInstabilityException : public std::runtime_error
  {public: NumericalInstabilityException(const std::string& what_arg) : std::runtime_error(what_arg) {}};
  class DualInfeasibleException : public NumericalInstabilityException
  {public: DualInfeasibleException(const std::string& what_arg)
    : NumericalInstabilityException(what_arg) {}};
  enum ERunTests { RUN_TESTS, NO_TESTS };
  enum ESolveReturn {
    OPTIMAL_SOLUTION
    ,MAX_ITER_EXCEEDED
    ,MAX_RUNTIME_EXEEDED_FAIL
    ,MAX_RUNTIME_EXEEDED_DUAL_FEAS
    ,MAX_ALLOWED_STORAGE_EXCEEDED
    ,INFEASIBLE_CONSTRAINTS
    ,NONCONVEX_QP
    ,DUAL_INFEASIBILITY
    ,SUBOPTIMAL_POINT
  };
  enum EOutputLevel {
     NO_OUTPUT          = 0
    ,OUTPUT_BASIC_INFO      = 1
    ,OUTPUT_ITER_SUMMARY    = 2
    ,OUTPUT_ITER_STEPS      = 3
    ,OUTPUT_ACT_SET       = 4
    ,OUTPUT_ITER_QUANTITIES   = 5
  };
  static value_type DEGENERATE_MULT;



  STANDARD_COMPOSITION_MEMBERS( MatrixSymAddDelUpdateableWithOpNonsingular, schur_comp );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( size_type, max_iter );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, max_real_runtime );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, feas_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, loose_feas_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, dual_infeas_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, huge_primal_step );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, huge_dual_step );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, warning_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, error_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( size_type, iter_refine_min_iter );
    
  STANDARD_MEMBER_COMPOSITION_MEMBERS( size_type, iter_refine_max_iter );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, iter_refine_opt_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, iter_refine_feas_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( bool, iter_refine_at_solution );
    
  STANDARD_MEMBER_COMPOSITION_MEMBERS( bool, salvage_init_schur_comp );

  void pivot_tols( MSADU::PivotTolerances pivot_tols );
  MSADU::PivotTolerances pivot_tols() const;

  virtual ~QPSchur() {}

  QPSchur(
    const schur_comp_ptr_t&   schur_comp           = Teuchos::null
    ,size_type                max_iter             = 100
    ,value_type               max_real_runtime     = 1e+20
    ,value_type               feas_tol             = 1e-8
    ,value_type               loose_feas_tol       = 1e-6
    ,value_type               dual_infeas_tol      = 1e-12
    ,value_type               huge_primal_step     = 1e+20
    ,value_type               huge_dual_step       = 1e+20
    ,value_type               warning_tol          = 1e-10
    ,value_type               error_tol            = 1e-5
    ,size_type                iter_refine_min_iter = 1
    ,size_type                iter_refine_max_iter = 3
    ,value_type               iter_refine_opt_tol  = 1e-12
    ,value_type               iter_refine_feas_tol = 1e-12
    ,bool                     iter_refine_at_solution = true
    ,bool                     salvage_init_schur_comp = true
    ,MSADU::PivotTolerances   pivot_tols = MSADU::PivotTolerances( 1e-8,1e-11,1e-11 )
    );

  virtual ESolveReturn solve_qp(
    QP& qp
    ,size_type num_act_change, const int ij_act_change[], const EBounds bnds[]
    ,std::ostream *out, EOutputLevel output_level, ERunTests test_what
    ,DVectorSlice* x, SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve
    ,size_type* iter, size_type* num_adds, size_type* num_drops
    );


  class U_hat_t : public MatrixOpSerial {
  public:
    U_hat_t();
    void initialize( 
       const MatrixSymOp    *G
      ,const MatrixOp     *A
      ,const MatrixOp     *A_bar
      ,const GenPermMatrixSlice *Q_R
      ,const GenPermMatrixSlice *P_XF_hat
      ,const GenPermMatrixSlice *P_plus_hat
      );
    const MatrixSymOp& G() const
    { return *G_; }
    const MatrixOp* A() const
    { return A_;  }
    const MatrixOp& A_bar() const
    { return *A_bar_; }
    const GenPermMatrixSlice& Q_R() const
    { return *Q_R_; }
    const GenPermMatrixSlice& P_XF_hat() const
    { return *P_XF_hat_;  }
    const GenPermMatrixSlice& P_plus_hat() const
    { return *P_plus_hat_;  }
    


    size_type rows() const;
    size_type cols() const;



    void Vp_StMtV(
      DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
      ,const DVectorSlice& vs_rhs2, value_type beta
      ) const;
    void Vp_StMtV(
      DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
      ,const SpVectorSlice& sv_rhs2, value_type beta
      ) const;

    
  private:
    const MatrixSymOp     *G_;
    const MatrixOp          *A_;
    const MatrixOp          *A_bar_;
    const GenPermMatrixSlice  *Q_R_;
    const GenPermMatrixSlice  *P_XF_hat_;
    const GenPermMatrixSlice  *P_plus_hat_;

  };  // end class U_hat_t

  class ActiveSet {
  public:

    // /////////////////////
    // Public types

    typedef QPSchurPack::QP            QP;
    typedef MatrixSymAddDelUpdateable  MSADU;

    // /////////////////////
    // Public interface

    STANDARD_COMPOSITION_MEMBERS( MatrixSymAddDelUpdateableWithOpNonsingular, schur_comp );

    STANDARD_MEMBER_COMPOSITION_MEMBERS( MSADU::PivotTolerances, pivot_tols );

    ActiveSet(
      const schur_comp_ptr_t   &schur_comp
      ,MSADU::PivotTolerances  pivot_tols = MSADU::PivotTolerances( 1e-6,1e-8,1e-8 )
      );


    void initialize( 
      QP& qp, size_type num_act_change, const int ij_act_change[]
      ,const EBounds bnds[], bool test, bool salvage_init_schur_comp
      ,std::ostream *out, EOutputLevel output_level );

    void refactorize_schur_comp();

    bool add_constraint(
      size_type ja, EBounds bnd_ja, bool update_steps
      ,std::ostream *out, EOutputLevel output_level
      ,bool force_refactorization = true
      ,bool allow_any_cond = false );

    bool drop_constraint(
      int jd, std::ostream *out, EOutputLevel output_level
      ,bool force_refactorization = true, bool allow_any_cond = false );

    bool drop_add_constraints(
      int jd, size_type ja, EBounds bnd_ja, bool update_steps
      ,std::ostream *out, EOutputLevel output_level );



    QP& qp();
    const QP& qp() const;



    size_type q_hat() const;

    size_type q_plus_hat() const;

    size_type q_F_hat() const;

    size_type q_C_hat() const;

    size_type q_D_hat() const;

    int ij_map( size_type s ) const;

    size_type s_map( int ij ) const;

    value_type constr_norm( size_type s ) const;

    EBounds bnd( size_type s ) const;

    size_type l_fxfx( size_type k ) const;

    const U_hat_t& U_hat() const;
    const MatrixSymOpNonsing& S_hat() const;
    const GenPermMatrixSlice& P_XF_hat() const;
    const GenPermMatrixSlice& P_FC_hat() const;
    const GenPermMatrixSlice& P_plus_hat() const;
    const GenPermMatrixSlice& Q_XD_hat() const;
    const DVectorSlice d_hat() const;
    DVectorSlice z_hat();
    const DVectorSlice z_hat() const;
    DVectorSlice p_z_hat();
    const DVectorSlice p_z_hat() const;
    DVectorSlice mu_D_hat();
    const DVectorSlice mu_D_hat() const;
    DVectorSlice p_mu_D_hat();
    const DVectorSlice p_mu_D_hat() const;

    bool is_init_fixed( size_type j ) const;

    bool all_dof_used_up() const;


  private:

    // ///////////////////////////
    // Private types

    typedef std::vector<int>      ij_map_t;
    typedef std::map<int,size_type>   s_map_t;
    typedef std::vector<EBounds>    bnds_t;
    typedef std::vector<int>      l_fxfx_t;
    typedef std::vector<size_type>    P_row_t;
    typedef std::vector<size_type>    P_col_t;

    // ///////////////////////////
    // Private data members

    bool        initialized_;
    bool        test_;
    QP*         qp_;  // QP being solved.
    const QP::x_init_t  *x_init_;
    size_type     n_;
    size_type     n_R_;
    size_type     m_;
    size_type     m_breve_;
    size_type     q_plus_hat_;
    size_type     q_F_hat_;
    size_type     q_C_hat_;
    ij_map_t      ij_map_;
//    s_map_t       s_map_;
    DVector       constr_norm_;
    bnds_t        bnds_;
    l_fxfx_t            l_fxfx_;
    U_hat_t       U_hat_;
    //
    // for s = 1...q_hat
    //
    //                     /  e(i)    if i > 0 (where: i = -ij_map(s))
    // [P_XF_hat](:,s)   = |
    //                     \  0       otherwise
    //                     
    GenPermMatrixSlice  P_XF_hat_;    // \hat{P}^{XF} \in \Re^{n \times \hat{q}}
    P_row_t       P_XF_hat_row_;  // i
    P_row_t       P_XF_hat_col_;  // s
    //
    // for s = 1...q_hat
    //
    //                   /  e(sd)   if 0 < j <= n && is_init_fixed(j)
    //                   |          (where: j = ij_map(s), sd = s_map(-j))
    // [P_FC_hat](:,s) = |
    //                   \  0       otherwise
    //
    GenPermMatrixSlice  P_FC_hat_;    // {\tilde{P}^{F}}^{T} \hat{P}^{C} \in \Re^{\hat{q} \times \hat{q}}
    P_row_t       P_FC_hat_row_;  // sd
    P_row_t       P_FC_hat_col_;  // s
    //
    // for s = 1...q_hat
    //
    //                     /  e(j)    if j > 0 && !is_init_fixed(j) (where: j = ij_map(s))
    // [P_plus_hat](:,s) = |
    //                     \  0       otherwise
    //
    GenPermMatrixSlice  P_plus_hat_;  // \hat{P}^{(+)} \in \Re^{n+\breve{m} \times \hat{q}^{D}}
    P_row_t       P_plus_hat_row_;  // j
    P_row_t       P_plus_hat_col_;  // s
    //
    // for k = 1...q_D_hat
    //
    // [Q_XD_hat](:,k) = e(i)  (where is_init_fixed(i) && s_map(-i) == 0)
    //
    GenPermMatrixSlice  Q_XD_hat_;    // \hat{Q}^{XD} \in \Re^{n_X \times \hat{q}^{D}}
    P_row_t       Q_XD_hat_row_;  // i
    P_row_t       Q_XD_hat_col_;  // k
    //
    DVector       d_hat_;     // \hat{d}
    DVector       z_hat_;     // \hat{z}
    DVector       p_z_hat_;
    DVector       mu_D_hat_;    // \hat{\mu}^{D}
    DVector       p_mu_D_hat_;  // p^{\hat{\mu}^{D}}

    // ///////////////////////////
    // Private member functions

    //
    void assert_initialized() const;

    // Assert in range.
    void assert_s( size_type s) const;

    // Reinitialize P_XF_hat, P_plus_hat, Q_XD_hat, and U_hat
    void reinitialize_matrices(bool test);
    
    // Remove an element from the augmented KKT system.
    // This does not update P_plus_hat, P_XF_hat or any
    // of the dimensions.  Returns true if *out was
    // written to.
    bool remove_augmented_element(
      size_type sd, bool force_refactorization
      ,MatrixSymAddDelUpdateable::EEigenValType eigen_val_drop
      ,std::ostream *out, EOutputLevel output_level
      ,bool allow_any_cond );

    // not defined and not to be called.
    ActiveSet();

  };  // end class ActiveSet

  const ActiveSet& act_set() const;

  static void dump_act_set_quantities( const ActiveSet& act_set, std::ostream& out
    , bool print_S_hat = true );

protected:

  // /////////////////////////
  // Protected types

  enum EPDSteps { PICK_VIOLATED_CONSTRAINT, UPDATE_ACTIVE_SET, COMPUTE_SEARCH_DIRECTION
    , COMPUTE_STEP_LENGTHS, TAKE_STEP };

  // ///////////////////////////
  // Protected Member functions

  virtual
  ESolveReturn qp_algo(
    EPDSteps first_step
    ,std::ostream *out, EOutputLevel output_level, ERunTests test_what
    ,const DVectorSlice& vo, ActiveSet* act_set, DVectorSlice* v
    ,DVectorSlice* x, size_type* iter, size_type* num_adds, size_type* num_drops
    ,size_type* iter_refine_num_resid, size_type* iter_refine_num_solves
    ,StopWatchPack::stopwatch* timer
    );

  virtual void set_x( const ActiveSet& act_set, const DVectorSlice& v, DVectorSlice* x );

  virtual void set_multipliers(
    const ActiveSet& act_set, const DVectorSlice& v
    ,SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve );

  bool timeout_return( StopWatchPack::stopwatch*timer, std::ostream *out, EOutputLevel output_level ) const;

  enum EIterRefineReturn {
    ITER_REFINE_NOT_PERFORMED    // Did not even perform it (iter_refine_max_iter == 0)
    ,ITER_REFINE_ONE_STEP        // Only performed one step and the status is not known.
    ,ITER_REFINE_NOT_NEEDED      // Convergence tolerance was already satisfied
    ,ITER_REFINE_IMPROVED        // Did not converge but it was improved
    ,ITER_REFINE_NOT_IMPROVED    // Tried iterative refinement but no improvement
    ,ITER_REFINE_CONVERGED       // Performed iterative refinement and converged!
  };
    EIterRefineReturn iter_refine(
    const ActiveSet      &act_set
    ,std::ostream        *out
    ,EOutputLevel        output_level
    ,const value_type    ao  // Only used if bo != NULL
    ,const DVectorSlice   *bo // If NULL then assumed to be zero!
    ,const value_type    aa  // Only used if q_hat > 0
    ,const DVectorSlice   *ba // If NULL then assumed to be zero!  Not accessed if q_hat > 0
    ,DVectorSlice         *v
    ,DVectorSlice         *z  // Can be NULL if q_hat > 0
    ,size_type           *iter_refine_num_resid
    ,size_type           *iter_refine_num_solves
    );

private:

  // /////////////////////////
  // Private data members

  ActiveSet  act_set_;  // The active set.

};  // end class QPSchur

} // end namespace ConstrainedOptPack 

#endif  // QPSCHUR_H