ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_QPSchur.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef QPSCHUR_H
43 #define QPSCHUR_H
44 
45 #include <ostream>
46 #include <map>
47 #include <vector>
48 
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"
62 
63 namespace ConstrainedOptPack {
64 
65 namespace QPSchurPack {
66 
68 template < class T >
69 class vector_one_based_checked : public std::vector<T>
70 {
72 public:
74  T& operator()( typename this_t::size_type i )
75  {
76 #ifdef LINALGPACK_CHECK_RANGE
77  return this->at(i-1);
78 #else
79  return this->operator[](i-1);
80 #endif
81  }
83  T operator()( typename this_t::size_type i ) const
84  {
85 #ifdef LINALGPACK_CHECK_RANGE
86  return this->at(i-1);
87 #else
88  return this->operator[](i-1);
89 #endif
90  }
91 }; // end class vector_one_based
92 
93 class Constraints;
94 class QP;
95 
166 class QP {
167 public:
168 
169  // /////////////////
170  // Public Types
171 
178 
181 
182  // /////////////////
183  // Public Interface
184 
186  virtual ~QP()
187  {}
188 
189  // ///////////////////////////////////////
190  // Initial active set independent members
191 
193  virtual size_type n() const = 0;
195  virtual size_type m() const = 0;
197  virtual const DVectorSlice g() const = 0;
199  virtual const MatrixSymOp& G() const = 0;
201  virtual const MatrixOp& A() const = 0;
202 
203  // /////////////////////////////////////
204  // Initial active set specific members
205 
207  virtual size_type n_R() const = 0;
208 
219  virtual const x_init_t& x_init() const = 0;
220 
232  virtual const l_x_X_map_t& l_x_X_map() const = 0;
233 
243  virtual const i_x_X_map_t& i_x_X_map() const = 0;
244 
246  /* The bounds of the initially fixed variables.
247  *
248  * For 1 <= l <= n_X:
249  *
250  \verbatim
251  / xl(i_x_X_map(l)) : if x_init(i_x_X_map(l)) == LOWER
252  b_X(l) = | xu(i_x_X_map(l)) : if x_init(i_x_X_map(l)) == UPPER
253  \ xl(i_x_X_map(l)) = xu(i_x_X_map(l)) : if x_init(i_x_X_map(l)) == EQUALITY
254  \endverbatim
255  *
256  */
257  virtual const DVectorSlice b_X() const = 0;
258 
260  virtual const GenPermMatrixSlice& Q_R() const = 0;
261 
263  virtual const GenPermMatrixSlice& Q_X() const = 0;
264 
266  virtual const MatrixSymOpNonsing& Ko() const = 0;
267 
269  virtual const DVectorSlice fo() const = 0;
270 
271  // //////////////////////////////////////////////////////////
272  // Additional constaints for cl_bar <= A_bar'*x <= cu_bar
273 
275  virtual Constraints& constraints() = 0;
276 
278  virtual const Constraints& constraints() const = 0;
279 
284  virtual void dump_qp( std::ostream& out );
285 
286 }; // end class QP
287 
328 class Constraints {
329 public:
330 
332  enum EPickPolicy { ANY_VIOLATED, MOST_VIOLATED };
333 
335  virtual ~Constraints() {}
336 
338  virtual size_type n() const = 0;
339 
341  virtual size_type m_breve() const = 0;
342 
344  virtual const MatrixOp& A_bar() const = 0;
345 
347  virtual void pick_violated_policy( EPickPolicy pick_policy ) = 0;
349  virtual EPickPolicy pick_violated_policy() const = 0;
350 
363  virtual void pick_violated(
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
366  ) const = 0;
367 
370  virtual void ignore( size_type j ) = 0;
371 
382  virtual value_type get_bnd( size_type j, EBounds bnd ) const = 0;
383 
384 }; // end class Constraints
385 
386 } // end namespace QPSchurPack
387 
395 class QPSchur {
396 public:
397 
400 
404  typedef MatrixSymAddDelUpdateable MSADU;
406  class TestFailed : public std::logic_error
407  {public: TestFailed(const std::string& what_arg) : std::logic_error(what_arg) {}};
409  class InconsistantConstraintsException : public std::logic_error
410  {public: InconsistantConstraintsException(const std::string& what_arg) : std::logic_error(what_arg) {}};
412  class NumericalInstabilityException : public std::runtime_error
413  {public: NumericalInstabilityException(const std::string& what_arg) : std::runtime_error(what_arg) {}};
416  {public: DualInfeasibleException(const std::string& what_arg)
417  : NumericalInstabilityException(what_arg) {}};
419  enum ERunTests { RUN_TESTS, NO_TESTS };
422  OPTIMAL_SOLUTION
423  ,MAX_ITER_EXCEEDED
424  ,MAX_RUNTIME_EXEEDED_FAIL
425  ,MAX_RUNTIME_EXEEDED_DUAL_FEAS
426  ,MAX_ALLOWED_STORAGE_EXCEEDED
427  ,INFEASIBLE_CONSTRAINTS
428  ,NONCONVEX_QP
429  ,DUAL_INFEASIBILITY
430  ,SUBOPTIMAL_POINT
431  };
434  NO_OUTPUT = 0
435  ,OUTPUT_BASIC_INFO = 1
436  ,OUTPUT_ITER_SUMMARY = 2
437  ,OUTPUT_ITER_STEPS = 3
438  ,OUTPUT_ACT_SET = 4
439  ,OUTPUT_ITER_QUANTITIES = 5
440  };
442  static value_type DEGENERATE_MULT;
443 
445 
448 
451 
455 
458  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, max_real_runtime );
459 
462  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, feas_tol );
463 
466  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, loose_feas_tol );
467 
471  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, dual_infeas_tol );
472 
477  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, huge_primal_step );
478 
483  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, huge_dual_step );
484 
487  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, warning_tol );
488 
491  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, error_tol );
492 
496  STANDARD_MEMBER_COMPOSITION_MEMBERS( size_type, iter_refine_min_iter );
497 
501  STANDARD_MEMBER_COMPOSITION_MEMBERS( size_type, iter_refine_max_iter );
502 
506  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, iter_refine_opt_tol );
507 
511  STANDARD_MEMBER_COMPOSITION_MEMBERS( value_type, iter_refine_feas_tol );
512 
516  STANDARD_MEMBER_COMPOSITION_MEMBERS( bool, iter_refine_at_solution );
517 
521  STANDARD_MEMBER_COMPOSITION_MEMBERS( bool, salvage_init_schur_comp );
522 
525  void pivot_tols( MSADU::PivotTolerances pivot_tols );
527  MSADU::PivotTolerances pivot_tols() const;
528 
530  virtual ~QPSchur() {}
531 
533  QPSchur(
534  const schur_comp_ptr_t& schur_comp = Teuchos::null
535  ,size_type max_iter = 100
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
544  ,size_type iter_refine_min_iter = 1
545  ,size_type iter_refine_max_iter = 3
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 )
551  );
552 
611  virtual ESolveReturn solve_qp(
612  QP& qp
613  ,size_type num_act_change, const int ij_act_change[], const EBounds bnds[]
614  ,std::ostream *out, EOutputLevel output_level, ERunTests test_what
615  ,DVectorSlice* x, SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve
616  ,size_type* iter, size_type* num_adds, size_type* num_drops
617  );
618 
620 
628  class U_hat_t : public MatrixOpSerial {
629  public:
631  U_hat_t();
633  void initialize(
634  const MatrixSymOp *G
635  ,const MatrixOp *A
636  ,const MatrixOp *A_bar
637  ,const GenPermMatrixSlice *Q_R
638  ,const GenPermMatrixSlice *P_XF_hat
639  ,const GenPermMatrixSlice *P_plus_hat
640  );
642  const MatrixSymOp& G() const
643  { return *G_; }
645  const MatrixOp* A() const
646  { return A_; }
648  const MatrixOp& A_bar() const
649  { return *A_bar_; }
651  const GenPermMatrixSlice& Q_R() const
652  { return *Q_R_; }
654  const GenPermMatrixSlice& P_XF_hat() const
655  { return *P_XF_hat_; }
657  const GenPermMatrixSlice& P_plus_hat() const
658  { return *P_plus_hat_; }
659 
661 
662 
664  size_type rows() const;
666  size_type cols() const;
667 
669 
672 
674  void Vp_StMtV(
675  DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
676  ,const DVectorSlice& vs_rhs2, value_type beta
677  ) const;
679  void Vp_StMtV(
680  DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
681  ,const SpVectorSlice& sv_rhs2, value_type beta
682  ) const;
683 
685 
686  private:
687  const MatrixSymOp *G_;
688  const MatrixOp *A_;
689  const MatrixOp *A_bar_;
690  const GenPermMatrixSlice *Q_R_;
691  const GenPermMatrixSlice *P_XF_hat_;
692  const GenPermMatrixSlice *P_plus_hat_;
693 
694  }; // end class U_hat_t
695 
727  class ActiveSet {
728  public:
729 
730  // /////////////////////
731  // Public types
732 
736  typedef MatrixSymAddDelUpdateable MSADU;
737 
738  // /////////////////////
739  // Public interface
740 
747 
750  STANDARD_MEMBER_COMPOSITION_MEMBERS( MSADU::PivotTolerances, pivot_tols );
751 
753  ActiveSet(
754  const schur_comp_ptr_t &schur_comp
755  ,MSADU::PivotTolerances pivot_tols = MSADU::PivotTolerances( 1e-6,1e-8,1e-8 )
756  );
757 
760 
768  void initialize(
769  QP& qp, size_type num_act_change, const int ij_act_change[]
770  ,const EBounds bnds[], bool test, bool salvage_init_schur_comp
771  ,std::ostream *out, EOutputLevel output_level );
772 
777  void refactorize_schur_comp();
778 
795  bool add_constraint(
796  size_type ja, EBounds bnd_ja, bool update_steps
797  ,std::ostream *out, EOutputLevel output_level
798  ,bool force_refactorization = true
799  ,bool allow_any_cond = false );
800 
808  bool drop_constraint(
809  int jd, std::ostream *out, EOutputLevel output_level
810  ,bool force_refactorization = true, bool allow_any_cond = false );
811 
820  int jd, size_type ja, EBounds bnd_ja, bool update_steps
821  ,std::ostream *out, EOutputLevel output_level );
822 
824 
827 
829  QP& qp();
831  const QP& qp() const;
832 
834 
837 
842  size_type q_hat() const;
843 
847  size_type q_plus_hat() const;
848 
853  size_type q_F_hat() const;
854 
859  size_type q_C_hat() const;
860 
865  size_type q_D_hat() const;
866 
874  int ij_map( size_type s ) const;
875 
891  size_type s_map( int ij ) const;
892 
899  value_type constr_norm( size_type s ) const;
900 
903  EBounds bnd( size_type s ) const;
904 
910  size_type l_fxfx( size_type k ) const;
911 
913  const U_hat_t& U_hat() const;
915  const MatrixSymOpNonsing& S_hat() const;
917  const GenPermMatrixSlice& P_XF_hat() const;
919  const GenPermMatrixSlice& P_FC_hat() const;
921  const GenPermMatrixSlice& P_plus_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;
931  DVectorSlice p_z_hat();
933  const DVectorSlice p_z_hat() const;
935  DVectorSlice mu_D_hat();
937  const DVectorSlice mu_D_hat() const;
939  DVectorSlice p_mu_D_hat();
941  const DVectorSlice p_mu_D_hat() const;
942 
952  bool is_init_fixed( size_type j ) const;
953 
955  bool all_dof_used_up() const;
956 
958 
959  private:
960 
961  // ///////////////////////////
962  // Private types
963 
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;
976 
977  // ///////////////////////////
978  // Private data members
979 
980  bool initialized_;
981  bool test_;
982  QP* qp_; // QP being solved.
983  const QP::x_init_t *x_init_;
984  size_type n_;
985  size_type n_R_;
986  size_type m_;
987  size_type m_breve_;
988  size_type q_plus_hat_;
989  size_type q_F_hat_;
990  size_type q_C_hat_;
991  ij_map_t ij_map_;
992 // s_map_t s_map_;
993  DVector constr_norm_;
994  bnds_t bnds_;
995  l_fxfx_t l_fxfx_;
996  U_hat_t U_hat_;
997  //
998  // for s = 1...q_hat
999  //
1000  // / e(i) if i > 0 (where: i = -ij_map(s))
1001  // [P_XF_hat](:,s) = |
1002  // \ 0 otherwise
1003  //
1004  GenPermMatrixSlice P_XF_hat_; // \hat{P}^{XF} \in \Re^{n \times \hat{q}}
1005  P_row_t P_XF_hat_row_; // i
1006  P_row_t P_XF_hat_col_; // s
1007  //
1008  // for s = 1...q_hat
1009  //
1010  // / e(sd) if 0 < j <= n && is_init_fixed(j)
1011  // | (where: j = ij_map(s), sd = s_map(-j))
1012  // [P_FC_hat](:,s) = |
1013  // \ 0 otherwise
1014  //
1015  GenPermMatrixSlice P_FC_hat_; // {\tilde{P}^{F}}^{T} \hat{P}^{C} \in \Re^{\hat{q} \times \hat{q}}
1016  P_row_t P_FC_hat_row_; // sd
1017  P_row_t P_FC_hat_col_; // s
1018  //
1019  // for s = 1...q_hat
1020  //
1021  // / e(j) if j > 0 && !is_init_fixed(j) (where: j = ij_map(s))
1022  // [P_plus_hat](:,s) = |
1023  // \ 0 otherwise
1024  //
1025  GenPermMatrixSlice P_plus_hat_; // \hat{P}^{(+)} \in \Re^{n+\breve{m} \times \hat{q}^{D}}
1026  P_row_t P_plus_hat_row_; // j
1027  P_row_t P_plus_hat_col_; // s
1028  //
1029  // for k = 1...q_D_hat
1030  //
1031  // [Q_XD_hat](:,k) = e(i) (where is_init_fixed(i) && s_map(-i) == 0)
1032  //
1033  GenPermMatrixSlice Q_XD_hat_; // \hat{Q}^{XD} \in \Re^{n_X \times \hat{q}^{D}}
1034  P_row_t Q_XD_hat_row_; // i
1035  P_row_t Q_XD_hat_col_; // k
1036  //
1037  DVector d_hat_; // \hat{d}
1038  DVector z_hat_; // \hat{z}
1039  DVector p_z_hat_;
1040  DVector mu_D_hat_; // \hat{\mu}^{D}
1041  DVector p_mu_D_hat_; // p^{\hat{\mu}^{D}}
1042 
1043  // ///////////////////////////
1044  // Private member functions
1045 
1046  //
1047  void assert_initialized() const;
1048 
1049  // Assert in range.
1050  void assert_s( size_type s) const;
1051 
1052  // Reinitialize P_XF_hat, P_plus_hat, Q_XD_hat, and U_hat
1053  void reinitialize_matrices(bool test);
1054 
1055  // Remove an element from the augmented KKT system.
1056  // This does not update P_plus_hat, P_XF_hat or any
1057  // of the dimensions. Returns true if *out was
1058  // written to.
1059  bool remove_augmented_element(
1060  size_type sd, bool force_refactorization
1061  ,MatrixSymAddDelUpdateable::EEigenValType eigen_val_drop
1062  ,std::ostream *out, EOutputLevel output_level
1063  ,bool allow_any_cond );
1064 
1065  // not defined and not to be called.
1066  ActiveSet();
1067 
1068  }; // end class ActiveSet
1069 
1071  const ActiveSet& act_set() const;
1072 
1074  static void dump_act_set_quantities( const ActiveSet& act_set, std::ostream& out
1075  , bool print_S_hat = true );
1076 
1077 protected:
1078 
1079  // /////////////////////////
1080  // Protected types
1081 
1083  enum EPDSteps { PICK_VIOLATED_CONSTRAINT, UPDATE_ACTIVE_SET, COMPUTE_SEARCH_DIRECTION
1084  , COMPUTE_STEP_LENGTHS, TAKE_STEP };
1085 
1086  // ///////////////////////////
1087  // Protected Member functions
1088 
1095  virtual
1097  EPDSteps first_step
1098  ,std::ostream *out, EOutputLevel output_level, ERunTests test_what
1099  ,const DVectorSlice& vo, ActiveSet* act_set, DVectorSlice* v
1100  ,DVectorSlice* x, size_type* iter, size_type* num_adds, size_type* num_drops
1101  ,size_type* iter_refine_num_resid, size_type* iter_refine_num_solves
1102  ,StopWatchPack::stopwatch* timer
1103  );
1104 
1107  virtual void set_x( const ActiveSet& act_set, const DVectorSlice& v, DVectorSlice* x );
1108 
1110  virtual void set_multipliers(
1111  const ActiveSet& act_set, const DVectorSlice& v
1112  ,SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve );
1113 
1115  bool timeout_return( StopWatchPack::stopwatch*timer, std::ostream *out, EOutputLevel output_level ) const;
1116 
1119  ITER_REFINE_NOT_PERFORMED // Did not even perform it (iter_refine_max_iter == 0)
1120  ,ITER_REFINE_ONE_STEP // Only performed one step and the status is not known.
1121  ,ITER_REFINE_NOT_NEEDED // Convergence tolerance was already satisfied
1122  ,ITER_REFINE_IMPROVED // Did not converge but it was improved
1123  ,ITER_REFINE_NOT_IMPROVED // Tried iterative refinement but no improvement
1124  ,ITER_REFINE_CONVERGED // Performed iterative refinement and converged!
1125  };
1135  const ActiveSet &act_set
1136  ,std::ostream *out
1137  ,EOutputLevel output_level
1138  ,const value_type ao // Only used if bo != NULL
1139  ,const DVectorSlice *bo // If NULL then assumed to be zero!
1140  ,const value_type aa // Only used if q_hat > 0
1141  ,const DVectorSlice *ba // If NULL then assumed to be zero! Not accessed if q_hat > 0
1142  ,DVectorSlice *v
1143  ,DVectorSlice *z // Can be NULL if q_hat > 0
1144  ,size_type *iter_refine_num_resid
1145  ,size_type *iter_refine_num_solves
1146  );
1147 
1148 private:
1149 
1150  // /////////////////////////
1151  // Private data members
1152 
1153  ActiveSet act_set_; // The active set.
1154 
1155 }; // end class QPSchur
1156 
1157 } // end namespace ConstrainedOptPack
1158 
1159 #endif // QPSCHUR_H
STANDARD_COMPOSITION_MEMBERS(MatrixSymAddDelUpdateableWithOpNonsingular, schur_comp)
Schur complement matrix object S_hat.
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
T & operator()(typename this_t::size_type i)
one based indexing
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.
size_type q_F_hat() const
Return the number of variables that where initially fixed but are currently free or fixed to another ...
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.
vector_one_based_checked< size_type > l_x_X_map_t
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
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).
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.
size_t size_type
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)
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).
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.
vector_one_based_checked< EBounds > x_init_t
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 ...
virtual void dump_qp(std::ostream &out)
Dump the definition of the QP to a stream.
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
Transp
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