105 namespace ConstrainedOptPack {
163 ,
bool maintain_original
164 ,
bool maintain_inverse
168 initial_setup(m,maintain_original,maintain_inverse,auto_rescaling);
173 bool maintain_original,
174 bool maintain_inverse,
180 !maintain_original && !maintain_inverse, std::invalid_argument
181 ,
"MatrixSymPosDefLBFGS::initial_setup(...) : "
182 "Error, both maintain_original and maintain_inverse can not both be false!" );
184 m < 1, std::invalid_argument
185 ,
"MatrixSymPosDefLBFGS::set_num_updates_stored(m) : "
186 "Error, the number of storage locations must be > 0" );
193 auto_rescaling_ = auto_rescaling;
206 out <<
"\n*** Limited Memory BFGS matrix\n";
207 out <<
"\nConversion to dense =\n";
209 out <<
"\nStored quantities\n"
213 <<
"\ngamma_k = " <<
gamma_k_ << std::endl;
215 out <<
"\nS =\n" << *
S()
218 <<
"\nlower(S'S) \\ zero diagonal \\ upper(Y'Y) =\n"
220 <<
"\nQ updated? = " <<
Q_updated_ << std::endl;
222 out <<
"\nCholesky of schur complement of Q, QJ =\n" <<
QJ_(1,
m_bar_,1,
m_bar_);
231 if( p_m ==
this )
return *
this;
233 auto_rescaling_ = p_m->auto_rescaling_;
253 true,std::invalid_argument
254 ,
"MatrixSymPosDefLBFGS::operator=(const MatrixOp& mwo) : Error, "
255 "The concrete type of mwo \'" <<
typeName(mwo) <<
"\' is not "
256 "as subclass of MatrixSymPosDefLBFGS as required" );
273 typedef VectorDenseEncap vde;
274 typedef VectorDenseMutableEncap vdme;
322 Workspace<value_type> t1_ws(wss,2*mb);
324 Workspace<value_type> t2_ws(wss,2*mb);
327 VectorSpace::vec_mut_ptr_t
328 t = S->space_rows().create_member();
334 t1(1,mb) = vde(*t)();
336 t1(mb+1,2*mb) = vde(*t)();
342 (vdme(*t)() = t2(1,mb));
346 (vdme(*t)() = t2(mb+1,2*mb));
367 typedef VectorDenseEncap vde;
368 typedef VectorDenseMutableEncap vdme;
419 Workspace<value_type> t1_ws(wss,2*mb);
421 Workspace<value_type> t2_ws(wss,mb);
423 Workspace<value_type> t3_ws(wss,mb);
425 Workspace<value_type> t4_ws(wss,mb);
427 Workspace<value_type> t5_ws(wss,mb);
430 VectorSpace::vec_mut_ptr_t
431 t = S->space_rows().create_member();
442 t1(1,mb) = vde(*t)();
444 t1(mb+1,2*mb) = vde(*t)();
450 V_mV( &t3, t1(mb+1,2*mb) );
451 V_InvMtV( &t3, R, BLAS_Cpp::trans, t3 );
460 V_InvMtV( &t5, R, BLAS_Cpp::trans, t4 );
481 alpha <= 0.0, std::invalid_argument
482 ,
"MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
483 "alpha = " << alpha <<
" <= 0 is not allowed!" );
515 VectorMutable* s, VectorMutable* y, VectorMutable* Bs
529 std::ostringstream omsg;
530 if( !
BFGS_sTy_suff_p_d(*s,*y,&sTy,&omsg,
"MatrixSymPosDefLBFGS::secant_update(...)" ) ) {
531 throw UpdateSkippedException( omsg.str() );
542 S_->col(k) =
S_->col(k+1);
543 Y_->col(k) =
Y_->col(k+1);
581 VectorSpace::vec_mut_ptr_t
582 t = S->space_rows().create_member();
610 Workspace<value_type> work_ws(wss,
m_bar_);
615 work = VectorDenseEncap(*t)();
644 work = VectorDenseEncap(*t)();
681 DVectorSlice::const_iterator
682 d_itr =
STY_.diag(0).begin(),
684 DVectorSlice::iterator
687 while( y_itr != y->end() )
688 *y_itr++ += (*d_itr++) * (*x_itr++);
749 if(
QJ_.rows() <
m_ )
779 comp_Cb(
STY_(2,mb,1,mb-1),
STY_.diag(0)(1,mb-1), &
C(2,mb,2,mb) );
797 true, UpdateFailedException
798 ,
"Error, the factorization of Q which should be s.p.d. failed with"
799 " the error message: {" << fe.what() <<
"}";
832 x2 = (*x)(mb+1,2*mb);
844 DVectorSlice::const_iterator
845 d_itr =
STY_.diag(0).begin(),
847 DVectorSlice::iterator
849 while( x1_itr != x1.end() )
850 *x1_itr++ = *y2_itr++ / *d_itr++;
867 std::copy( x1a.rbegin(), x1a.rend(), x1b.rbegin() );
903 DVectorSlice::const_iterator
904 d_itr =
STY_.diag(0).begin();
905 DVectorSlice::iterator
907 while( x2_itr != x2.end() )
908 *x2_itr++ /= *d_itr++;
915 throw std::logic_error(
"MatrixSymPosDefLBFGS::assert_initialized() : "
916 "Error, matrix not initialized" );
950 const size_type p = Db_diag.dim();
952 for( size_type i = 1; i <= p; ++i ) {
953 for( size_type j = i; j <= p; ++j ) {
954 value_type &c = (*Cb)(i,j) = 0.0;
955 for( size_type k = 1; k <= i; ++k ) {
956 c += Lb(i,k) * Lb(j,k) / Db_diag(k);
void V_invQtV(DVectorSlice *y, const DVectorSlice &x) const
y = inv(Q) * x
const DMatrixSliceTri Lb() const
Strictly lower triangular part of L.
void init_identity(const VectorSpace &space_diag, value_type alpha)
Teuchos::Ordinal size_type
Typedef for the size type of elements that are used by the library.
AbstractLinAlgPack::size_type size_type
bool maintain_original() const
std::string typeName(const T &t)
void Vp_StV(DVectorSlice *vs_lhs, value_type alpha, const DVectorSlice &vs_rhs)
vs_lhs += alpha * vs_rhs (BLAS xAXPY)
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)
v_lhs *= alpha
const VectorSpace & space_cols() const
Exception for factorization error.
bool original_is_updated_
VectorSpace::space_ptr_t vec_spc_
const DMatrixSliceTriEle tri_ele(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
void update_Q() const
Update Q.
VectorSliceTmpl< value_type > DVectorSlice
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
v_lhs = alpha * V_rhs.
std::ostream & output(std::ostream &o, const COOMatrix &coom)
Output stream function for COOMatrix.
MatrixSymPosDefLBFGS(size_type m=10, bool maintain_original=true, bool maintain_inverse=true, bool auto_rescaling=false)
Calls this->initial_setup()
void V_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = alpha * op(M_rhs1) * V_rhs2.
const DMatrixSliceSym sym(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo)
const DMatrixSliceTri tri(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo, BLAS_Cpp::Diag diag)
void init_diagonal(const Vector &diag)
Actually this calls init_identity( diag.space(), diag.norm_inf() ).
size_type num_secant_updates_
void V_InvMtV(DVector *v_lhs, const DMatrixSliceTri &tri_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2)
void Vp_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs += op(M_rhs1) * V_rhs2.
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = - V_rhs.
DMatrixSliceSym nonconst_sym(DMatrixSlice gms, BLAS_Cpp::Uplo uplo)
Return a symmetric matrix.
void copy(const f_int &N, const f_dbl_prec *X, const f_int &INCX, f_dbl_prec *Y, const f_int &INCY)
const multi_vec_ptr_t S() const
DenseLinAlgPack::DMatrixSliceTriEle DMatrixSliceTriEle
const DMatrixSliceTri R() const
const multi_vec_ptr_t Y() const
void potrf(DMatrixSliceTriEle *A)
Calls xPOTRF to compute the cholesky factorization of a symmetric positive definte matrix...
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)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_dbl_prec const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_dbl_prec LAPACK_C_Decl::f_dbl_prec * C
std::ostream & output(std::ostream &out) const
friend void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
bool BFGS_sTy_suff_p_d(const Vector &s, const Vector &y, const value_type *sTy=NULL, std::ostream *out=NULL, const char func_name[]=NULL)
Check that s'*y is sufficiently positive and print the result if it is not.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
void assert_initialized() const
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
Implementation of limited Memory BFGS matrix for arbitrary vector spaces.
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
AbstractLinAlgPack::value_type value_type
DenseLinAlgPack::DMatrixSliceTri DMatrixSliceTri
void Mp_StM(DMatrixSliceTriEle *tri_lhs, value_type alpha, const DMatrixSliceTriEle &tri_rhs)
tri_lhs += alpha * tri_rhs (BLAS xAXPY)
void Vp_MtV_assert_sizes(size_type v_lhs_size, size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type v_rhs2_size)
v_lhs += op(m_rhs1) * v_rhs2
FortranTypes::f_dbl_prec value_type
Typedef for the value type of elements that is used for the library.
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, const DMatrixSlice &gms_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta=1.0)
vs_lhs = alpha * op(gms_rhs1) * vs_rhs2 + beta * vs_lhs (BLAS xGEMV)
bool maintain_inverse() const
MatrixOp & operator=(const MatrixOp &mwo)
DenseLinAlgPack::DMatrixSlice DMatrixSlice
DenseLinAlgPack::DMatrixSliceSym DMatrixSliceSym
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)
v_lhs += V_rhs.
void Vp_DtV(DVectorSlice *y, const DVectorSlice &x) const
y += D * x