ConstrainedOptPack_ConstraintsRelaxedStd.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.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#ifndef QP_SCHUR_CONSTRAINTS_RELAXED_STD_H
#define QP_SCHUR_CONSTRAINTS_RELAXED_STD_H

#include <list>

#include "ConstrainedOptPack_QPSchur.hpp"
#include "AbstractLinAlgPack_MatrixOp.hpp"
#include "AbstractLinAlgPack_VectorSpaceBlocked.hpp"

namespace ConstrainedOptPack {
namespace QPSchurPack {

class ConstraintsRelaxedStd : public Constraints {
public:

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

  class MatrixConstraints : public MatrixOp {
  public:

    MatrixConstraints();

    void initialize(
      const VectorSpace::space_ptr_t   &space_d_eta
      ,const size_type                 m_in
      ,const size_type                 m_eq
      ,const MatrixOp              *E
      ,BLAS_Cpp::Transp                trans_E
      ,const Vector              *b
      ,const MatrixOp              *F
      ,BLAS_Cpp::Transp                trans_F
      ,const Vector              *f
      ,size_type                       m_undecomp
      ,const size_type                 j_f_undecomp[]
      );


    size_type     nd() const
    { return nd_; }
    size_type     m_in() const
    { return m_in_; }
    size_type     m_eq() const
    { return m_eq_; }
    const MatrixOp* E() const
    { return E_;  }
    BLAS_Cpp::Transp  trans_E() const
    { return trans_E_;  }
    const Vector* b() const
    { return b_;  }
    const MatrixOp* F() const
    { return F_;  }
    BLAS_Cpp::Transp  trans_F() const
    { return trans_F_;  }
    const Vector* f() const
    { return f_;  }
    const GenPermMatrixSlice& P_u() const
    { return P_u_;  }



    const VectorSpace& space_cols() const;
    const VectorSpace& space_rows() const;
    MatrixOp& operator=(const MatrixOp& m);
//    ///
//    void Mp_StPtMtP(
//      DMatrixSlice* gms_lhs, value_type alpha
//      ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
//      ,BLAS_Cpp::Transp M_trans
//      ,const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
//      ) const ;
    void Vp_StMtV(
      VectorMutable* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
      ,const Vector& vs_rhs2, value_type beta
      ) const;
    void Vp_StPtMtV(
      VectorMutable* vs_lhs, value_type alpha
      ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
      ,BLAS_Cpp::Transp M_rhs2_trans
      ,const SpVectorSlice& sv_rhs3, value_type beta
      ) const;

    
  private:
    typedef std::vector<size_type>    row_i_t;
    typedef std::vector<size_type>    col_j_t;
    size_type     nd_;  // # unknowns d
    size_type     m_in_;  // # op(E)*d inequality constraints
    size_type     m_eq_;  // # op(F)*d equality constraints
    const MatrixOp  *E_;  // If NULL then no general inequalities
    BLAS_Cpp::Transp  trans_E_;
    const Vector  *b_;
    const MatrixOp  *F_;  // If NULL then no general equalities
    BLAS_Cpp::Transp  trans_F_;
    const Vector  *f_;
    GenPermMatrixSlice  P_u_;
    row_i_t             P_u_row_i_;
    col_j_t             P_u_col_j_;
    VectorSpace::space_ptr_t   space_cols_;
    VectorSpaceBlocked         space_rows_;
  };  // end class MatrixConstraints

  enum EInequalityPickPolicy {
    ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY
    ,ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY
    ,ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY
  };


  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, bounds_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, inequality_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, equality_tol );

  STANDARD_MEMBER_COMPOSITION_MEMBERS( EInequalityPickPolicy, inequality_pick_policy );

  ConstraintsRelaxedStd();

  void initialize(
    const VectorSpace::space_ptr_t   &space_d_eta
    ,value_type                      etaL
    ,const Vector              *dL
    ,const Vector              *dU
    ,const MatrixOp              *E
    ,BLAS_Cpp::Transp                trans_E
    ,const Vector              *b
    ,const Vector              *eL
    ,const Vector              *eU
    ,const MatrixOp              *F
    ,BLAS_Cpp::Transp                trans_F
    ,const Vector              *f
    ,size_type                       m_undecomp
    ,const size_type                 j_f_undecomp[]
    ,VectorMutable             *Ed
    ,bool                            check_F           = true
    ,value_type                      bounds_tol        = 1e-10
    ,value_type                      inequality_tol    = 1e-10
    ,value_type                      equality_tol      = 1e-10
    );

  const MatrixConstraints& A_bar_relaxed() const;



  size_type n() const;
  size_type m_breve() const;
  const MatrixOp& A_bar() const;
  void pick_violated_policy( EPickPolicy pick_policy );
  EPickPolicy pick_violated_policy() const;
  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;
  void ignore( size_type j );
  value_type get_bnd( size_type j, EBounds bnd ) const;


private:

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

  typedef std::list<size_type>  passed_over_equalities_t;

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

  MatrixConstraints   A_bar_;
  value_type          etaL_;
  const Vector  *dL_; // If NULL then no simple bounds
  const Vector  *dU_;
  const Vector  *eL_;
  const Vector  *eU_;
  VectorMutable *Ed_;
  bool                check_F_;
  mutable size_type   last_added_j_;          // Remember the last bound added so that
  mutable value_type  last_added_bound_;      // we can save our selfs some work.
  mutable EBounds     last_added_bound_type_; // ...
  mutable size_type   next_undecomp_f_k_;
       // Determines the next constraint [P_u'*op(F)*d + (1 - eta) * P_u'*f](next_undecomp_f_k)
      // to be checked to see if it is violated.  If next_undecomp_f_k > P_u.nz() then all
      // of the constriants have been checked.
  mutable passed_over_equalities_t passed_over_equalities_;
      // This is a list that keeps track of those equality constraints that were checked
      // for being violated but were within tolerance and therefore passed over and not added.
      // This list can be traversed again and again to check these constraints.  Specifically, the
      // indexes of f(k) are sorted, not the indexes in P_u'.

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

  void cache_last_added(
    size_type last_added_j, value_type last_added_bound
    ,EBounds last_added_bound_type
    ) const;

}; // end class ConstraintsRelaxedStd

} // end namespace QPSchurPack 
} // end namespace ConstrainedOptPack 

#endif // QP_SCHUR_CONSTRAINTS_RELAXED_STD_H