49 #include "ConstrainedOptPack_Types.hpp"
50 #include "ConstrainedOptPack_MatrixSymAddDelUpdateableWithOpNonsingular.hpp"
51 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
52 #include "AbstractLinAlgPack_SpVectorClass.hpp"
53 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
54 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
55 #include "AbstractLinAlgPack_MatrixOp.hpp"
56 #include "AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
57 #include "AbstractLinAlgPack_MatrixOpSerial.hpp"
58 #include "DenseLinAlgPack_DMatrixClass.hpp"
61 #include "StopWatchPack_stopwatch.hpp"
63 namespace ConstrainedOptPack {
65 namespace QPSchurPack {
76 #ifdef LINALGPACK_CHECK_RANGE
79 return this->operator[](i-1);
85 #ifdef LINALGPACK_CHECK_RANGE
88 return this->operator[](i-1);
197 virtual const DVectorSlice
g()
const = 0;
199 virtual const MatrixSymOp&
G()
const = 0;
201 virtual const MatrixOp&
A()
const = 0;
257 virtual const DVectorSlice
b_X()
const = 0;
260 virtual const GenPermMatrixSlice&
Q_R()
const = 0;
263 virtual const GenPermMatrixSlice&
Q_X()
const = 0;
266 virtual const MatrixSymOpNonsing&
Ko()
const = 0;
269 virtual const DVectorSlice
fo()
const = 0;
284 virtual void dump_qp( std::ostream& out );
344 virtual const MatrixOp&
A_bar()
const = 0;
364 const DVectorSlice& x,
size_type* j_viol, value_type* constr_val
365 ,value_type* viol_bnd_val, value_type* norm_2_constr, EBounds* bnd,
bool* can_ignore
404 typedef MatrixSymAddDelUpdateable
MSADU;
407 {
public:
TestFailed(
const std::string& what_arg) : std::logic_error(what_arg) {}};
424 ,MAX_RUNTIME_EXEEDED_FAIL
425 ,MAX_RUNTIME_EXEEDED_DUAL_FEAS
426 ,MAX_ALLOWED_STORAGE_EXCEEDED
427 ,INFEASIBLE_CONSTRAINTS
435 ,OUTPUT_BASIC_INFO = 1
436 ,OUTPUT_ITER_SUMMARY = 2
437 ,OUTPUT_ITER_STEPS = 3
439 ,OUTPUT_ITER_QUANTITIES = 5
534 const schur_comp_ptr_t& schur_comp = Teuchos::null
536 ,value_type max_real_runtime = 1e+20
537 ,value_type feas_tol = 1e-8
538 ,value_type loose_feas_tol = 1e-6
539 ,value_type dual_infeas_tol = 1e-12
540 ,value_type huge_primal_step = 1e+20
541 ,value_type huge_dual_step = 1e+20
542 ,value_type warning_tol = 1e-10
543 ,value_type error_tol = 1e-5
546 ,value_type iter_refine_opt_tol = 1e-12
547 ,value_type iter_refine_feas_tol = 1e-12
548 ,
bool iter_refine_at_solution =
true
549 ,
bool salvage_init_schur_comp =
true
550 ,MSADU::PivotTolerances
pivot_tols = MSADU::PivotTolerances( 1e-8,1e-11,1e-11 )
613 ,
size_type num_act_change,
const int ij_act_change[],
const EBounds bnds[]
615 ,DVectorSlice* x, SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve
636 ,
const MatrixOp *
A_bar
637 ,
const GenPermMatrixSlice *
Q_R
642 const MatrixSymOp&
G()
const
645 const MatrixOp*
A()
const
651 const GenPermMatrixSlice&
Q_R()
const
655 {
return *P_XF_hat_; }
658 {
return *P_plus_hat_; }
676 ,
const DVectorSlice& vs_rhs2, value_type beta
681 ,
const SpVectorSlice& sv_rhs2, value_type beta
687 const MatrixSymOp *G_;
689 const MatrixOp *A_bar_;
690 const GenPermMatrixSlice *Q_R_;
691 const GenPermMatrixSlice *P_XF_hat_;
692 const GenPermMatrixSlice *P_plus_hat_;
736 typedef MatrixSymAddDelUpdateable
MSADU;
754 const schur_comp_ptr_t &schur_comp
755 ,MSADU::PivotTolerances
pivot_tols = MSADU::PivotTolerances( 1e-6,1e-8,1e-8 )
770 ,
const EBounds bnds[],
bool test,
bool salvage_init_schur_comp
796 size_type ja, EBounds bnd_ja,
bool update_steps
798 ,
bool force_refactorization =
true
799 ,
bool allow_any_cond =
false );
810 ,
bool force_refactorization =
true,
bool allow_any_cond =
false );
820 int jd,
size_type ja, EBounds bnd_ja,
bool update_steps
831 const QP&
qp()
const;
915 const MatrixSymOpNonsing&
S_hat()
const;
917 const GenPermMatrixSlice&
P_XF_hat()
const;
919 const GenPermMatrixSlice&
P_FC_hat()
const;
923 const GenPermMatrixSlice&
Q_XD_hat()
const;
925 const DVectorSlice
d_hat()
const;
927 DVectorSlice
z_hat();
929 const DVectorSlice
z_hat()
const;
933 const DVectorSlice
p_z_hat()
const;
937 const DVectorSlice
mu_D_hat()
const;
965 typedef std::vector<int> ij_map_t;
967 typedef std::map<int,size_type> s_map_t;
969 typedef std::vector<EBounds> bnds_t;
971 typedef std::vector<int> l_fxfx_t;
973 typedef std::vector<size_type> P_row_t;
975 typedef std::vector<size_type> P_col_t;
993 DVector constr_norm_;
1004 GenPermMatrixSlice P_XF_hat_;
1005 P_row_t P_XF_hat_row_;
1006 P_row_t P_XF_hat_col_;
1015 GenPermMatrixSlice P_FC_hat_;
1016 P_row_t P_FC_hat_row_;
1017 P_row_t P_FC_hat_col_;
1025 GenPermMatrixSlice P_plus_hat_;
1026 P_row_t P_plus_hat_row_;
1027 P_row_t P_plus_hat_col_;
1033 GenPermMatrixSlice Q_XD_hat_;
1034 P_row_t Q_XD_hat_row_;
1035 P_row_t Q_XD_hat_col_;
1041 DVector p_mu_D_hat_;
1047 void assert_initialized()
const;
1053 void reinitialize_matrices(
bool test);
1059 bool remove_augmented_element(
1060 size_type sd,
bool force_refactorization
1061 ,MatrixSymAddDelUpdateable::EEigenValType eigen_val_drop
1063 ,
bool allow_any_cond );
1075 ,
bool print_S_hat =
true );
1083 enum EPDSteps { PICK_VIOLATED_CONSTRAINT, UPDATE_ACTIVE_SET, COMPUTE_SEARCH_DIRECTION
1084 , COMPUTE_STEP_LENGTHS, TAKE_STEP };
1099 ,
const DVectorSlice& vo, ActiveSet*
act_set, DVectorSlice* v
1107 virtual void set_x(
const ActiveSet&
act_set,
const DVectorSlice& v, DVectorSlice* x );
1111 const ActiveSet&
act_set,
const DVectorSlice& v
1112 ,SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve );
1119 ITER_REFINE_NOT_PERFORMED
1120 ,ITER_REFINE_ONE_STEP
1121 ,ITER_REFINE_NOT_NEEDED
1122 ,ITER_REFINE_IMPROVED
1123 ,ITER_REFINE_NOT_IMPROVED
1124 ,ITER_REFINE_CONVERGED
1138 ,
const value_type ao
1139 ,
const DVectorSlice *bo
1140 ,
const value_type aa
1141 ,
const DVectorSlice *ba
STANDARD_COMPOSITION_MEMBERS(MatrixSymAddDelUpdateableWithOpNonsingular, schur_comp)
Schur complement matrix object S_hat.
MatrixSymAddDelUpdateable MSADU
virtual size_type n() const =0
bool drop_constraint(int jd, std::ostream *out, EOutputLevel output_level, bool force_refactorization=true, bool allow_any_cond=false)
Drop a constraint from the active set then refactorize the schur complement (if forced).
static value_type DEGENERATE_MULT
Value for near degenerate lagrange multipliers.
virtual const DVectorSlice fo() const =0
size_type s_map(int ij) const
Map from a constraint or initially fixed variable to a row and column in the schur complement S_bar...
virtual const MatrixOp & A_bar() const =0
int ij_map(size_type s) const
Returns -i for row & column of S_bar for an initially fixed variable left out of Ko that became free ...
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))
const GenPermMatrixSlice & Q_R() const
MatrixSymAddDelUpdateable MSADU
Thrown if there is some numerical instability.
T & operator()(typename this_t::size_type i)
one based indexing
const MatrixSymOp & G() const
virtual const x_init_t & x_init() const =0
Return the status of a variable initially.
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)
Initialize.
EBounds bnd(size_type s) const
Return which bound is active for the active constraint.
virtual value_type get_bnd(size_type j, EBounds bnd) const =0
Return the bound for a constraint.
const MatrixSymOpNonsing & S_hat() const
size_type q_F_hat() const
Return the number of variables that where initially fixed but are currently free or fixed to another ...
Represents the matrix U_hat.
virtual void set_multipliers(const ActiveSet &act_set, const DVectorSlice &v, SpVector *mu, DVectorSlice *lambda, SpVector *lambda_breve)
Map from the active set to the sparse multipliers for the inequality constraints. ...
MSADU::PivotTolerances pivot_tols() const
size_type q_C_hat() const
Return the number of variables that where initially fixed but are currently fixed to another bound...
const GenPermMatrixSlice & P_FC_hat() const
STANDARD_MEMBER_COMPOSITION_MEMBERS(size_type, max_iter)
Set the maximum number of primal-dual QP iterations to take.
bool is_init_fixed(size_type j) const
Determine if a constriant was an initially fixed variable.
void refactorize_schur_comp()
Reinitialize the schur complement factorization for the current active set.
static void dump_act_set_quantities(const ActiveSet &act_set, std::ostream &out, bool print_S_hat=true)
Dump all the active set quantities for debugging.
DVectorSlice p_mu_D_hat()
vector_one_based_checked< size_type > l_x_X_map_t
EOutputLevel
Output level.
const GenPermMatrixSlice & P_plus_hat() const
virtual const l_x_X_map_t & l_x_X_map() const =0
Map from full x(i) to initially fixed x_X(l).
virtual EPickPolicy pick_violated_policy() const =0
Represents and manages the active set for the QPSchur algorithm.
T operator()(typename this_t::size_type i) const
one based indexing
virtual const DVectorSlice g() const =0
virtual const GenPermMatrixSlice & Q_X() const =0
(Q_X().ordered_by() == BY_ROW)
const GenPermMatrixSlice & P_XF_hat() const
U_hat_t()
Construct uninitialized.
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)
Add a constraint to the active set then refactorize the schur complemnt (if forced).
const U_hat_t & U_hat() const
const MatrixOp * A() const
virtual const MatrixSymOpNonsing & Ko() const =0
virtual size_type m_breve() const =0
const GenPermMatrixSlice & P_plus_hat() const
Interface for updating a symmetric matrix and its factorization by adding and deleting rows and colum...
Represents the QP to be solved by QPSchur {abstract}.
virtual void set_x(const ActiveSet &act_set, const DVectorSlice &v, DVectorSlice *x)
Set the values in x for all the variables.
value_type constr_norm(size_type s) const
Returns ||a(j)||2 where j = ij_map(s).
virtual const MatrixOp & A() const =0
If m == 0 then don't call this, it may throw an exception or worse.
bool timeout_return(StopWatchPack::stopwatch *timer, std::ostream *out, EOutputLevel output_level) const
Determine if time has run out and if we should return.
size_type q_D_hat() const
Return the number of variables that where initially fixed and are still currently fixed to their inti...
virtual const GenPermMatrixSlice & Q_R() const =0
(Q_R().ordered_by() == BY_ROW)
const MatrixOp & A_bar() const
virtual Constraints & constraints()=0
STANDARD_COMPOSITION_MEMBERS(MatrixSymAddDelUpdateableWithOpNonsingular, schur_comp)
«std comp» members for schur complement matrix S_hat.
size_type q_plus_hat() const
Return the number of constraints from A_bar added to the active set.
virtual void ignore(size_type j)=0
Inform to ignore the jth constraint the next time pick_violated(...) is called.
virtual const i_x_X_map_t & i_x_X_map() const =0
Map from initially fixed x_X(l) to full x(i).
ESolveReturn
solve_qp return values
bool drop_add_constraints(int jd, size_type ja, EBounds bnd_ja, bool update_steps, std::ostream *out, EOutputLevel output_level)
Drop a constraint from, then add a constraint to the active set and refactorize the schur complement...
size_type q_hat() const
Return the total size of the schur complement.
Thrown if constraints are inconsistant (no feasible region)
bool all_dof_used_up() const
Returns true if all the degrees of freedom of the QP are used up.
const ActiveSet & act_set() const
Return a reference to the active set object.
virtual size_type m() const =0
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)
Solve a QP.
QPSchurPack::Constraints Constraints
vector_one_based_checked< EBounds > x_init_t
Utility class for a ranged check vector.
size_type l_fxfx(size_type k) const
Returns the indice of x_X(l) of the initially fixed variables that are still fixed at their original ...
const DVectorSlice d_hat() const
virtual void dump_qp(std::ostream &out)
Dump the definition of the QP to a stream.
virtual size_type n() const =0
const GenPermMatrixSlice & Q_XD_hat() const
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)
Run the algorithm from a dual feasible point.
ERunTests
Enumeration for if to run internal tests or not.
Represents the extra constraints in the QP to be satisfied by the schur complement QP solver QPSchur ...
const GenPermMatrixSlice & P_XF_hat() const
virtual const MatrixSymOp & G() const =0
EIterRefineReturn iter_refine(const ActiveSet &act_set, std::ostream *out, EOutputLevel output_level, const value_type ao, const DVectorSlice *bo, const value_type aa, const DVectorSlice *ba, DVectorSlice *v, DVectorSlice *z, size_type *iter_refine_num_resid, size_type *iter_refine_num_solves)
Perform iterative refinement on the augmented KKT system for the current active set.
vector_one_based_checked< size_type > i_x_X_map_t
Solves a Quadratic Program with a dual QP method using a schur complement factorization.
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)
Initialize with an additional active set.
Thrown if during the course of the primal-dual iteration a non-dual feasible point if found...
STANDARD_MEMBER_COMPOSITION_MEMBERS(MSADU::PivotTolerances, pivot_tols)
Set the tolerances to use when updating the schur complement.
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
Pick a violated constraint.
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
virtual const DVectorSlice b_X() const =0
virtual size_type n_R() const =0