45 #include "SymInvCholMatrixOp.hpp"
52 namespace LinAlgOpPack {
61 namespace ConstrainedOptPack {
74 return MatrixWithOpConcreteEncap<SymInvCholMatrix>::operator=(m);
78 {
return out <<
"\n*** Inverse Cholesky factor:\n" << m().UInv(); }
172 , B.rows(), B.cols(),
trans_not(B_trans) );
174 , B.rows(), B.cols(), B_trans
175 , B.rows(), B.cols(),
trans_not(B_trans) );
205 std::ostringstream omsg;
206 omsg <<
"MatrixSymPosDefInvCholFactor::init_identity(...) : Error, alpha = " << alpha
207 <<
" <= 0.0 and therefore this is not a positive definite matrix.";
208 throw UpdateSkippedException( omsg.str() );
212 m().UInv().diag() = 1.0 /
::sqrt( alpha );
217 DVectorSlice::const_iterator
218 min_ele_ptr = std::min_element( diag.begin(), diag.end() );
219 if( *min_ele_ptr <= 0.0 ) {
220 std::ostringstream omsg;
221 omsg <<
"MatrixSymPosDefInvCholFactor::init_diagonal(...) : Error, "
222 <<
"diag(" << min_ele_ptr - diag.begin() + 1 <<
" ) = "
223 << (*min_ele_ptr) <<
" <= 0.0.\n"
224 <<
"Therefore this is not a positive definite matrix.";
225 throw UpdateSkippedException( omsg.str() );
231 DVectorSlice::const_iterator
232 diag_itr = diag.begin();
233 DVectorSlice::iterator
234 inv_fact_diag_itr = m().UInv().diag().begin();
236 while( diag_itr != diag.end() )
237 *inv_fact_diag_itr++ = 1.0 / ::
sqrt( *diag_itr++ );
247 ConstrainedOptPack::BFGS_update(&m(),s,y,&Bs());
250 ConstrainedOptPack::BFGS_update(&m(),s,y,_Bs);
253 catch(
const BFGSUpdateSkippedException& excpt) {
254 throw UpdateSkippedException( excpt.what() );
void M_StMtM(MatrixOp *M_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2)
M_lhs = alpha * op(M_rhs1) * op(M_rhs2).
AbstractLinAlgPack::size_type size_type
value_type transVtInvMtV(const Vector &v_rhs1, const MatrixNonsing &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
result = v_rhs1' * inv(op(M_rhs2)) * v_rhs3
void assign(DMatrix *gm_lhs, value_type alpha)
gm_lhs = alpha (elementwise)
void MtM_assert_sizes(size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type m_rhs2_rows, size_type m_rhs2_cols, BLAS_Cpp::Transp trans_rhs2)
op(m_lhs) += op(m_rhs1)
void M_StMtInvMtM(DMatrixSliceSym *sym_gms_lhs, value_type alpha, const MatrixOp &mwo, BLAS_Cpp::Transp mwo_trans, EMatrixDummyArg) const
const DMatrixSliceTriEle tri_ele(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo)
void init_diagonal(const DVectorSlice &diag)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return rows of a possible transposed matrix.
MatrixOp & operator=(const MatrixOp &m)
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
void serialize(std::ostream &out) const
void Mp_MtM_assert_sizes(size_type m_lhs_rows, size_type m_lhs_cols, BLAS_Cpp::Transp trans_lhs, size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type m_rhs2_rows, size_type m_rhs2_cols, BLAS_Cpp::Transp trans_rhs2)
op(m_lhs) += op(m_rhs1) * op(m_rhs2)
void secant_update(DVectorSlice *s, DVectorSlice *y, DVectorSlice *Bs)
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
void syrk(BLAS_Cpp::Transp trans, value_type alpha, const DMatrixSlice &gms_rhs, value_type beta, DMatrixSliceSym *sym_lhs)
value_type transVtMtV(const DVectorSlice &vs_rhs1, BLAS_Cpp::Transp trans_rhs2, const DVectorSlice &vs_rhs3) const
const DMatrixSliceTri tri(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo, BLAS_Cpp::Diag diag)
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
result = v_rhs1' * op(M_rhs2) * v_rhs3
std::ostream & output(std::ostream &out) const
value_type transVtInvMtV(const DVectorSlice &vs_rhs1, BLAS_Cpp::Transp trans_rhs2, const DVectorSlice &vs_rhs3) const
void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
Perform a rank-k update of a symmetric matrix of the form:
DenseLinAlgPack::DMatrixSliceTriEle DMatrixSliceTriEle
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
v_lhs = inv(op(M_rhs1)) * v_rhs2
f_dbl_prec f_dbl_prec f_dbl_prec * S
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)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
void extract_inv_chol(DMatrixSliceTriEle *InvChol) const
const f_int f_dbl_prec a[]
void init_identity(size_type n, value_type alpha)
void sqrt(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs)
vs_lhs = sqrt(vs_rhs)
SparseVectorSlice< SparseElement< index_type, value_type > > SpVectorSlice
Transp trans_not(Transp _trans)
Return the opposite of the transpose argument.
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
void V_InvMtV(DVector *v_lhs, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2) const
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
void unserialize(std::istream &in)
void M_StInvMtM(MatrixOp *m_lhs, value_type alpha, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2)
m_lhs = alpha * inv(op(mwo_rhs1)) * op(mwo_rhs2) (right)
DenseLinAlgPack::DMatrixSliceSym DMatrixSliceSym
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void M_StInvMtM(DMatrix *gm_lhs, value_type alpha, const DMatrixSliceTri &tri_rhs1, BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice &gms_rhs2, BLAS_Cpp::Transp trans_rhs2)