76 #include "ConstrainedOptPack_MatrixSymPosDefLBFGS.hpp"
77 #include "ConstrainedOptPack/src/AbstractLinAlgPack_BFGS_helpers.hpp"
78 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
79 #include "DenseLinAlgPack_DMatrixOut.hpp"
80 #include "DenseLinAlgLAPack.hpp"
84 using DenseLinAlgPack::DVectorSlice;
94 void comp_Cb(
const DMatrixSlice& Lb,
const DVectorSlice& Db_diag
99 namespace ConstrainedOptPack {
105 const DMatrixSliceTri MatrixSymPosDefLBFGS::R()
const
107 return DenseLinAlgPack::tri( STY_(1,m_bar_,1,m_bar_), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
111 const DMatrixSliceTri MatrixSymPosDefLBFGS::Lb()
const
113 return DenseLinAlgPack::tri( STY_(2,m_bar_,1,m_bar_-1), BLAS_Cpp::lower, BLAS_Cpp::nonunit );
117 DMatrixSlice MatrixSymPosDefLBFGS::STY()
119 return STY_(1,m_bar_,1,m_bar_);
123 const DMatrixSlice MatrixSymPosDefLBFGS::STY()
const
125 return STY_(1,m_bar_,1,m_bar_);
129 DMatrixSliceSym MatrixSymPosDefLBFGS::STS()
131 return DenseLinAlgPack::nonconst_sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
135 const DMatrixSliceSym MatrixSymPosDefLBFGS::STS()
const
137 return DenseLinAlgPack::sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
141 DMatrixSliceSym MatrixSymPosDefLBFGS::YTY()
143 return DenseLinAlgPack::nonconst_sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
147 const DMatrixSliceSym MatrixSymPosDefLBFGS::YTY()
const
149 return DenseLinAlgPack::sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
158 ,
bool maintain_original
159 ,
bool maintain_inverse
163 initial_setup(max_size,m,maintain_original,maintain_inverse,auto_rescaling);
169 ,
bool maintain_original
170 ,
bool maintain_inverse
175 if( !maintain_original && !maintain_inverse )
176 throw std::invalid_argument(
177 "MatrixSymPosDefLBFGS::initial_setup(...) : "
178 "Error, both maintain_original and maintain_inverse can not both be false!" );
180 throw std::invalid_argument(
181 "MatrixSymPosDefLBFGS::set_num_updates_stored(m) : "
182 "Error, the number of storage locations must be > 0" );
188 num_secant_updates_= 0;
202 assert_initialized();
203 out <<
"*** Limited Memory BFGS matrix.\n"
204 <<
"Conversion to dense =\n";
205 MatrixOp::output(out);
206 out <<
"\n*** Stored quantities\n"
207 <<
"\ngamma_k = " << gamma_k_ << std::endl;
209 out <<
"\nS =\n" <<
S()
211 <<
"\nS'Y =\n" << STY_(1,m_bar_,1,m_bar_)
212 <<
"\nlower(S'S) \\ zero diagonal \\ upper(Y'Y) =\n"
213 << STSYTY_(1,m_bar_+1,1,m_bar_+1)
214 <<
"\nQ updated? = " << Q_updated_ << std::endl
215 <<
"\nCholesky of schur complement of Q, QJ =\n" << QJ_(1,m_bar_,1,m_bar_);
224 if( p_m ==
this )
return *
this;
226 auto_rescaling_ = p_m->auto_rescaling_;
227 maintain_original_ = p_m->maintain_original_;
228 original_is_updated_ = p_m->original_is_updated_;
229 maintain_inverse_ = p_m->maintain_inverse_;
230 inverse_is_updated_ = p_m->inverse_is_updated_;
231 n_max_ = p_m->n_max_;
234 m_bar_ = p_m->m_bar_;
235 k_bar_ = p_m->k_bar_;
236 num_secant_updates_ = p_m->num_secant_updates_;
237 gamma_k_ = p_m->gamma_k_;
241 STSYTY_ = p_m->STSYTY_;
242 Q_updated_ = p_m->Q_updated_;
246 throw std::invalid_argument(
"MatrixSymPosDefLBFGS::operator=(const MatrixOp& m)"
247 " : The concrete type of m is not a subclass of MatrixSymPosDefLBFGS as expected" );
256 ,
const DVectorSlice& x, value_type beta)
const
258 using LinAlgOpPack::V_StMtV;
259 using LinAlgOpPack::V_MtV;
265 assert_initialized();
289 invgk = 1.0 / gamma_k_;
307 if( work_.size() < 4 * m_ )
308 work_.resize( 4 * m_ );
320 t1 = work_( t1s, t1s + t1n - 1 ),
321 t2 = work_( t2s, t2s + t2n - 1 );
352 ,
const DVectorSlice& x )
const
359 using LinAlgOpPack::V_MtV;
360 using LinAlgOpPack::V_StMtV;
361 using LinAlgOpPack::Vp_MtV;
364 assert_initialized();
397 V_StV( y, gamma_k_, x );
407 if( work_.size() < 6*m_ )
408 work_.resize( 6*m_ );
423 t1 = work_( t1s, t1s + t1n - 1 ),
424 t2 = work_( t2s, t2s + t2n - 1 ),
425 t3 = work_( t3s, t3s + t3n - 1 ),
426 t4 = work_( t4s, t4s + t4n - 1 ),
427 t5 = work_( t5s, t5s + t5n - 1 );
433 const DMatrixSliceTri
436 const DMatrixSliceSym
448 V_mV( &t3, t1(mb+1,2*mb) );
477 std::ostringstream omsg;
479 <<
"MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
480 <<
"alpha = " << alpha <<
" <= 0 is not allowed!";
481 throw std::invalid_argument( omsg.str() );
486 else if( n > n_max_ ) {
487 std::ostringstream omsg;
489 <<
"MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
490 <<
"n = " << n <<
" > max_size = " << n_max_;
491 throw std::invalid_argument( omsg.str() );
495 S_.resize( n_max_, m_ );
496 Y_.resize( n_max_, m_ );
497 STY_.resize( m_, m_ );
498 STSYTY_.resize( m_+1, m_+1 );
499 STSYTY_.diag(0) = 0.0;
501 gamma_k_ = 1.0/alpha;
508 original_is_updated_ =
true;
509 inverse_is_updated_ =
true;
510 num_secant_updates_ = 0;
515 using DenseLinAlgPack::norm_inf;
520 DVectorSlice* s, DVectorSlice* y, DVectorSlice* Bs)
523 using DenseLinAlgPack::norm_2;
525 using LinAlgOpPack::V_MtV;
527 assert_initialized();
532 std::ostringstream omsg;
533 if( !BFGS_sTy_suff_p_d(*s,*y,&sTy,&omsg,
"MatrixSymPosDefLBFGS::secant_update(...)" ) ) {
534 throw UpdateSkippedException( omsg.str() );
546 {
for(
size_type k = 1; k <= m_-1; ++k ) {
547 S_.col(k) = S_.col(k+1);
548 Y_.col(k) = Y_.col(k+1);
549 STY_.col(k)(1,m_-1) = STY_.col(k+1)(2,m_);
550 STSYTY_.col(k)(k+1,m_) = STSYTY_.col(k+1)(k+2,m_+1);
551 STSYTY_.col(k+1)(1,k) = STSYTY_.col(k+2)(2,k+1);
563 S_.col(k_bar_)(1,n_) = *s;
564 Y_.col(k_bar_)(1,n_) = *y;
613 if( work_.size() < m_ )
620 STSYTY_.row(k_bar_+1)(1,k_bar_) = work_(1,k_bar_);
622 STSYTY_.col(k_bar_)(k_bar_+1,m_bar_+1) = work_(k_bar_,m_bar_);
648 STSYTY_.col(k_bar_+1)(1,k_bar_) = work_(1,k_bar_);
650 STSYTY_.row(k_bar_)(k_bar_+1,m_bar_+1) = work_(k_bar_,m_bar_);
657 gamma_k_ = STY_(k_bar_,k_bar_) / STSYTY_(k_bar_,k_bar_+1);
663 num_secant_updates_++;
683 std::ostringstream omsg;
685 <<
"MatrixSymPosDefLBFGS::initialize(alpha,max_size) : Error, "
686 <<
"alpha = " << alpha <<
" <= 0 is not allowed!";
687 throw std::invalid_argument( omsg.str() );
694 const DMatrixSliceSym &A
696 ,
bool force_factorization
698 ,PivotTolerances pivot_tols
701 throw std::runtime_error(
702 "MatrixSymPosDefLBFGS::initialize(A,max_size,force_refactorization,inertia) : Error, "
703 "This function is undefined for this subclass. I am so sorry for this terrible hack!" );
713 return Inertia(0,0,n_);
722 const DVectorSlice *t
724 ,
bool force_refactorization
725 ,EEigenValType add_eigen_val
726 ,PivotTolerances pivot_tols
729 assert_initialized();
731 std::ostringstream omsg;
733 <<
"MatrixSymPosDefLBFGS::augment_update(...) : Error, "
734 <<
"this->rows() = " << n_ <<
" == this->max_size() = " << n_max_
735 <<
" so we can't allow the matrix to grow!";
736 throw std::invalid_argument( omsg.str() );
739 throw std::invalid_argument(
740 "MatrixSymPosDefLBFGS::augment_update(...) : Error, "
741 "t must be NULL in this implemention. Sorry for this hack" );
744 std::ostringstream omsg;
746 <<
"MatrixSymPosDefLBFGS::augment_update(...) : Error, "
747 <<
"alpha = " << alpha <<
" <= 0 is not allowed!";
748 throw std::invalid_argument( omsg.str() );
750 if( add_eigen_val == MatrixSymAddDelUpdateable::EIGEN_VAL_NEG ) {
751 std::ostringstream omsg;
753 <<
"MatrixSymPosDefLBFGS::augment_update(...) : Error, "
754 <<
"add_eigen_val == EIGEN_VAL_NEG is not allowed!";
755 throw std::invalid_argument( omsg.str() );
770 S_.row(n_+1)(1,m_bar_) = 0.0;
771 Y_.row(n_+1)(1,m_bar_) = 0.0;
777 ,
bool force_refactorization
778 ,EEigenValType drop_eigen_val
779 ,PivotTolerances pivot_tols
782 assert_initialized();
823 DMatrixSlice S = this->
S();
824 DMatrixSlice
Y = this->
Y();
825 DMatrixSlice STY = this->STY();
826 DMatrixSliceSym STS = this->STS();
827 DMatrixSliceSym YTY = this->YTY();
829 DenseLinAlgPack::ger( -1.0, S.row(jd), Y.row(jd), &STY );
831 DenseLinAlgPack::syr( -1.0, S.row(jd), &STS );
833 DenseLinAlgPack::syr( -1.0, Y.row(jd), &YTY );
837 {
for(
size_type k = 1; k <= m_bar_; ++k ) {
838 value_type *ptr = S.col_ptr(k);
839 std::copy( ptr + jd, ptr + n_, ptr + jd - 1 );
841 {
for(
size_type k = 1; k <= m_bar_; ++k ) {
842 value_type *ptr = Y.col_ptr(k);
843 std::copy( ptr + jd, ptr + n_, ptr + jd - 1 );
853 void MatrixSymPosDefLBFGS::Vp_DtV( DVectorSlice* y,
const DVectorSlice& x )
const
855 DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), m_bar_, m_bar_
858 DVectorSlice::const_iterator
859 d_itr = STY_.diag(0).begin(),
861 DVectorSlice::iterator
864 while( y_itr != y->end() )
865 *y_itr++ += (*d_itr++) * (*x_itr++);
893 void MatrixSymPosDefLBFGS::update_Q()
const
895 using DenseLinAlgPack::tri;
896 using DenseLinAlgPack::tri_ele;
926 if( QJ_.rows() < m_ )
927 QJ_.resize( m_, m_ );
956 comp_Cb( STY_(2,mb,1,mb-1), STY_.diag(0)(1,mb-1), &C(2,mb,2,mb) );
960 const DMatrixSliceSym &STS = this->STS();
961 Mp_StM( &C, (1/gamma_k_), tri( STS.gms(), STS.uplo(), BLAS_Cpp::nonunit )
968 DMatrixSliceTriEle C_upper = tri_ele( C, BLAS_Cpp::upper );
969 DenseLinAlgLAPack::potrf( &C_upper );
974 void MatrixSymPosDefLBFGS::V_invQtV( DVectorSlice* x,
const DVectorSlice& y )
const
976 using DenseLinAlgPack::sym;
977 using DenseLinAlgPack::tri;
982 using LinAlgOpPack::V_MtV;
1000 x2 = (*x)(mb+1,2*mb);
1012 DVectorSlice::const_iterator
1013 d_itr = STY_.diag(0).begin(),
1014 y2_itr = y2.begin();
1015 DVectorSlice::iterator
1016 x1_itr = x1.begin();
1017 while( x1_itr != x1.end() )
1018 *x1_itr++ = *y2_itr++ / *d_itr++;
1035 std::copy( x1a.rbegin(), x1a.rend(), x1b.rbegin() );
1036 V_MtV( &x1b, Lb(), BLAS_Cpp::no_trans, x1b );
1044 const DMatrixSliceTri J = tri( QJ_(1,mb,1,mb), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
1048 V_InvMtV( &x1, J, BLAS_Cpp::no_trans, x1 );
1071 DVectorSlice::const_iterator
1072 d_itr = STY_.diag(0).begin();
1073 DVectorSlice::iterator
1074 x2_itr = x2.begin();
1075 while( x2_itr != x2.end() )
1076 *x2_itr++ /= *d_itr++;
1080 void MatrixSymPosDefLBFGS::assert_initialized()
const
1083 throw std::logic_error(
"MatrixSymPosDefLBFGS::assert_initialized() : "
1084 "Error, matrix not initialized" );
1091 void comp_Cb(
const DMatrixSlice& Lb,
const DVectorSlice& Db_diag
1092 , DMatrixSlice* Cb )
1112 typedef DenseLinAlgPack::value_type value_type;
1116 const size_type p = Db_diag.size();
1118 for( size_type i = 1; i <= p; ++i ) {
1119 for( size_type j = i; j <= p; ++j ) {
1120 value_type &c = (*Cb)(i,j) = 0.0;
1121 for( size_type k = 1; k <= i; ++k ) {
1122 c += Lb(i,k) * Lb(j,k) / Db_diag(k);
void init_identity(const VectorSpace &space_diag, value_type alpha)
bool maintain_original() const
void secant_update(VectorMutable *s, VectorMutable *y, VectorMutable *Bs)
void initial_setup(size_type m=10, bool maintain_original=true, bool maintain_inverse=true, bool auto_rescaling=false)
Initial setup for the matrix.
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
void initialize(value_type alpha, size_type max_size)
This is fine as long as alpha > 0.0.
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
MatrixSymPosDefLBFGS(size_type m=10, bool maintain_original=true, bool maintain_inverse=true, bool auto_rescaling=false)
Calls this->initial_setup()
void delete_update(size_type jd, bool force_refactorization, EEigenValType drop_eigen_val, PivotTolerances pivot_tols)
Should always succeed unless user gives a wrong value for drop_eigen_val.
void init_diagonal(const Vector &diag)
Actually this calls init_identity( diag.space(), diag.norm_inf() ).
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
const multi_vec_ptr_t S() const
const multi_vec_ptr_t Y() const
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
std::ostream & output(std::ostream &out) const
friend friend void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
Inertia inertia() const
Returns (0,0,rows())
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
bool maintain_inverse() const
MatrixOp & operator=(const MatrixOp &mwo)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void augment_update(const DVectorSlice *t, value_type alpha, bool force_refactorization, EEigenValType add_eigen_val, PivotTolerances pivot_tols)
Augment the matrix to add a row and column.
void set_uninitialized()
Will set rows() == 0.
size_type max_size() const
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)