73 #ifdef HAVE_MOOCHO_FORTRAN
74 # define ALAP_DCHUD_DECL FORTRAN_FUNC_DECL_UL( void, DCHUD, dchud )
75 # define ALAP_DCHUD_CALL FORTRAN_FUNC_CALL_UL( DCHUD, dchud )
77 # define ALAP_DCHUD_DECL void dchud_c
78 # define ALAP_DCHUD_CALL dchud_c
95 T my_max(
const T& v1,
const T& v2 ) {
return v1 > v2 ? v1 : v2; }
99 T my_min(
const T& v1,
const T& v2 ) {
return v1 < v2 ? v1 : v2; }
102 namespace AbstractLinAlgPack {
106 MatrixSymPosDefCholFactor::MatrixSymPosDefCholFactor()
107 : maintain_original_(true), maintain_factor_(false)
108 , factor_is_updated_(false), allocates_storage_(true)
109 , max_size_(0), M_size_(0), M_l_r_(1), M_l_c_(1)
110 , U_l_r_(1), U_l_c_(1), scale_(0.0), is_diagonal_(false)
111 , pivot_tols_(0.0,0.0,0.0)
114 MatrixSymPosDefCholFactor::MatrixSymPosDefCholFactor(
116 ,
const release_resource_ptr_t& release_resource_ptr
118 ,
bool maintain_original
119 ,
bool maintain_factor
124 : pivot_tols_(0.0,0.0,0.0)
126 init_setup(MU_store,release_resource_ptr,max_size,maintain_original
127 ,maintain_factor,allow_factor,set_full_view,scale);
130 void MatrixSymPosDefCholFactor::init_setup(
132 ,
const release_resource_ptr_t& release_resource_ptr
134 ,
bool maintain_original
135 ,
bool maintain_factor
142 if( MU_store == NULL ) {
143 maintain_original_ = maintain_original;
144 maintain_factor_ = maintain_factor;
145 factor_is_updated_ = maintain_factor;
146 allocates_storage_ =
true;
149 max_size_ = max_size;
152 if( !maintain_factor && !allow_factor )
157 is_diagonal_ =
false;
161 allocates_storage_ =
false;
162 MU_store_.bind(*MU_store);
163 release_resource_ptr_ = release_resource_ptr;
164 max_size_ = my_min( MU_store->rows(), MU_store->cols() ) - 1;
165 if( set_full_view ) {
169 ,scale,maintain_original,1,1
170 ,maintain_factor, allow_factor ? 1: 0, allow_factor ? 1 : 0
176 ,0.0,maintain_original,1,1
177 ,maintain_factor, allow_factor ? 1: 0, allow_factor ? 1 : 0
183 void MatrixSymPosDefCholFactor::set_view(
186 ,
bool maintain_original
189 ,
bool maintain_factor
196 allocate_storage(max_size_);
198 allocate_storage( my_max( M_l_r + M_size, M_l_c + M_size ) - 1 );
200 if( maintain_original ) {
204 if( maintain_factor ) {
209 maintain_original_ = maintain_original;
210 maintain_factor_ = maintain_factor;
211 is_diagonal_ =
false;
213 max_size_ = my_min( MU_store_.rows() - U_l_r, MU_store_.cols() - U_l_c );
220 factor_is_updated_ = maintain_factor;
230 factor_is_updated_ = maintain_factor;
234 void MatrixSymPosDefCholFactor::pivot_tols( PivotTolerances pivot_tols )
236 pivot_tols_ = pivot_tols;
239 MatrixSymAddDelUpdateable::PivotTolerances MatrixSymPosDefCholFactor::pivot_tols()
const
253 void MatrixSymPosDefCholFactor::zero_out()
255 this->init_identity(this->space_cols(),0.0);
263 if( maintain_original_ ) {
265 <<
"Unfactored symmetric matrix stored as lower triangle (ignore upper nonzeros):\n"
268 if( factor_is_updated_ ) {
270 <<
"Matrix scaling M = scale*U'*U, scale = " << scale_ << std::endl
271 <<
"Upper cholesky factor U (ignore lower nonzeros):\n"
286 MatrixSymOpGetGMSSymMutable
287 *symwo_gms_lhs =
dynamic_cast<MatrixSymOpGetGMSSymMutable*
>(m_lhs);
290 MatrixDenseSymMutableEncap sym_lhs(symwo_gms_lhs);
306 !maintain_original_, std::logic_error
307 ,
"MatrixSymPosDefCholFactor::Mp_StM(alpha,M_rhs,trans_rhs): Error, Current implementation "
308 "can not perform this operation unless the original matrix is being maintained." );
311 bool diag_op =
false;
312 if(
const MatrixSymOpGetGMSSym *symwo_gms_rhs = dynamic_cast<const MatrixSymOpGetGMSSym*>(&M_rhs)) {
314 MatrixDenseSymEncap sym_rhs(*symwo_gms_rhs);
323 else if(
const MatrixSymDiag *symwo_diag_rhs = dynamic_cast<const MatrixSymDiag*>(&M_rhs)) {
325 VectorDenseEncap sym_rhs_diag(symwo_diag_rhs->diag());
330 else if(
const MatrixSymOp *symwo_rhs = dynamic_cast<const MatrixSymOp*>(&M_rhs)) {
336 if( diag_op && is_diagonal_ )
337 this->init_diagonal(VectorMutableDense(this->
M().gms().diag(),
Teuchos::null));
339 this->initialize(this->
M());
346 const MatrixOp &mwo_rhs
352 MatrixDenseSymMutableEncap sym_gms_lhs(
this);
353 const MatrixOpGetGMS *mwo_rhs_gms =
dynamic_cast<const MatrixOpGetGMS*
>(&mwo_rhs);
363 const VectorSpace &spc = mwo_rhs.space_rows();
365 VectorSpace::multi_vec_mut_ptr_t
S = spc.create_members(m);
371 LinAlgOpPack::Vp_StV( &sym_gms_lhs().gms().
col(j)(j,m), alpha, VectorDenseEncap(*S->col(j)->sub_view(j,m))() );
376 LinAlgOpPack::Vp_StV( &sym_gms_lhs().gms().
col(j)(1,j), alpha, VectorDenseEncap(*S->col(j)->sub_view(1,j))() );
394 #ifdef PROFILE_HACK_ENABLED
397 assert_initialized();
402 if( maintain_original_ ) {
428 if( a != 1.0 || scale_ != 1.0 )
454 #ifdef PROFILE_HACK_ENABLED
457 assert_initialized();
458 if( maintain_original_ ) {
471 for( SpVectorSlice::const_iterator x_itr = x.begin(); x_itr != x.end(); ++x_itr ) {
472 const size_type i = x_itr->index() + x.offset();
488 #ifdef PROFILE_HACK_ENABLED
491 assert_initialized();
500 #ifdef PROFILE_HACK_ENABLED
503 assert_initialized();
504 if( maintain_original_ ) {
531 , EMatRhsPlaceHolder dummy_place_holder
537 #ifdef PROFILE_HACK_ENABLED
540 assert_initialized();
541 if( !maintain_original_ ) {
557 #ifdef PROFILE_HACK_ENABLED
560 assert_initialized();
573 update_factorization();
590 #ifdef PROFILE_HACK_ENABLED
593 assert_initialized();
604 #ifdef PROFILE_HACK_ENABLED
621 assert_initialized();
622 update_factorization();
645 ,B.rows(), B.cols(), B_trans
663 M, M.rows(), maintain_factor_
664 , MatrixSymAddDelUpdateable::Inertia()
665 , MatrixSymAddDelUpdateable::PivotTolerances()
674 !maintain_original_, std::logic_error
675 ,
"MatrixSymPosDefCholFactor::get_sym_gms_view(): Error, maintain_original must be "
676 "true in order to call this method!" );
683 !maintain_original_, std::logic_error
684 ,
"MatrixSymPosDefCholFactor::free_sym_gms_view(...): Error, maintain_original must be "
685 "true in order to call this method!" );
694 !maintain_original_, std::logic_error
695 ,
"MatrixSymPosDefCholFactor::get_sym_gms_view(): Error, maintain_original must be "
696 "true in order to call this method!" );
703 !maintain_original_, std::logic_error
704 ,
"MatrixSymPosDefCholFactor::commit_sym_gms_view(...): Error, maintain_original must be "
705 "true in order to call this method!" );
706 this->initialize(*sym_gms_view);
713 assert_initialized();
714 update_factorization();
735 scale_ < 0.0, std::logic_error
736 ,
"MatrixSymPosDefCholFactor::extract_inv_chol(...) : "
737 "Error, we can not compute the inverse cholesky factor "
738 "af a negative definite matrix." );
740 InvChol->gms().diag() = 1.0;
742 &InvChol->gms(), 1.0 /
std::sqrt(scale_), U()
749 void MatrixSymPosDefCholFactor::init_identity(
const VectorSpace& space_diag,
value_type alpha )
752 allocate_storage( max_size_ ? max_size_ : n );
760 const value_type scale = alpha > 0.0 ? +1.0: -1.0;
761 resize_and_zero_off_diagonal(n,scale);
762 if( maintain_original_ ) {
763 M().gms().diag()(1,n) = alpha;
765 if( maintain_factor_ ) {
766 U().gms().diag()(1,n) =
std::sqrt(std::fabs(alpha));
767 factor_is_updated_ =
true;
772 void MatrixSymPosDefCholFactor::init_diagonal(
const Vector& diag_in )
774 VectorDenseEncap diag_encap(diag_in);
777 allocate_storage( max_size_ ? max_size_ : diag.dim() );
783 if( diag.dim() == 0 ) {
787 const value_type scale = diag(1) > 0.0 ? +1.0: -1.0;
788 resize_and_zero_off_diagonal(diag.dim(),scale);
789 if( maintain_original_ ) {
790 M().gms().diag() = diag;
793 if( maintain_factor_ ) {
800 factor_is_updated_ =
true;
805 void MatrixSymPosDefCholFactor::secant_update(
806 VectorMutable* s_in, VectorMutable* y_in, VectorMutable* Bs_in
814 namespace rcp = MemMngPack;
816 assert_initialized();
823 VectorDenseMutableEncap s_encap(*s_in);
824 VectorDenseMutableEncap y_encap(*y_in);
825 VectorDenseMutableEncap Bs_encap( Bs_in ? *Bs_in : *y_in);
829 *Bs = ( Bs_in ? &Bs_encap() : NULL );
834 sTy_scale = sTy*scale_;
835 std::ostringstream omsg;
839 ,&sTy_scale,&omsg,
"\nMatrixSymPosDefCholFactor::secant_update(...)"
843 throw UpdateSkippedException( omsg.str() );
850 Bs_view.bind( Bs_store() );
856 scale_*sTBs <= 0.0 && scale_ > 0.0, UpdateFailedException
857 ,
"MatrixSymPosDefCholFactor::secant_update(...) : "
858 "Error, B can't be positive definite if s'*Bs <= 0.0" );
860 scale_*sTBs <= 0.0 && scale_ <= 0.0, UpdateFailedException
861 ,
"MatrixSymPosDefCholFactor::secant_update(...) : "
862 "Error, B can't be negative definite if s'*Bs >= 0.0" );
863 if( maintain_original_ ) {
875 if( maintain_factor_ ) {
911 a, &v, u, v.dim() > 1 ? &U.gms().diag(-1) : NULL
915 factor_is_updated_ =
false;
917 is_diagonal_ =
false;
922 void MatrixSymPosDefCholFactor::initialize(
927 #ifdef PROFILE_HACK_ENABLED
930 allocate_storage(max_size);
933 throw SingularUpdateException(
934 "MatrixSymPosDefCholFactor::initialize(...): "
935 "Error, alpha == 0.0, matrix is singular.", 0.0 );
938 if( maintain_original_ ) {
947 max_size_ = my_min( MU_store_.rows(), MU_store_.cols() ) - 1;
948 scale_ = alpha > 0.0 ? +1.0 : -1.0;
950 if( maintain_original_ ) {
951 M().gms()(1,1) = alpha;
954 U().gms()(1,1) =
std::sqrt( scale_ * alpha );
955 factor_is_updated_ =
true;
957 is_diagonal_ =
false;
960 void MatrixSymPosDefCholFactor::initialize(
963 ,
bool force_factorization
965 ,PivotTolerances pivot_tols
968 #ifdef PROFILE_HACK_ENABLED
971 typedef MatrixSymAddDelUpdateable::Inertia Inertia;
973 allocate_storage(max_size);
980 TEUCHOS_TEST_FOR_EXCEPT( !( (inertia.neg_eigens == Inertia::UNKNOWN && inertia.pos_eigens == Inertia::UNKNOWN ) )
981 || ( inertia.neg_eigens == n && inertia.pos_eigens == 0 )
982 || ( inertia.neg_eigens == 0 && inertia.pos_eigens == n )
990 const int sign_a_11 = ( a_11 == 0.0 ? 0 : ( a_11 > 0.0 ? +1 : -1 ) );
991 if( sign_a_11 == 0.0 )
993 "MatrixSymPosDefCholFactor::initialize(...) : "
994 "Error, A can not be positive definite or negative definete "
995 "if A(1,1) == 0.0" );
996 if( inertia.pos_eigens == n && sign_a_11 < 0 )
998 "MatrixSymPosDefCholFactor::initialize(...) : "
999 "Error, A can not be positive definite "
1000 "if A(1,1) < 0.0" );
1001 if( inertia.neg_eigens == n && sign_a_11 > 0 )
1003 "MatrixSymPosDefCholFactor::initialize(...) : "
1004 "Error, A can not be negative definite "
1005 "if A(1,1) > 0.0" );
1011 n,scale,maintain_original_,1,1
1012 ,maintain_factor_, U_l_r_ ? 1 : 0, U_l_r_ ? 1 : 0
1015 if( maintain_original_ ) {
1022 if( maintain_factor_ || force_factorization ) {
1030 factor_is_updated_ =
true;
1034 throw SingularUpdateException( excpt.what(), 0.0 );
1046 DVectorSlice::iterator
1047 U_itr = U_ele.gms().diag().begin(),
1048 U_end = U_ele.gms().diag().end();
1049 while( U_itr != U_end ) {
1051 if(U_abs < min_diag) min_diag = U_abs;
1052 if(U_abs > max_diag) max_diag = U_abs;
1054 const value_type gamma = (min_diag*min_diag)/(max_diag*max_diag);
1056 PivotTolerances use_pivot_tols = pivot_tols_;
1057 if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
1058 use_pivot_tols.warning_tol = pivot_tols.warning_tol;
1059 if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
1060 use_pivot_tols.singular_tol = pivot_tols.singular_tol;
1061 if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
1062 use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
1063 const bool throw_exception = (
1065 || gamma <= use_pivot_tols.warning_tol
1066 || gamma <= use_pivot_tols.singular_tol
1069 std::ostringstream onum_msg, omsg;
1070 if(throw_exception) {
1072 <<
"gamma = (min(|diag(U)(k)|)/|max(|diag(U)(k)|))^2 = (" << min_diag <<
"/"
1073 << max_diag <<
")^2 = " << gamma;
1075 <<
"MatrixSymPosDefCholFactor::initialize(...): ";
1076 if( gamma <= use_pivot_tols.singular_tol ) {
1079 <<
"Singular update!\n" << onum_msg.str() <<
" <= singular_tol = "
1080 << use_pivot_tols.singular_tol;
1081 throw SingularUpdateException( omsg.str(), gamma );
1083 else if( gamma <= use_pivot_tols.warning_tol ) {
1085 <<
"Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
1086 <<
" < " << onum_msg.str() <<
" <= warning_tol = "
1087 << use_pivot_tols.warning_tol;
1088 throw WarnNearSingularUpdateException( omsg.str(), gamma );
1096 factor_is_updated_ =
false;
1100 size_type MatrixSymPosDefCholFactor::max_size()
const
1105 MatrixSymAddDelUpdateable::Inertia
1106 MatrixSymPosDefCholFactor::inertia()
const
1108 typedef MatrixSymAddDelUpdateable MSADU;
1109 typedef MSADU::Inertia Inertia;
1112 ? Inertia(0,0,M_size_)
1113 : Inertia(M_size_,0,0) )
1117 void MatrixSymPosDefCholFactor::set_uninitialized()
1122 void MatrixSymPosDefCholFactor::augment_update(
1125 ,
bool force_refactorization
1126 ,EEigenValType add_eigen_val
1127 ,PivotTolerances pivot_tols
1130 #ifdef PROFILE_HACK_ENABLED
1136 typedef MatrixSymAddDelUpdateable MSADU;
1138 assert_initialized();
1142 implicit_ref_cast<const MatrixBase>(*this).rows() >= max_size(),
1143 MaxSizeExceededException,
1144 "MatrixSymPosDefCholFactor::augment_update(...) : "
1145 "Error, the maximum size would be exceeded." );
1147 t && t->dim() != M_size_, std::length_error
1148 ,
"MatrixSymPosDefCholFactor::augment_update(...): "
1149 "Error, t.dim() must be equal to this->rows()." );
1150 if( !(add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN || add_eigen_val != MSADU::EIGEN_VAL_ZERO ) )
1151 throw SingularUpdateException(
1152 "MatrixSymPosDefCholFactor::augment_update(...): "
1153 "Error, the client has specified a singular update in add_eigen_val.", -1.0 );
1155 throw SingularUpdateException(
1156 "MatrixSymPosDefCholFactor::augment_update(...): "
1157 "Error, alpha == 0.0 so the matrix is not positive definite "
1158 "or negative definite.", -1.0 );
1159 if( scale_ > 0.0 && alpha < 0.0 )
1160 throw WrongInertiaUpdateException(
1161 "MatrixSymPosDefCholFactor::augment_update(...): "
1162 "Error, alpha < 0.0 so the matrix is not postivie definite ", -1.0 );
1163 if( scale_ < 0.0 && alpha > 0.0 )
1164 throw WrongInertiaUpdateException(
1165 "MatrixSymPosDefCholFactor::augment_update(...): "
1166 "Error, alpha > 0.0 so the matrix is not negative definite ", -1.0 );
1167 if( scale_ > 0.0 && !( add_eigen_val == MSADU::EIGEN_VAL_POS || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
1168 throw WrongInertiaUpdateException(
1169 "MatrixSymPosDefCholFactor::augment_update(...): "
1170 "Error, alpha > 0.0 but the user specified a non-positive definite new matrix.", -1.0 );
1171 if( scale_ < 0.0 && !( add_eigen_val == MSADU::EIGEN_VAL_NEG || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
1172 throw WrongInertiaUpdateException(
1173 "MatrixSymPosDefCholFactor::augment_update(...): "
1174 "Error, alpha > 0.0 but the user specified a non-positive definite new matrix.", -1.0 );
1176 bool throw_exception =
false;
1177 std::ostringstream omsg;
1179 if( maintain_factor_ ) {
1212 DVectorSlice u12 = MU_store_.col(U_l_c_+M_size_+1)(U_l_r_,U_l_r_+M_size_-1);
1213 value_type &u22 = MU_store_(U_l_r_+M_size_,U_l_c_+M_size_+1);
1224 u22sqr = (1/scale_) * alpha - ( t ?
dot( u12, u12 ) : 0.0 ),
1226 nrm_U_diag =
norm_inf(U().gms().diag()),
1227 sqr_nrm_U_diag = nrm_U_diag * nrm_U_diag;
1229 gamma = u22sqrabs / sqr_nrm_U_diag;
1232 correct_inertia = u22sqr > 0.0;
1233 PivotTolerances use_pivot_tols = pivot_tols_;
1234 if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
1235 use_pivot_tols.warning_tol = pivot_tols.warning_tol;
1236 if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
1237 use_pivot_tols.singular_tol = pivot_tols.singular_tol;
1238 if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
1239 use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
1242 || correct_inertia && gamma <= use_pivot_tols.warning_tol
1243 || correct_inertia && gamma <= use_pivot_tols.singular_tol
1247 std::ostringstream onum_msg;
1248 if(throw_exception) {
1250 <<
"gamma = |u22^2|/(||diag(U11)||inf)^2 = |" << u22sqr <<
"|/("
1251 << nrm_U_diag <<
")^2 = " << gamma;
1253 <<
"MatrixSymPosDefCholFactor::augment_update(...): ";
1254 if( correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
1256 <<
"Singular update!\n" << onum_msg.str() <<
" <= singular_tol = "
1257 << use_pivot_tols.singular_tol;
1258 throw SingularUpdateException( omsg.str(), gamma );
1260 else if( !correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
1262 <<
"Singular update!\n" << onum_msg.str() <<
" <= wrong_inertia_tol = "
1263 << use_pivot_tols.wrong_inertia_tol;
1264 throw SingularUpdateException( omsg.str(), gamma );
1267 else if( !correct_inertia ) {
1269 <<
"Indefinite update!\n" << onum_msg.str() <<
" >= wrong_inertia_tol = "
1270 << use_pivot_tols.wrong_inertia_tol;
1271 throw WrongInertiaUpdateException( omsg.str(), gamma );
1273 else if( correct_inertia && gamma <= use_pivot_tols.warning_tol ) {
1275 <<
"Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
1276 <<
" < " << onum_msg.str() <<
" <= warning_tol = "
1277 << use_pivot_tols.warning_tol;
1288 factor_is_updated_ =
false;
1291 if( maintain_original_ ) {
1296 DVectorSlice M12 = MU_store_.row(M_l_r_+M_size_+1)(M_l_c_,M_l_c_+M_size_-1);
1301 MU_store_(M_l_r_+M_size_+1,M_l_c_+M_size_) = alpha;
1304 is_diagonal_ =
false;
1305 if( throw_exception )
1306 throw WarnNearSingularUpdateException(omsg.str(),gamma);
1309 void MatrixSymPosDefCholFactor::delete_update(
1311 ,
bool force_refactorization
1312 ,EEigenValType drop_eigen_val
1313 ,PivotTolerances pivot_tols
1316 #ifdef PROFILE_HACK_ENABLED
1319 typedef MatrixSymAddDelUpdateable MSADU;
1322 jd < 1 || M_size_ < jd, std::out_of_range
1323 ,
"MatrixSymPosDefCholFactor::delete_update(jd,...): "
1324 "Error, the indice jd must be 1 <= jd <= rows()" );
1327 || (scale_ > 0.0 && drop_eigen_val == MSADU::EIGEN_VAL_POS)
1328 || (scale_ < 0.0 && drop_eigen_val == MSADU::EIGEN_VAL_NEG)
1331 if( maintain_original_ ) {
1353 if( maintain_factor_ ) {
1386 if(work_.dim() < work_size)
1387 work_.resize(work_size);
1397 c = work_(rng+_size),
1398 s = work_(rng+2*_size);
1401 u23 = U.row(jd)(U_rng);
1405 U33.col_ptr(1), U33.max_rows()
1406 ,U_rng.size(), u23.start_ptr(), &dummy, 1, 0, &dummy
1407 ,&dummy, c.start_ptr(), s.start_ptr() );
1414 factor_is_updated_ =
false;
1426 void MatrixSymPosDefCholFactor::serialize( std::ostream &out )
const
1429 out << build_serialization_string() << std::endl;
1431 out.precision(std::numeric_limits<value_type>::digits10+4);
1433 out << M_size_ << std::endl;
1436 if( maintain_original_ ) {
1438 write_matrix( M.gms(), M.uplo(),
out );
1442 write_matrix( U.gms(), U.uplo(),
out );
1448 void MatrixSymPosDefCholFactor::unserialize( std::istream &in )
1451 std::string keywords;
1452 std::getline( in, keywords,
'\n' );
1454 const std::string this_keywords = build_serialization_string();
1456 this_keywords != keywords, std::logic_error
1457 ,
"MatrixSymPosDefCholFactor::unserialize(...): Error, the matrix type being read in from file of "
1458 "\'"<<keywords<<
"\' does not equal the type expected of \'"<<this_keywords<<
"\'!"
1463 M_size_ < 0, std::logic_error
1464 ,
"MatrixSymPosDefCholFactor::unserialize(...): Error, read in a size of M_size = "<<M_size_<<
" < 0!"
1466 allocate_storage(M_size_);
1468 if(maintain_original_) {
1469 M_l_r_ = M_l_c_ = 1;
1474 U_l_r_ = U_l_c_ = 1;
1482 void MatrixSymPosDefCholFactor::assert_storage()
const
1487 void MatrixSymPosDefCholFactor::allocate_storage(
size_type max_size)
const
1489 namespace rcp = MemMngPack;
1490 if( allocates_storage_ && MU_store_.rows() < max_size + 1 ) {
1495 const_cast<MatrixSymPosDefCholFactor*
>(
this)->release_resource_ptr_ =
Teuchos::rcp(
new ptr_t(MU_store));
1496 const_cast<MatrixSymPosDefCholFactor*
>(
this)->MU_store_.bind( (*MU_store)() );
1497 const_cast<MatrixSymPosDefCholFactor*
>(
this)->max_size_ = max_size;
1504 void MatrixSymPosDefCholFactor::assert_initialized()
const
1509 void MatrixSymPosDefCholFactor::resize_and_zero_off_diagonal(
size_type n,
value_type scale)
1518 n, scale, maintain_original_, 1, 1, maintain_factor_
1519 ,U_l_r_ ? 1 : 0, U_l_r_ ? 1 : 0 );
1520 if( maintain_original_ ) {
1521 if(!is_diagonal_ && n > 1)
1525 if(!is_diagonal_ && n > 1)
1530 void MatrixSymPosDefCholFactor::update_factorization()
const
1532 if( factor_is_updated_ )
return;
1534 U_l_r_ == 0, std::logic_error
1535 ,
"MatrixSymPosDefCholFactor::update_factorization() : "
1536 "Error, U_l_r == 0 was set in MatrixSymPosDefCholFactor::set_view(...) "
1537 "and therefore we can not update the factorization in the provided storage space." );
1539 #ifdef PROFILE_HACK_ENABLED
1542 MatrixSymPosDefCholFactor
1543 *nc_this =
const_cast<MatrixSymPosDefCholFactor*
>(
this);
1547 #ifdef PROFILE_HACK_ENABLED
1552 nc_this->factor_is_updated_ =
true;
1555 std::string MatrixSymPosDefCholFactor::build_serialization_string()
const
1557 std::string str =
"SYMMETRIC POS_DEF";
1558 if( !maintain_original_ )
1559 str.append(
" CHOL_FACTOR");
1560 if( maintain_original_ )
1561 str.append(
" LOWER");
1563 str.append(
" UPPER");
1569 const int Q_dim = Q.rows();
1571 for(
int i = 1; i <= Q_dim; ++i ) {
1572 for(
int j = 1; j <= i; ++j ) {
1573 out <<
" " << Q(i,j);
1579 for(
int i = 1; i <= Q_dim; ++i ) {
1580 for(
int j = i; j <= Q_dim; ++j ) {
1581 out <<
" " << Q(i,j);
1591 const int Q_dim = Q.
rows();
1593 for(
int i = 1; i <= Q_dim; ++i ) {
1594 for(
int j = 1; j <= i; ++j ) {
1595 #ifdef TEUCHOS_DEBUG
1597 in.eof(), std::logic_error
1598 ,
"MatrixSymPosDefCholFactor::read_matrix(in,lower,Q_out): Error, not finished reading in matrix yet (i="<<i<<
",j="<<j<<
")!"
1606 for(
int i = 1; i <= Q_dim; ++i ) {
1607 for(
int j = i; j <= Q_dim; ++j ) {
1608 #ifdef TEUCHOS_DEBUG
1610 in.eof(), std::logic_error
1611 ,
"MatrixSymPosDefCholFactor::read_matrix(in,upper,Q_out): Error, not finished reading in matrix yet (i="<<i<<
",j="<<j<<
")!"
f_dbl_prec const f_int f_dbl_prec * Y
void V_MtV(DVector *v_lhs, const DMatrixSliceTri &tri_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2)
AbstractLinAlgPack::size_type size_type
void Vp_StV(DVectorSlice *vs_lhs, value_type alpha, const DVectorSlice &vs_rhs)
vs_lhs += alpha * vs_rhs (BLAS xAXPY)
RTOp_index_type index_type
FortranTypes::f_int f_int
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)
const f_int const f_int const f_dbl_prec const f_dbl_prec * X
Exception for factorization error.
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs2, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta=1.0)
v_lhs = alpha * op(P_rhs1) * op(M_rhs2) * v_rhs3 + beta * v_rhs
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 abs(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs)
vs_lhs = abs(vs_rhs)
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return rows of a possible transposed matrix.
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
std::ostream & output(std::ostream &o, const COOMatrix &coom)
Output stream function for COOMatrix.
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)
Helper class that takes care of timing.
void M_MtM(MatrixOp *M_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2)
M_lhs = op(M_rhs1) * op(M_rhs2).
bool assert_print_nan_inf(const value_type &val, const std::string &name, bool throw_excpt, std::ostream *out)
TypeTo & implicit_ref_cast(TypeFrom &t)
void syrk(BLAS_Cpp::Transp trans, value_type alpha, const DMatrixSlice &gms_rhs, value_type beta, DMatrixSliceSym *sym_lhs)
const DMatrixSliceTri tri(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo, BLAS_Cpp::Diag diag)
void V_InvMtV(DVector *v_lhs, const DMatrixSliceTri &tri_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2)
void Mp_StPtMtP(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs, const GenPermMatrixSlice &P_rhs2, BLAS_Cpp::Transp P_rhs2_trans)
mwo_lhs += alpha * op(P_rhs1) * op(M_rhs) * op(P_rhs2).
value_type norm_2(const DVectorSlice &vs_rhs)
result = ||vs_rhs||2 (BLAS xNRM2)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void syr(value_type alpha, const DVectorSlice &vs_rhs, DMatrixSliceSym *sym_lhs)
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:
void Mt_S(DMatrixSlice *gms_lhs, value_type alpha)
gms_lhs *= alpha (BLAS xSCAL)
DenseLinAlgPack::DMatrixSliceTriEle DMatrixSliceTriEle
void rank_2_chol_update(const value_type a, DVectorSlice *u, const DVectorSlice &v, DVectorSlice *w, DMatrixSliceTriEle *R, BLAS_Cpp::Transp R_trans)
Update the cholesky factor of symmetric positive definite matrix for a rank-2 update.
const LAPACK_C_Decl::f_int & M
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 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)
DVectorSlice col(DMatrixSlice &gms, BLAS_Cpp::Transp trans, size_type j)
value_type norm_inf(const SparseVectorSlice< T_Ele > &sv_rhs)
result = ||sv_rhs||inf (BLAS IxAMAX)
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.
DMatrixSliceTriEle nonconst_tri_ele(DMatrixSlice gms, BLAS_Cpp::Uplo uplo)
Return a triangular element-wise matrix.
const f_int f_dbl_prec a[]
Template class that implements ReleaseResource interface for a RCP<T> object.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
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
size_type rows() const
Return the number of rows.
void Vt_S(DVectorSlice *vs_lhs, value_type alpha)
vs_lhs *= alpha (BLAS xSCAL) (*** Note that alpha == 0.0 is handeled as vs_lhs = 0.0)
int read_matrix(const std::string &filename, Epetra_Comm &Comm, const Epetra_Map *rowmap, Epetra_Map *colmap, const Epetra_Map *rangemap, const Epetra_Map *domainmap, Epetra_CrsMatrix *&mat)
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
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)
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)
value_type norm_inf(const DVectorSlice &vs_rhs)
result = ||vs_rhs||infinity (BLAS IxAMAX)
void delete_row_col(size_type kd, DMatrixSliceTriEle *M)
DenseLinAlgPack::DMatrixSlice DMatrixSlice
void M_StMtInvMtM(MatrixSymOp *sym_gms_lhs, value_type alpha, const MatrixOp &mwo, BLAS_Cpp::Transp mwo_trans, const MatrixSymNonsing &mswof, MatrixSymNonsing::EMatrixDummyArg mwo_rhs)
sym_gms_lhs = alpha * op(mwo) * inv(mswof) * op(mwo)'
DenseLinAlgPack::DMatrixSliceSym DMatrixSliceSym
FortranTypes::f_dbl_prec f_dbl_prec
DenseLinAlgPack::DMatrix DMatrix
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return columns of a possible transposed matrix.
RangePack::Range1D Range1D
const f_int f_dbl_prec const f_int f_int const f_int f_int const f_dbl_prec & u
#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)
value_type dot(const DVectorSlice &vs_rhs1, const DVectorSlice &vs_rhs2)
result = vs_rhs1' * vs_rhs2 (BLAS xDOT)