AbstractLinAlgPack_MatrixSymPosDefCholFactor.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <assert.h>
00030 
00031 #include <limits>
00032 #include <vector>
00033 
00034 #include "AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp"
00035 #include "AbstractLinAlgPack_BFGS_helpers.hpp"
00036 #include "AbstractLinAlgPack_rank_2_chol_update.hpp"
00037 #include "AbstractLinAlgPack_VectorMutableDense.hpp"
00038 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00039 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00040 #include "AbstractLinAlgPack_MatrixOpGetGMS.hpp"
00041 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00042 #include "AbstractLinAlgPack_VectorSpace.hpp"
00043 #include "AbstractLinAlgPack_VectorMutable.hpp"
00044 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00045 #include "AbstractLinAlgPack_MatrixSymDiag.hpp"
00046 #include "AbstractLinAlgPack_VectorSpaceFactory.hpp"
00047 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00048 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00049 #include "DenseLinAlgLAPack.hpp"
00050 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00051 #include "DenseLinAlgPack_DMatrixOp.hpp"
00052 #include "DenseLinAlgPack_DMatrixOut.hpp"
00053 #include "DenseLinAlgPack_delete_row_col.hpp"
00054 #include "DenseLinAlgPack_assert_print_nan_inf.hpp"
00055 #include "ReleaseResource_ref_count_ptr.hpp"
00056 #include "ProfileHackPack_profile_hack.hpp"
00057 #include "Teuchos_TestForException.hpp"
00058 #include "Teuchos_FancyOStream.hpp"
00059 
00060 // Helper functions
00061 extern "C" {
00062   FORTRAN_FUNC_DECL_UL( void, DCHUD, dchud ) ( FortranTypes::f_dbl_prec* R
00063     , const FortranTypes::f_int& LDR, const FortranTypes::f_int& P
00064     , FortranTypes::f_dbl_prec* X, FortranTypes::f_dbl_prec* Z
00065     , const FortranTypes::f_int& LDZ, const FortranTypes::f_int& NZ
00066     , FortranTypes::f_dbl_prec* Y, FortranTypes::f_dbl_prec*  RHO
00067     , FortranTypes::f_dbl_prec* C, FortranTypes::f_dbl_prec* S );
00068 } // end extern "C"
00069 
00070 namespace {
00071 //
00072 template< class T >
00073 inline
00074 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00075 //
00076 template< class T >
00077 inline
00078 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00079 } // end namespace
00080 
00081 namespace AbstractLinAlgPack {
00082 
00083 // Constructors/initalizers
00084 
00085 MatrixSymPosDefCholFactor::MatrixSymPosDefCholFactor()
00086   : maintain_original_(true), maintain_factor_(false)
00087   , factor_is_updated_(false), allocates_storage_(true)
00088   , max_size_(0), M_size_(0), M_l_r_(1), M_l_c_(1)
00089   , U_l_r_(1), U_l_c_(1), scale_(0.0), is_diagonal_(false)
00090   , pivot_tols_(0.0,0.0,0.0) // This is what DPOTRF(...) uses!
00091 {}
00092 
00093 MatrixSymPosDefCholFactor::MatrixSymPosDefCholFactor(
00094   DMatrixSlice                      *MU_store
00095   ,const release_resource_ptr_t&    release_resource_ptr
00096   ,size_type                        max_size
00097   ,bool                             maintain_original
00098   ,bool                             maintain_factor
00099   ,bool                             allow_factor
00100   ,bool                             set_full_view
00101   ,value_type                       scale
00102   )
00103   : pivot_tols_(0.0,0.0,0.0) // This is what DPOTRF(...) uses!
00104 {
00105   init_setup(MU_store,release_resource_ptr,max_size,maintain_original
00106          ,maintain_factor,allow_factor,set_full_view,scale);
00107 }
00108 
00109 void MatrixSymPosDefCholFactor::init_setup(
00110   DMatrixSlice                      *MU_store
00111   ,const release_resource_ptr_t&    release_resource_ptr
00112   ,size_type                        max_size
00113   ,bool                             maintain_original
00114   ,bool                             maintain_factor
00115   ,bool                             allow_factor
00116   ,bool                             set_full_view
00117   ,value_type                       scale
00118   )
00119 {
00120   TEST_FOR_EXCEPT( !(  maintain_original || maintain_factor  ) );
00121   if( MU_store == NULL ) {
00122     maintain_original_ = maintain_original;
00123     maintain_factor_   = maintain_factor;
00124     factor_is_updated_ = maintain_factor;
00125     allocates_storage_ = true; // We will be able to allocate our own storage!
00126     release_resource_ptr_ = Teuchos::null; // Free any bound resource
00127     MU_store_.bind( DMatrixSlice(NULL,0,0,0,0) ); // Unbind this!
00128     max_size_ = max_size;
00129     M_size_ = 0;
00130     M_l_r_ = M_l_c_ = 1;
00131     if( !maintain_factor && !allow_factor )
00132       U_l_r_ = 0;  // Do not allow the factor to be computed
00133     else
00134       U_l_r_ = 1;  // Allow the factor to be computed
00135     scale_ = +1.0;
00136     is_diagonal_ = false;
00137   }
00138   else {
00139     TEST_FOR_EXCEPT( !(  MU_store->rows()  ) );
00140     allocates_storage_ = false; // The client allocated the storage!
00141     MU_store_.bind(*MU_store);
00142     release_resource_ptr_ = release_resource_ptr;
00143     max_size_ = my_min( MU_store->rows(), MU_store->cols() ) - 1;
00144     if( set_full_view ) {
00145       TEST_FOR_EXCEPT( !(  scale != 0.0  ) );
00146       this->set_view(
00147         max_size_
00148         ,scale,maintain_original,1,1
00149         ,maintain_factor, allow_factor ? 1: 0, allow_factor ? 1 : 0
00150         );
00151     }
00152     else {
00153       this->set_view(
00154         0
00155         ,0.0,maintain_original,1,1
00156         ,maintain_factor, allow_factor ? 1: 0, allow_factor ? 1 : 0
00157         );
00158     }
00159   }
00160 }
00161 
00162 void MatrixSymPosDefCholFactor::set_view(
00163   size_t            M_size
00164   ,value_type       scale
00165   ,bool             maintain_original
00166   ,size_t           M_l_r
00167   ,size_t           M_l_c
00168   ,bool             maintain_factor
00169   ,size_t           U_l_r
00170   ,size_t           U_l_c
00171   )
00172 {
00173   TEST_FOR_EXCEPT( !(  maintain_original || maintain_factor  ) );
00174   if( max_size_ )
00175     allocate_storage(max_size_);
00176   else
00177     allocate_storage( my_max( M_l_r + M_size, M_l_c + M_size ) - 1 );
00178   // Check the preconditions
00179   if( maintain_original ) {
00180     TEST_FOR_EXCEPT( !(  1 <= M_l_r && M_l_r <= M_l_c  ) );
00181     TEST_FOR_EXCEPT( !(  M_l_r+M_size <= MU_store_.rows()  ) );
00182   }
00183   if( maintain_factor ) {
00184     TEST_FOR_EXCEPT( !(  1 <= U_l_r && U_l_r >= U_l_c  ) );
00185     TEST_FOR_EXCEPT( !(  U_l_c+M_size <= MU_store_.cols()  ) );
00186   }
00187   // Set the members
00188   maintain_original_    = maintain_original;
00189   maintain_factor_      = maintain_factor;
00190   is_diagonal_          = false;
00191   if( M_size ) {
00192     max_size_             = my_min( MU_store_.rows() - U_l_r, MU_store_.cols() - U_l_c );
00193     M_size_               = M_size;
00194     scale_                = scale;
00195     M_l_r_                = M_l_r;
00196     M_l_c_                = M_l_c;
00197     U_l_r_                = U_l_r;
00198     U_l_c_                = U_l_c;
00199     factor_is_updated_    = maintain_factor;
00200   }
00201   else {
00202     max_size_             = 0;
00203     M_size_               = 0;
00204     scale_                = 0.0;
00205     M_l_r_                = 0;
00206     M_l_c_                = 0;
00207     U_l_r_                = U_l_r;
00208     U_l_c_                = U_l_r;
00209     factor_is_updated_    = maintain_factor;
00210   }
00211 }
00212 
00213 void MatrixSymPosDefCholFactor::pivot_tols( PivotTolerances pivot_tols )
00214 {
00215   pivot_tols_ = pivot_tols;
00216 }
00217 
00218 MatrixSymAddDelUpdateable::PivotTolerances MatrixSymPosDefCholFactor::pivot_tols() const
00219 {
00220   return pivot_tols_;
00221 }
00222 
00223 // Overridden from MatrixBase
00224 
00225 size_type MatrixSymPosDefCholFactor::rows() const
00226 {
00227   return M_size_;
00228 }
00229 
00230 // Overridden from MatrixOp
00231 
00232 void MatrixSymPosDefCholFactor::zero_out()
00233 {
00234   this->init_identity(this->space_cols(),0.0);
00235 }
00236 
00237 std::ostream& MatrixSymPosDefCholFactor::output(std::ostream& out_arg) const
00238 {
00239   Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
00240   Teuchos::OSTab tab(out);
00241   if( M_size_ ) {
00242     if( maintain_original_ ) {
00243       *out
00244         << "Unfactored symmetric matrix stored as lower triangle (ignore upper nonzeros):\n"
00245         << M().gms();
00246     }
00247     if( factor_is_updated_ ) {
00248       *out
00249         << "Matrix scaling M = scale*U'*U, scale = " << scale_ << std::endl
00250         << "Upper cholesky factor U (ignore lower nonzeros):\n"
00251         << U().gms();
00252     }
00253   }
00254   else {
00255     *out << "0 0\n";
00256   }
00257   return out_arg;
00258 }
00259 
00260 bool MatrixSymPosDefCholFactor::Mp_StM(
00261   MatrixOp* m_lhs, value_type alpha
00262   ,BLAS_Cpp::Transp trans_rhs
00263   ) const
00264 {
00265   MatrixSymOpGetGMSSymMutable
00266     *symwo_gms_lhs = dynamic_cast<MatrixSymOpGetGMSSymMutable*>(m_lhs);
00267   if(!symwo_gms_lhs)
00268     return false;
00269   MatrixDenseSymMutableEncap sym_lhs(symwo_gms_lhs);
00270   const DMatrixSliceSym M       = this->M();
00271   DenseLinAlgPack::Mp_StM(
00272     &DMatrixSliceTriEle(sym_lhs().gms(),sym_lhs().uplo())
00273     ,alpha
00274     ,DMatrixSliceTriEle(M.gms(),M.uplo())
00275     );
00276 
00277   return true;
00278 }
00279 
00280 bool MatrixSymPosDefCholFactor::Mp_StM(
00281   value_type alpha,const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs
00282   )
00283 {
00284   TEST_FOR_EXCEPTION(
00285     !maintain_original_, std::logic_error
00286     ,"MatrixSymPosDefCholFactor::Mp_StM(alpha,M_rhs,trans_rhs): Error, Current implementation "
00287     "can not perform this operation unless the original matrix is being maintained." );
00288   // Perform the operation
00289   bool did_op  = false;
00290   bool diag_op = false;
00291   if(const MatrixSymOpGetGMSSym *symwo_gms_rhs = dynamic_cast<const MatrixSymOpGetGMSSym*>(&M_rhs)) {
00292     DMatrixSliceSym               M = this->M();
00293     MatrixDenseSymEncap   sym_rhs(*symwo_gms_rhs);
00294     DenseLinAlgPack::Mp_StM(
00295       &DMatrixSliceTriEle(M.gms(),M.uplo())
00296       ,alpha
00297       ,DMatrixSliceTriEle(sym_rhs().gms(),sym_rhs().uplo())
00298       );
00299     did_op  = true;
00300     diag_op = false;
00301   }
00302   else if(const MatrixSymDiag *symwo_diag_rhs = dynamic_cast<const MatrixSymDiag*>(&M_rhs)) {
00303     DMatrixSliceSym            M = this->M();
00304     VectorDenseEncap   sym_rhs_diag(symwo_diag_rhs->diag());
00305     LinAlgOpPack::Vp_StV( &M.gms().diag(), alpha, sym_rhs_diag() );
00306     did_op  = true;
00307     diag_op = true;
00308   }
00309   else if(const MatrixSymOp *symwo_rhs = dynamic_cast<const MatrixSymOp*>(&M_rhs)) {
00310     // ToDo: Implement this!
00311   }
00312   // Properly update the state of *this.
00313   // If only the original is updated
00314   if(did_op) {
00315     if( diag_op && is_diagonal_ )
00316       this->init_diagonal(VectorMutableDense(this->M().gms().diag(),Teuchos::null));
00317     else
00318       this->initialize(this->M());
00319     return true;
00320   }
00321   return false;
00322 }
00323 
00324 bool MatrixSymPosDefCholFactor::syrk(
00325   const MatrixOp      &mwo_rhs
00326   ,BLAS_Cpp::Transp   M_trans
00327   ,value_type         alpha
00328   ,value_type         beta
00329   )
00330 {
00331   MatrixDenseSymMutableEncap  sym_gms_lhs(this);
00332   const MatrixOpGetGMS *mwo_rhs_gms = dynamic_cast<const MatrixOpGetGMS*>(&mwo_rhs);
00333   if(mwo_rhs_gms) {
00334     TEST_FOR_EXCEPT(true); // ToDo: Implement
00335     return true;
00336   }
00337   else {
00338     // Here we will give up on symmetry and just compute the whole product S = op(mwo_rhs)*op(mwo_rhs')
00339     DenseLinAlgPack::DMatrixSliceTriEle tri_ele_gms_lhs = DenseLinAlgPack::tri_ele(sym_gms_lhs().gms(),sym_gms_lhs().uplo());
00340     if(beta==0.0)        DenseLinAlgPack::assign( &tri_ele_gms_lhs, 0.0 );
00341     else if(beta!=1.0)   DenseLinAlgPack::Mt_S( &tri_ele_gms_lhs, beta );
00342     const VectorSpace                 &spc = mwo_rhs.space_rows();
00343     const index_type                  m    = spc.dim();
00344     VectorSpace::multi_vec_mut_ptr_t  S    = spc.create_members(m);
00345     S->zero_out();
00346     LinAlgOpPack::M_MtM( S.get(), mwo_rhs, M_trans, mwo_rhs, BLAS_Cpp::trans_not(M_trans) );
00347     // Copy S into sym_ghs_lhs
00348     if( sym_gms_lhs().uplo() == BLAS_Cpp::lower ) {
00349       for( index_type j = 1; j <= m; ++j ) {
00350         LinAlgOpPack::Vp_StV( &sym_gms_lhs().gms().col(j)(j,m), alpha, VectorDenseEncap(*S->col(j)->sub_view(j,m))() );
00351       }
00352     }
00353     else {
00354       for( index_type j = 1; j <= m; ++j ) {
00355         LinAlgOpPack::Vp_StV( &sym_gms_lhs().gms().col(j)(1,j), alpha, VectorDenseEncap(*S->col(j)->sub_view(1,j))() );
00356       }
00357     }
00358     return true;
00359   }
00360   return false;
00361 }
00362 
00363 // Overridden from MatrixOpSerial
00364 
00365 void MatrixSymPosDefCholFactor::Vp_StMtV(
00366   DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00367   ,const DVectorSlice& x, value_type b
00368   ) const
00369 {
00370   using BLAS_Cpp::no_trans;
00371   using BLAS_Cpp::trans;
00372 #ifdef PROFILE_HACK_ENABLED
00373   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StMtV(...DVectorSlice...)" );
00374 #endif
00375   assert_initialized();
00376 
00377   DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), rows(), cols(), no_trans, x.dim() );
00378 
00379   if( maintain_original_ ) {
00380     //
00381     // M = symmetric
00382     //
00383     // y = b*y + a*M*x
00384     //
00385     DenseLinAlgPack::Vp_StMtV( y, a, M(), no_trans, x, b );
00386   }
00387   else {
00388     //
00389     // M = scale*U'*U
00390     //
00391     // y = b*y + a*op(M)*x
00392     //   = b*y = scale*a*U'*U*x
00393     //
00394     const DMatrixSliceTri
00395       U = this->U();
00396 
00397     if( b == 0.0 ) {
00398       // No temporary needed
00399       //
00400       // y = U*x
00401       DenseLinAlgPack::V_MtV( y, U, no_trans, x ); 
00402       // y = U'*y
00403       DenseLinAlgPack::V_MtV( y, U, trans, *y ); 
00404       // y *= scale*a
00405       if( a != 1.0 || scale_ != 1.0 )
00406         DenseLinAlgPack::Vt_S( y, scale_*a );
00407     }
00408     else {
00409       // We need a temporary
00410       DVector t;
00411       // t = U*x
00412       DenseLinAlgPack::V_MtV( &t, U, no_trans, x );
00413       // t = U'*t
00414       DenseLinAlgPack::V_MtV( &t(), U, trans, t() );
00415       // y *= b
00416       if(b != 1.0)
00417         DenseLinAlgPack::Vt_S( y, b );
00418       // y += scale*a*t
00419       DenseLinAlgPack::Vp_StV( y, scale_*a, t() );
00420     }
00421   }
00422 }
00423 
00424 void MatrixSymPosDefCholFactor::Vp_StMtV(
00425   DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00426   ,const SpVectorSlice& x, value_type b
00427   ) const
00428 {
00429   using BLAS_Cpp::no_trans;
00430   using BLAS_Cpp::trans;
00431 #ifdef PROFILE_HACK_ENABLED
00432   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StMtV(...SpVectorSlice...)" );
00433 #endif
00434   assert_initialized();
00435   if( maintain_original_ ) {
00436     const DMatrixSlice M = this->M().gms(); // This is lower triangular!
00437     const size_type n = M.rows();
00438     DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), n, n, no_trans, x.dim() );
00439     DenseLinAlgPack::Vt_S(y,b); // y = b*y
00440     //
00441     // Compute product column by column corresponding to x_itr->index() + x.offset()
00442     //
00443     // y += a * M * e(i) * x(i)
00444     //
00445     // [ y(1:i-1) ] += a * x(i) * [ ...  M(1:i-1,i) ... ]  stored as M(i,1:i-1) 
00446     // [ y(i:n)   ]               [ ...  M(i:n,i)   ... ]  stored as M(i:n,i)
00447     //
00448     for( SpVectorSlice::const_iterator x_itr = x.begin(); x_itr != x.end(); ++x_itr ) {
00449       const size_type i = x_itr->index() + x.offset();
00450       if( i > 1 )
00451         DenseLinAlgPack::Vp_StV( &(*y)(1,i-1), a * x_itr->value(), M.row(i)(1,i-1) );
00452       DenseLinAlgPack::Vp_StV( &(*y)(i,n), a * x_itr->value(), M.col(i)(i,n) );
00453     }
00454   }
00455   else {
00456     MatrixOpSerial::Vp_StMtV(y,a,M_trans,x,b); // ToDo: Specialize when needed!
00457   }
00458 }
00459 
00460 void MatrixSymPosDefCholFactor::Vp_StPtMtV(
00461   DVectorSlice* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00462   , BLAS_Cpp::Transp H_trans, const DVectorSlice& x, value_type b
00463   ) const
00464 {
00465 #ifdef PROFILE_HACK_ENABLED
00466   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StPtMtV(...DVectorSlice...)" );
00467 #endif
00468   assert_initialized();
00469   MatrixOpSerial::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // ToDo: Specialize when needed!
00470 }
00471 
00472 void MatrixSymPosDefCholFactor::Vp_StPtMtV(
00473   DVectorSlice* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00474   , BLAS_Cpp::Transp H_trans, const SpVectorSlice& x, value_type b
00475   ) const
00476 {
00477 #ifdef PROFILE_HACK_ENABLED
00478   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StPtMtV(...SpVectorSlice...)" );
00479 #endif
00480   assert_initialized();
00481   if( maintain_original_ ) {
00482     DenseLinAlgPack::Vt_S(y,b); // y = b*y
00483     const DMatrixSlice M = this->M().gms(); // This is lower triangular!
00484     // Compute product column by corresponding to x_itr->index() + x.offset()
00485     /*
00486     if( P.is_identity() ) {
00487       TEST_FOR_EXCEPT(true); // ToDo: Implement
00488     }
00489     else {
00490       for( SpVectorSlice::const_iterator x_itr = x.begin(); x_itr != x.end(); ++x_itr ) {
00491         const size_type i = x_itr->index() + x.offset();
00492         
00493 
00494       }
00495     }
00496     */
00497     MatrixOpSerial::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // ToDo: Specialize when needed!
00498   }
00499   else {
00500     MatrixOpSerial::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // ToDo: Specialize when needed!
00501   }
00502 }
00503 
00504 // Overridden from MatrixSymOpSerial
00505 
00506 void MatrixSymPosDefCholFactor::Mp_StPtMtP(
00507   DMatrixSliceSym* S, value_type a
00508   , EMatRhsPlaceHolder dummy_place_holder
00509   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00510   , value_type b ) const
00511 {
00512   using BLAS_Cpp::no_trans;
00513   using BLAS_Cpp::trans;
00514 #ifdef PROFILE_HACK_ENABLED
00515   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Mp_StPtMtP(...)" );
00516 #endif
00517   assert_initialized();
00518   if( !maintain_original_ ) {
00519     MatrixSymOpSerial::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b);
00520   }
00521   else {
00522     MatrixSymOpSerial::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b); // ToDo: Override when needed!
00523   }
00524 }
00525 
00526 // Overridden from MatrixNonsingSerial
00527 
00528 void MatrixSymPosDefCholFactor::V_InvMtV(
00529   DVectorSlice* y, BLAS_Cpp::Transp M_trans, const DVectorSlice& x
00530   ) const
00531 {
00532   using BLAS_Cpp::no_trans;
00533   using BLAS_Cpp::trans;
00534 #ifdef PROFILE_HACK_ENABLED
00535   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::V_InvMtV(...DVectorSlice...)" );
00536 #endif
00537   assert_initialized();
00538 
00539   //
00540   // M = scale*U'*U
00541   //
00542   // y = inv(op(M))*x
00543   // =>
00544   // op(M)*y = x
00545   // =>
00546   // scale*U'*U*y = x
00547   // =>
00548   // y = (1/scale)*inv(U)*inv(U')*x
00549   //
00550   update_factorization();
00551   const DMatrixSliceTri
00552     U = this->U();
00553   DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), U.rows(), U.cols(), no_trans, x.dim() );
00554   // y = inv(U')*x
00555   DenseLinAlgPack::V_InvMtV( y, U, trans, x );
00556   // y = inv(U)*y
00557   DenseLinAlgPack::V_InvMtV( y, U, no_trans, *y );
00558   // y *= (1/scale)
00559   if( scale_ != 1.0 )
00560     DenseLinAlgPack::Vt_S( y, 1.0/scale_ );
00561 }
00562 
00563 void MatrixSymPosDefCholFactor::V_InvMtV(
00564   DVectorSlice* y, BLAS_Cpp::Transp M_trans, const SpVectorSlice& x
00565   ) const
00566 {
00567 #ifdef PROFILE_HACK_ENABLED
00568   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::V_InvMtV(...SpVectorSlice...)" );
00569 #endif
00570   assert_initialized();
00571   MatrixNonsingSerial::V_InvMtV(y,M_trans,x);
00572 }
00573 
00574 // Overridden from MatrixSymNonsingSerial
00575 
00576 void MatrixSymPosDefCholFactor::M_StMtInvMtM(
00577     DMatrixSliceSym* S, value_type a, const MatrixOpSerial& B
00578   , BLAS_Cpp::Transp B_trans, EMatrixDummyArg dummy_arg
00579   ) const
00580 {
00581 #ifdef PROFILE_HACK_ENABLED
00582   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::M_StMtInvMtM(...)" );
00583 #endif
00584 
00585 //  // Uncomment to use the defalut implementation (for debugging)
00586 //  MatrixSymFactorized::M_StMtInvMtM(S,a,B,B_trans,dummy_arg); return;
00587 
00588   using BLAS_Cpp::trans;
00589   using BLAS_Cpp::no_trans;
00590   using BLAS_Cpp::trans_not;
00591   using BLAS_Cpp::upper;
00592   using BLAS_Cpp::nonunit;
00593   using DenseLinAlgPack::tri;
00594   using DenseLinAlgPack::syrk;
00595   using DenseLinAlgPack::M_StInvMtM;
00596   using LinAlgOpPack::assign;
00597 
00598   assert_initialized();
00599   update_factorization();
00600 
00601   //
00602   // S = a * op(B) * inv(M) * op(B)'
00603   // 
00604   // M = scale*U'*U
00605   // =>
00606   // inv(M) = scale*inv(U'*U) = scale*inv(U)*inv(U')
00607   // =>
00608   // S = scale*a * op(B) * inv(U) * inv(U') * op(B)'
00609   // 
00610   // T = op(B)'
00611   // 
00612   // T = inv(U') * T (inplace with BLAS)
00613   // 
00614   // S = scale*a * T' * T
00615   // 
00616   DenseLinAlgPack::MtM_assert_sizes( 
00617     rows(), cols(), no_trans
00618     ,B.rows(), B.cols(), trans_not(B_trans)
00619     );
00620   DenseLinAlgPack::Mp_MtM_assert_sizes(
00621     S->rows(), S->cols(), no_trans
00622     ,B.rows(), B.cols(), B_trans
00623     ,B.rows(), B.cols(), trans_not(B_trans)
00624     );
00625   // T = op(B)'
00626   DMatrix T;
00627   assign( &T, B, trans_not(B_trans) );
00628   // T = inv(U') * T (inplace with BLAS)
00629   M_StInvMtM( &T(), 1.0, this->U(), trans, T(), no_trans );
00630   // S = scale*a * T' * T
00631   syrk( trans, scale_*a, T(), 0.0, S );
00632 }
00633 
00634 // Overridden from MatrixSymDenseInitialize
00635 
00636 void MatrixSymPosDefCholFactor::initialize( const DMatrixSliceSym& M )
00637 {
00638   // Initialize without knowing the inertia but is must be p.d.
00639   this->initialize(
00640     M, M.rows(), maintain_factor_
00641     , MatrixSymAddDelUpdateable::Inertia()
00642     , MatrixSymAddDelUpdateable::PivotTolerances()
00643     );
00644 }
00645 
00646 // Overridden from MatrixSymOpGetGMSSym
00647 
00648 const DenseLinAlgPack::DMatrixSliceSym MatrixSymPosDefCholFactor::get_sym_gms_view() const
00649 {
00650   TEST_FOR_EXCEPTION(
00651     !maintain_original_, std::logic_error
00652     ,"MatrixSymPosDefCholFactor::get_sym_gms_view(): Error, maintain_original must be "
00653     "true in order to call this method!" );
00654   return this->M();
00655 }
00656 
00657 void MatrixSymPosDefCholFactor::free_sym_gms_view(const DenseLinAlgPack::DMatrixSliceSym* sym_gms_view) const
00658 {
00659   TEST_FOR_EXCEPTION(
00660     !maintain_original_, std::logic_error
00661     ,"MatrixSymPosDefCholFactor::free_sym_gms_view(...): Error, maintain_original must be "
00662     "true in order to call this method!" );
00663   // Nothing todo
00664 }
00665 
00666 // Overridden from MatrixSymOpGetGMSSymMutable
00667 
00668 DenseLinAlgPack::DMatrixSliceSym MatrixSymPosDefCholFactor::get_sym_gms_view()
00669 {
00670   TEST_FOR_EXCEPTION(
00671     !maintain_original_, std::logic_error
00672     ,"MatrixSymPosDefCholFactor::get_sym_gms_view(): Error, maintain_original must be "
00673     "true in order to call this method!" );
00674   return this->M();
00675 }
00676 
00677 void MatrixSymPosDefCholFactor::commit_sym_gms_view(DenseLinAlgPack::DMatrixSliceSym* sym_gms_view)
00678 {
00679   TEST_FOR_EXCEPTION(
00680     !maintain_original_, std::logic_error
00681     ,"MatrixSymPosDefCholFactor::commit_sym_gms_view(...): Error, maintain_original must be "
00682     "true in order to call this method!" );
00683   this->initialize(*sym_gms_view);
00684 }
00685 
00686 // Overridden from MatrixExtractInvCholFactor
00687 
00688 void MatrixSymPosDefCholFactor::extract_inv_chol( DMatrixSliceTriEle* InvChol ) const
00689 {
00690   assert_initialized();
00691   update_factorization();
00692   //
00693   // The matrix is represented as the upper cholesky factor:
00694   //   M = scale * U' * U
00695   //
00696   //   inv(M) = inv(scale*U*U') = (1/sqrt(scale))*inv(U)*(1/sqrt(scale))*inv(U')
00697   //          = UInv * UInv'
00698   // =>
00699   //   UInv = (1/sqrt(scale))*inv(U)
00700   //
00701    // Here scale > 0 or an exception will be thrown.
00702   //
00703   // Compute the inverse cholesky factor as:
00704   //
00705   // Upper cholesky:
00706   //    sqrt(scale) * U * UInv = I => InvChol = UInv = (1/sqrt(scale))*inv(U) * I
00707   //    
00708   // Lower cholesky:
00709   //    sqrt(scale) * L * LInv = I => InvChol = LInv = (1/sqrt(scale))*inv(U) * inv(U') * I
00710   //
00711   TEST_FOR_EXCEPTION(
00712     scale_ < 0.0, std::logic_error
00713     ,"MatrixSymPosDefCholFactor::extract_inv_chol(...) : "
00714     "Error, we can not compute the inverse cholesky factor "
00715     "af a negative definite matrix." );
00716   DenseLinAlgPack::assign( &InvChol->gms(), 0.0 );  // Set InvChol to identity first.
00717   InvChol->gms().diag() = 1.0;
00718   DenseLinAlgPack::M_StInvMtM(                      // Comput InvChol using Level-3 BLAS
00719       &InvChol->gms(), 1.0 / std::sqrt(scale_), U()
00720     , InvChol->uplo() == BLAS_Cpp::upper ? BLAS_Cpp::no_trans : BLAS_Cpp::trans
00721     , InvChol->gms(), BLAS_Cpp::no_trans );
00722 }
00723 
00724 // Overridden from MatrixSymSecantUpdateble
00725 
00726 void MatrixSymPosDefCholFactor::init_identity( const VectorSpace& space_diag, value_type alpha )
00727 {
00728   const size_type n = space_diag.dim();
00729   allocate_storage( max_size_ ? max_size_ : n );
00730   //
00731   // B = alpha*I = = alpha*I*I = scale*U'*U
00732   // =>
00733   // U = sqrt(|alpha|) * I
00734   //
00735   // Here we will set scale = sign(alpha)
00736   //
00737   const value_type scale = alpha > 0.0 ? +1.0: -1.0; // Explicitly set the scale
00738   resize_and_zero_off_diagonal(n,scale);
00739   if( maintain_original_ ) {
00740     M().gms().diag()(1,n) = alpha;
00741   }
00742   if( maintain_factor_ ) {
00743     U().gms().diag()(1,n) = std::sqrt(std::fabs(alpha));
00744     factor_is_updated_ = true;
00745   }
00746   is_diagonal_ = true;
00747 }
00748 
00749 void MatrixSymPosDefCholFactor::init_diagonal( const Vector& diag_in )
00750 {
00751   VectorDenseEncap diag_encap(diag_in);
00752   const DVectorSlice diag = diag_encap(); // When diag_encap is destroyed, bye-bye view!
00753 
00754   allocate_storage( max_size_ ? max_size_ : diag.dim() );
00755   //
00756   // M = scale * U' * U = scale * (1/scale)*diag^(1/2) * (1/scale)*diag^(1/2)
00757   //
00758   // Here we will set scale = sign(diag(1)) and validate the rest
00759   //
00760   if( diag.dim() == 0 ) {
00761     M_size_ = 0;
00762     return; // We are unsizing this thing
00763   }
00764   const value_type scale = diag(1) > 0.0 ? +1.0: -1.0;
00765   resize_and_zero_off_diagonal(diag.dim(),scale);
00766   if( maintain_original_ ) {
00767     M().gms().diag() = diag;
00768     // ToDo: validate that scale*diag > 0
00769   }
00770   if( maintain_factor_ ) {
00771     DVectorSlice U_diag = U().gms().diag();
00772     U_diag = diag;
00773     if( scale_ != 1.0 )
00774       DenseLinAlgPack::Vt_S( &U_diag, 1.0/scale_ );
00775     DenseLinAlgPack::sqrt( &U_diag, U_diag );
00776     DenseLinAlgPack::assert_print_nan_inf( U_diag, "(1/scale)*diag", true, NULL );
00777     factor_is_updated_ = true;
00778   }
00779   is_diagonal_ = true;
00780 }
00781 
00782 void MatrixSymPosDefCholFactor::secant_update(
00783   VectorMutable* s_in, VectorMutable* y_in, VectorMutable* Bs_in
00784   )
00785 {
00786   using BLAS_Cpp::no_trans;
00787   using BLAS_Cpp::trans;
00788   using DenseLinAlgPack::dot;
00789   using DenseLinAlgPack::norm_2;
00790   using DenseLinAlgPack::norm_inf;
00791   namespace rcp = MemMngPack;
00792 
00793   assert_initialized();
00794 
00795   // Validate the input
00796   TEST_FOR_EXCEPT( !(  s_in && y_in  ) );
00797   DenseLinAlgPack::Vp_MtV_assert_sizes( y_in->dim(), M_size_, M_size_, no_trans, s_in->dim() );
00798 
00799   // Get the serial vectors
00800   VectorDenseMutableEncap s_encap(*s_in);
00801   VectorDenseMutableEncap y_encap(*y_in);
00802   VectorDenseMutableEncap Bs_encap( Bs_in ? *Bs_in : *y_in); // Must pass something on
00803   DVectorSlice
00804     *s  = &s_encap(),   // When s_encap, y_encap and Bs_encap are destroyed
00805     *y  = &y_encap(),   // these views go bye-bye!
00806     *Bs = ( Bs_in ? &Bs_encap() : NULL );
00807 
00808   // Check skipping the BFGS update
00809   const value_type
00810     sTy       = dot(*s,*y),
00811     sTy_scale = sTy*scale_;
00812   std::ostringstream omsg;
00813   if( !BFGS_sTy_suff_p_d(
00814       VectorMutableDense((*s)(),Teuchos::null)
00815       ,VectorMutableDense((*y)(),Teuchos::null)
00816       ,&sTy_scale,&omsg,"\nMatrixSymPosDefCholFactor::secant_update(...)"
00817       )
00818     )
00819   {
00820     throw UpdateSkippedException( omsg.str() ); 
00821   }
00822   // Compute Bs if it was not passed in
00823   DVector Bs_store;
00824   DVectorSlice Bs_view;
00825   if( !Bs ) {
00826     LinAlgOpPack::V_MtV( &Bs_store, *this, no_trans, *s );
00827     Bs_view.bind( Bs_store() );
00828     Bs = &Bs_view;
00829   }
00830   // Check that s'*Bs is positive and if not then throw exception
00831   const value_type sTBs = dot(*s,*Bs);
00832   TEST_FOR_EXCEPTION(
00833     scale_*sTBs <= 0.0 && scale_ > 0.0, UpdateFailedException
00834     ,"MatrixSymPosDefCholFactor::secant_update(...) : "
00835     "Error, B can't be positive definite if s'*Bs <= 0.0" );
00836   TEST_FOR_EXCEPTION(
00837     scale_*sTBs <= 0.0 && scale_ <= 0.0, UpdateFailedException
00838     ,"MatrixSymPosDefCholFactor::secant_update(...) : "
00839     "Error, B can't be negative definite if s'*Bs >= 0.0" );
00840   if( maintain_original_ ) {
00841     //
00842     // Compute the BFGS update of the original, nonfactored matrix
00843     //
00844     // Just preform two symmetric updates.
00845     //
00846     // B = B + (-1/s'*Bs) * Bs*Bs' + (1/s'*y) * y*y'
00847     //
00848     DMatrixSliceSym M = this->M();
00849     DenseLinAlgPack::syr( -1.0/sTBs, *Bs, &M );
00850     DenseLinAlgPack::syr( 1.0/sTy, *y, &M );
00851   }
00852   if( maintain_factor_ ) {
00853     //
00854     // Compute the BFGS update for the cholesky factor
00855     //
00856     // If this implementation is based on the one in Section 9.2, page 198-201 of:
00857     //
00858     // Dennis, J.E., R.B. Schnabel
00859     // "Numerical Methods for Unconstrained Optimization"
00860     //
00861     // Given that we have B = scale*U'*U and the BFGS update:
00862     //
00863     // B_new = B + y*y'/(y'*s) - (B*s)*(B*s)'/(s'*B*s)
00864     //
00865     // We can rewrite it in the following form:
00866     //
00867     // B_new = scale*(U' + a*u*v')*(U + a*v*u') = scale*U_new'*U_new
00868     //
00869     // where:
00870     //     v = sqrt(y'*s/(s'*B*s))*U*s
00871     //     u = y - U'*v
00872     //     a = 1/(v'*v)
00873     //
00874     DMatrixSliceTri U = this->U();
00875     // v = sqrt(y'*s/(s'*B*s))*U*s
00876     DVectorSlice v = *s; // Reuse s as storage for v
00877     DenseLinAlgPack::V_MtV( &v, U, no_trans, v ); // Direct call to xSYMV(...)
00878     DenseLinAlgPack::Vt_S( &v, std::sqrt( sTy / sTBs ) );
00879     // u = (y - U'*v)
00880     DVectorSlice u = *y;  // Reuse y as storage for u
00881     DenseLinAlgPack::Vp_StMtV( &u, -1.0, U, trans, v );
00882     // a = 1/(v'*v)
00883     const value_type  a = 1.0/dot(v,v);
00884     // Perform Givens rotations to make Q*(U' + a*u*v') -> U_new upper triangular:
00885     //
00886     // B_new = scale*(U' + a*u*v')*Q'*Q*(U + a*v*u') = scale*U_new'*U_new
00887     rank_2_chol_update(
00888       a, &v, u, v.dim() > 1 ? &U.gms().diag(-1) : NULL
00889       , &DMatrixSliceTriEle(U.gms(),BLAS_Cpp::upper), no_trans );
00890   }
00891   else {
00892     factor_is_updated_ = false;
00893   }
00894   is_diagonal_ = false;
00895 }
00896 
00897 // Overridden from MatrixSymAddDelUpdateble
00898 
00899 void MatrixSymPosDefCholFactor::initialize(
00900   value_type         alpha
00901   ,size_type         max_size
00902   )
00903 {
00904 #ifdef PROFILE_HACK_ENABLED
00905   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::initialize(alpha,max_size)" );
00906 #endif
00907   allocate_storage(max_size);
00908 
00909   if( alpha == 0.0 ) 
00910     throw SingularUpdateException(
00911       "MatrixSymPosDefCholFactor::initialize(...): "
00912       "Error, alpha == 0.0, matrix is singular.", 0.0 );
00913 
00914   // Resize the view
00915   if( maintain_original_ ) {
00916     M_l_r_ = 1;
00917     M_l_c_ = 1;
00918   }
00919   if( U_l_r_ ) {
00920     U_l_r_ = 1;
00921     U_l_c_ = 1;
00922   }
00923   M_size_ = 1;
00924   max_size_ = my_min( MU_store_.rows(), MU_store_.cols() ) - 1;
00925   scale_ = alpha > 0.0 ? +1.0 : -1.0;
00926   // Update the matrix
00927   if( maintain_original_ ) {
00928     M().gms()(1,1) = alpha;
00929   }
00930   if( U_l_r_ ) {
00931     U().gms()(1,1) = std::sqrt( scale_ * alpha );
00932     factor_is_updated_ = true;
00933   }
00934   is_diagonal_ = false;
00935 }
00936 
00937 void MatrixSymPosDefCholFactor::initialize(
00938   const DMatrixSliceSym      &A
00939   ,size_type         max_size
00940   ,bool              force_factorization
00941   ,Inertia           inertia
00942   ,PivotTolerances   pivot_tols
00943   )
00944 {
00945 #ifdef PROFILE_HACK_ENABLED
00946   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::initialize(A,max_size...)" );
00947 #endif
00948   typedef MatrixSymAddDelUpdateable::Inertia   Inertia;
00949 
00950   allocate_storage(max_size);
00951 
00952   const size_type
00953     n = A.rows();
00954 
00955   // Validate proper usage of inertia parameter
00956   TEST_FOR_EXCEPT( !(  inertia.zero_eigens == Inertia::UNKNOWN || inertia.zero_eigens == 0  ) );
00957   TEST_FOR_EXCEPT( !(  (inertia.neg_eigens == Inertia::UNKNOWN && inertia.pos_eigens == Inertia::UNKNOWN ) )
00958       || ( inertia.neg_eigens == n && inertia.pos_eigens == 0 )
00959       || ( inertia.neg_eigens == 0 && inertia.pos_eigens == n )
00960     );
00961 
00962   // We can infer if the matrix is p.d. or n.d. by the sign of the diagonal
00963   // elements.  If a matrix is s.p.d. (s.n.d) then A(i,i) > 0 (A(i,i) < 0)
00964   // for all i = 1...n so we will just check the first i.
00965   const value_type
00966     a_11 = A.gms()(1,1);
00967   const int sign_a_11 = ( a_11 == 0.0 ? 0 : ( a_11 > 0.0 ? +1 : -1 ) );
00968   if( sign_a_11 == 0.0 )
00969     std::logic_error(
00970       "MatrixSymPosDefCholFactor::initialize(...) : "
00971       "Error, A can not be positive definite or negative definete "
00972       "if A(1,1) == 0.0" );
00973   if( inertia.pos_eigens == n && sign_a_11 < 0 )
00974     std::logic_error(
00975       "MatrixSymPosDefCholFactor::initialize(...) : "
00976       "Error, A can not be positive definite "
00977       "if A(1,1) < 0.0" );
00978   if( inertia.neg_eigens == n && sign_a_11 > 0 )
00979     std::logic_error(
00980       "MatrixSymPosDefCholFactor::initialize(...) : "
00981       "Error, A can not be negative definite "
00982       "if A(1,1) > 0.0" );
00983   // Now we have got the sign
00984   const value_type
00985     scale = (value_type)sign_a_11;
00986   // Setup the view
00987   set_view(
00988     n,scale,maintain_original_,1,1
00989     ,maintain_factor_, U_l_r_ ? 1 : 0, U_l_r_ ? 1 : 0
00990     );
00991   // Now set the matrix and update the factors
00992   if( maintain_original_ ) {
00993     // Set M = S
00994     DenseLinAlgPack::assign( 
00995       &DMatrixSliceTriEle( M().gms(), BLAS_Cpp::lower )
00996       ,DMatrixSliceTriEle( A.gms(), A.uplo() )
00997       );
00998   }
00999   if( maintain_factor_ || force_factorization ) {
01000     // Copy S into U for an inplace factorization.
01001     DMatrixSliceTriEle U_ele = DMatrixSliceTriEle( U().gms(), BLAS_Cpp::upper );
01002     DenseLinAlgPack::assign( &U_ele, DMatrixSliceTriEle( A.gms(), A.uplo() ) );
01003     if( sign_a_11 < 0 )
01004       DenseLinAlgPack::Mt_S( &U_ele, -1.0 );
01005     try {
01006       DenseLinAlgLAPack::potrf( &U_ele );
01007       factor_is_updated_ = true;
01008     }
01009     catch( const DenseLinAlgLAPack::FactorizationException& excpt ) {
01010       M_size_ = 0;  // set unsized
01011       throw SingularUpdateException( excpt.what(), 0.0 );
01012     }
01013     catch(...) {
01014       M_size_ = 0;
01015       throw;
01016     }
01017     // Validate that the tolerances are met and throw appropriate
01018     // exceptions.  We already know that the matrix is technically
01019     // p.d. or n.d.  Now we must determine gamma = (min(|diag|)/max(|diag|))^2
01020     value_type
01021       min_diag = std::numeric_limits<value_type>::max(),
01022       max_diag = 0.0;
01023     DVectorSlice::iterator
01024       U_itr = U_ele.gms().diag().begin(),
01025       U_end = U_ele.gms().diag().end();
01026     while( U_itr != U_end ) {
01027       const value_type U_abs = std::abs(*U_itr++);
01028       if(U_abs < min_diag) min_diag = U_abs;
01029       if(U_abs > max_diag) max_diag = U_abs;
01030     }
01031     const value_type gamma = (min_diag*min_diag)/(max_diag*max_diag);
01032     // Validate gamma
01033     PivotTolerances use_pivot_tols = pivot_tols_;
01034     if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
01035       use_pivot_tols.warning_tol = pivot_tols.warning_tol;
01036     if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
01037       use_pivot_tols.singular_tol = pivot_tols.singular_tol;
01038     if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
01039       use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
01040     const bool throw_exception = (
01041       gamma == 0.0
01042       || gamma <= use_pivot_tols.warning_tol
01043       || gamma <= use_pivot_tols.singular_tol
01044       );
01045     // Create message and throw exceptions
01046     std::ostringstream onum_msg, omsg;
01047     if(throw_exception) {
01048       onum_msg
01049         << "gamma = (min(|diag(U)(k)|)/|max(|diag(U)(k)|))^2 = (" << min_diag <<"/"
01050         << max_diag << ")^2 = " << gamma;
01051       omsg
01052         << "MatrixSymPosDefCholFactor::initialize(...): ";
01053       if( gamma <= use_pivot_tols.singular_tol ) {
01054         M_size_ = 0;  // The initialization failed!
01055         omsg
01056           << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
01057           << use_pivot_tols.singular_tol;
01058         throw SingularUpdateException( omsg.str(), gamma );
01059       }
01060       else if( gamma <= use_pivot_tols.warning_tol ) {
01061         omsg
01062           << "Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
01063           << " < " << onum_msg.str() << " <= warning_tol = "
01064           << use_pivot_tols.warning_tol;
01065         throw WarnNearSingularUpdateException( omsg.str(), gamma ); // The initialization still succeeded through
01066       }
01067       else {
01068         TEST_FOR_EXCEPT(true); // Only local programming error?
01069       }
01070     }
01071   }
01072   else {
01073     factor_is_updated_ = false; // The factor is not updated!
01074   }
01075 }
01076 
01077 size_type MatrixSymPosDefCholFactor::max_size() const
01078 {
01079   return max_size_;
01080 }
01081 
01082 MatrixSymAddDelUpdateable::Inertia
01083 MatrixSymPosDefCholFactor::inertia() const
01084 {
01085   typedef MatrixSymAddDelUpdateable MSADU;
01086   typedef MSADU::Inertia Inertia;
01087   return ( M_size_
01088        ? ( scale_ > 0.0
01089          ? Inertia(0,0,M_size_)
01090          : Inertia(M_size_,0,0) )
01091        : Inertia(0,0,0) );
01092 }
01093 
01094 void MatrixSymPosDefCholFactor::set_uninitialized()
01095 {
01096   M_size_ = 0;
01097 }
01098 
01099 void MatrixSymPosDefCholFactor::augment_update(
01100   const DVectorSlice  *t
01101   ,value_type        alpha
01102   ,bool              force_refactorization
01103   ,EEigenValType     add_eigen_val
01104   ,PivotTolerances   pivot_tols
01105   )
01106 {
01107 #ifdef PROFILE_HACK_ENABLED
01108   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::augment_udpate(...)" );
01109 #endif
01110   using DenseLinAlgPack::dot;
01111   using DenseLinAlgPack::norm_inf;
01112   typedef MatrixSymAddDelUpdateable  MSADU;
01113 
01114   assert_initialized();
01115 
01116   // Validate the input
01117   TEST_FOR_EXCEPTION(
01118     rows() >= max_size(), MaxSizeExceededException
01119     ,"MatrixSymPosDefCholFactor::augment_update(...) : "
01120     "Error, the maximum size would be exceeded." );
01121   TEST_FOR_EXCEPTION(
01122     t && t->dim() != M_size_, std::length_error
01123     ,"MatrixSymPosDefCholFactor::augment_update(...): "
01124     "Error, t.dim() must be equal to this->rows()." );
01125   if( !(add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN || add_eigen_val != MSADU::EIGEN_VAL_ZERO ) )
01126     throw SingularUpdateException(
01127       "MatrixSymPosDefCholFactor::augment_update(...): "
01128       "Error, the client has specified a singular update in add_eigen_val.", -1.0 );
01129   if( alpha == 0.0 )
01130     throw SingularUpdateException(
01131       "MatrixSymPosDefCholFactor::augment_update(...): "
01132       "Error, alpha == 0.0 so the matrix is not positive definite "
01133       "or negative definite.", -1.0 );
01134   if( scale_ > 0.0 && alpha < 0.0 )
01135     throw WrongInertiaUpdateException(
01136       "MatrixSymPosDefCholFactor::augment_update(...): "
01137       "Error, alpha < 0.0 so the matrix is not postivie definite ", -1.0 );
01138   if( scale_ < 0.0 && alpha > 0.0 )
01139     throw WrongInertiaUpdateException(
01140       "MatrixSymPosDefCholFactor::augment_update(...): "
01141       "Error, alpha > 0.0 so the matrix is not negative definite ", -1.0 );
01142   if( scale_ > 0.0 && !( add_eigen_val == MSADU::EIGEN_VAL_POS || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
01143     throw WrongInertiaUpdateException(
01144       "MatrixSymPosDefCholFactor::augment_update(...): "
01145       "Error, alpha > 0.0 but the user specified a non-positive definite new matrix.", -1.0 );
01146   if( scale_ < 0.0 && !( add_eigen_val == MSADU::EIGEN_VAL_NEG || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
01147     throw WrongInertiaUpdateException(
01148       "MatrixSymPosDefCholFactor::augment_update(...): "
01149       "Error, alpha > 0.0 but the user specified a non-positive definite new matrix.", -1.0 );
01150   // First try to augment the factor to verify that the matrix is still p.d. or n.d.
01151   bool throw_exception = false; // If true then throw exception
01152   std::ostringstream omsg;      // Will be set if an exception has to be thrown.
01153   value_type gamma;             // ...
01154   if( maintain_factor_ ) {
01155     //
01156     // The update is:
01157     //
01158     // B_new = [  B     t   ]
01159     //         [  t'  alpha ]
01160     //
01161     //       = scale * [     U'*U         (1/scale)*t   ]
01162     //                 [ (1/scale)*t'   (1/scale)*alpha ]
01163     //
01164     // We seek to find a new cholesky factor of the form:
01165     //
01166     // U_new = [ U11  u12  ]
01167     //         [  0   u22  ]
01168     //
01169     // B_new = scale*U_new'*U_new
01170     //
01171     //       = scale * [ U11'   0  ] * [ U11   u12 ]
01172     //                 [ u12'  u22 ]   [  0    u22 ]
01173     //
01174     //       = scale * [  U11'*U11          U11'*u12  ]
01175     //                 [  u12'*U11   u12'*u12 + u22^2 ]
01176     //
01177     // From the above we see that:
01178     // =>
01179     //   U11 = U
01180     //   u12 = inv(U') * (1/scale) * t
01181     //   u22 = sqrt( (1/scale)*alpha - u12'*u12 );
01182     //
01183     // We must compute gamma relative to the LU factorization
01184     // so we must square ||U11.diag()||inf.
01185 
01186     // Get references to the storage for the to-be-updated parts for the new factor.
01187     DVectorSlice u12 = MU_store_.col(U_l_c_+M_size_+1)(U_l_r_,U_l_r_+M_size_-1);
01188     value_type &u22 = MU_store_(U_l_r_+M_size_,U_l_c_+M_size_+1);
01189     // u12 = inv(U') * (1/scale) * t
01190     if(t) {
01191       DenseLinAlgPack::V_InvMtV( &u12, U(), BLAS_Cpp::trans, *t );
01192       if( scale_ != 1.0 ) DenseLinAlgPack::Vt_S( &u12, 1.0/scale_ );
01193     }
01194     else {
01195       u12 = 0.0;
01196     }
01197     // u22^2 = (1/scale)*alpha - u12'*u12;
01198     const value_type
01199       u22sqr          = (1/scale_) * alpha - ( t ? dot( u12, u12 ) : 0.0 ),
01200       u22sqrabs       = std::abs(u22sqr),
01201       nrm_U_diag      = norm_inf(U().gms().diag()),
01202       sqr_nrm_U_diag  = nrm_U_diag * nrm_U_diag;
01203     // Calculate gamma in proper context
01204     gamma = u22sqrabs / sqr_nrm_U_diag;
01205     // Check gamma
01206     const bool
01207       correct_inertia = u22sqr > 0.0;
01208     PivotTolerances use_pivot_tols = pivot_tols_;
01209     if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
01210       use_pivot_tols.warning_tol = pivot_tols.warning_tol;
01211     if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
01212       use_pivot_tols.singular_tol = pivot_tols.singular_tol;
01213     if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
01214       use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
01215     throw_exception = (
01216       gamma == 0.0
01217       || correct_inertia && gamma <= use_pivot_tols.warning_tol
01218       || correct_inertia && gamma <= use_pivot_tols.singular_tol
01219       || !correct_inertia
01220       );
01221     // Create message and throw exceptions
01222     std::ostringstream onum_msg;
01223     if(throw_exception) {
01224       onum_msg
01225         << "gamma = |u22^2|/(||diag(U11)||inf)^2 = |" << u22sqr <<"|/("
01226         << nrm_U_diag << ")^2 = " << gamma;
01227       omsg
01228         << "MatrixSymPosDefCholFactor::augment_update(...): ";
01229       if( correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
01230         omsg
01231           << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
01232           << use_pivot_tols.singular_tol;
01233         throw SingularUpdateException( omsg.str(), gamma );
01234       }
01235       else if( !correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
01236         omsg
01237           << "Singular update!\n" << onum_msg.str() << " <= wrong_inertia_tol = "
01238           << use_pivot_tols.wrong_inertia_tol;
01239         throw SingularUpdateException( omsg.str(), gamma );
01240       }
01241       
01242       else if( !correct_inertia ) {
01243         omsg
01244           << "Indefinite update!\n" << onum_msg.str() << " >= wrong_inertia_tol = "
01245           << use_pivot_tols.wrong_inertia_tol;
01246         throw WrongInertiaUpdateException( omsg.str(), gamma );
01247       }
01248       else if( correct_inertia && gamma <= use_pivot_tols.warning_tol ) {
01249         omsg
01250           << "Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
01251           << " < " << onum_msg.str() << " <= warning_tol = "
01252           << use_pivot_tols.warning_tol;
01253         // Don't throw the exception till the very end!
01254       }
01255       else {
01256         TEST_FOR_EXCEPT(true); // Only local programming error?
01257       }
01258     }
01259     // u22 = sqrt(u22^2)
01260     u22 = std::sqrt(u22sqrabs);
01261   }
01262   else {
01263     factor_is_updated_ = false;
01264   }
01265   // Now augment the original
01266   if( maintain_original_ ) {
01267     //
01268     // M_new = [  M    t    ]
01269     //         [  t'  alpha ]
01270     //
01271     DVectorSlice M12 = MU_store_.row(M_l_r_+M_size_+1)(M_l_c_,M_l_c_+M_size_-1);
01272     if(t)
01273       M12 = *t;
01274     else
01275       M12 = 0.0;
01276     MU_store_(M_l_r_+M_size_+1,M_l_c_+M_size_) = alpha;
01277   }
01278   ++M_size_; // Enlarge the matrix by one
01279   is_diagonal_ = false;
01280   if( throw_exception )
01281     throw WarnNearSingularUpdateException(omsg.str(),gamma);
01282 }
01283 
01284 void MatrixSymPosDefCholFactor::delete_update(
01285   size_type        jd
01286   ,bool            force_refactorization
01287   ,EEigenValType   drop_eigen_val
01288   ,PivotTolerances pivot_tols
01289   )
01290 {
01291 #ifdef PROFILE_HACK_ENABLED
01292   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::delete_udpate(...)" );
01293 #endif
01294   typedef MatrixSymAddDelUpdateable MSADU;
01295 
01296   TEST_FOR_EXCEPTION(
01297     jd < 1 || M_size_ < jd, std::out_of_range
01298     ,"MatrixSymPosDefCholFactor::delete_update(jd,...): "
01299     "Error, the indice jd must be 1 <= jd <= rows()" );
01300 
01301   TEST_FOR_EXCEPT( !( drop_eigen_val == MSADU::EIGEN_VAL_UNKNOWN
01302     || (scale_ > 0.0 && drop_eigen_val == MSADU::EIGEN_VAL_POS)
01303     || (scale_ < 0.0 && drop_eigen_val == MSADU::EIGEN_VAL_NEG)
01304     ) );
01305 
01306   if( maintain_original_ ) {
01307     //
01308     // Here we have the lower portion of M partitioned as:
01309     //
01310     //  1 |\
01311     //    |  \
01312     //    |    \
01313     //    | M11  \
01314     //    |________\ _
01315     // jd |_________|_|
01316     //    |         | |\
01317     //    |         | |  \
01318     //    |         | |    \
01319     //    |   M31   | | M33  \
01320     //  n |         | |        \
01321     //    ----------------------
01322     //    1         jd         n
01323     //
01324     // We want to move M31 up one row and M33 up and to the left.
01325     //
01326     DenseLinAlgPack::delete_row_col( jd, &DenseLinAlgPack::nonconst_tri_ele( M().gms(), BLAS_Cpp::lower ) );
01327   }
01328   if( maintain_factor_ ) {
01329     //
01330     // Here we have U partitioned as:
01331     //
01332     // 1              jd      n        
01333     // ------------------------- 1
01334     // \             | |       |
01335     //   \    U11    | |  U13  |
01336     //     \         | |       |
01337     //       \  u12->| |  u23' |
01338     //         \     | |   |   |
01339     //           \   | |  \./  |
01340     //             \ |_|_______|
01341     //        u22 -> |_|_______| jd  
01342     //                 \       |
01343     //                   \ U33 |
01344     //                     \   |
01345     //                       \ | n
01346     //
01347     // To preform the update we need to first update U33 as:
01348     //   U33'*U33 + u23'*u23 ==> U33'*U33
01349     // Then we need to move U12 (one column at a time) to the
01350     // left one column and overwrite u12.
01351     // Then we will move the updated U33 (one column at a time)
01352     // up and to the left one position to overwrite u22 and
01353     // the left part of u23'.  We then decrease the dimension
01354     // of U by one and we are finished updating the factorization.
01355     //
01356     // See RAB notes from 1/21/99 and 1/26/99 for details on this update.
01357     //
01358     const size_type n = M_size_;
01359     // Resize workspace if it has not been done yet.
01360     size_type work_size = 3 * max_size_;
01361     if(work_.dim() < work_size)
01362       work_.resize(work_size);
01363     // Update the factors
01364     {
01365       DMatrixSlice U = this->U().gms();
01366       // Update U33 where it sits.
01367       if(jd < n) {
01368         size_type _size = n-jd; // Set storage for u23, c and s
01369         Range1D rng(1,_size);
01370         DVectorSlice
01371           u23 = work_(rng),
01372           c   = work_(rng+_size),
01373           s   = work_(rng+2*_size);
01374         Range1D U_rng(jd+1,n);  // Set U33 and u23
01375         DMatrixSlice U33 = U(U_rng,U_rng);
01376         u23 = U.row(jd)(U_rng);
01377         // Update U33
01378         value_type dummy;
01379         FORTRAN_FUNC_CALL_UL( DCHUD, dchud ) (
01380           U33.col_ptr(1), U33.max_rows()
01381           ,U_rng.size(), u23.start_ptr(), &dummy, 1, 0, &dummy
01382           ,&dummy, c.start_ptr(), s.start_ptr() );
01383       }
01384       // Move U13 and U33 to delete row and column jd
01385       DenseLinAlgPack::delete_row_col( jd, &DenseLinAlgPack::nonconst_tri_ele( U, BLAS_Cpp::upper ) );
01386     }
01387   }
01388   else {
01389     factor_is_updated_ = false;
01390   }
01391   // Strink the size of M and U
01392   --M_size_;
01393 }
01394 
01395 // Overridden from Serializable
01396 
01397 // ToDo: Refactor this code and create an external utility matrix
01398 // serialization class that will convert from matrix type to matrix
01399 // type.
01400 
01401 void MatrixSymPosDefCholFactor::serialize( std::ostream &out ) const
01402 {
01403   // Write key words on top line
01404   out << build_serialization_string() << std::endl;
01405   // Set the precision (very important!)
01406   out.precision(std::numeric_limits<value_type>::digits10+4);
01407   // Write the dimmension
01408   out << M_size_ << std::endl;
01409   if(M_size_) {
01410     // Write the matrix values
01411     if( maintain_original_ ) {
01412       const DMatrixSliceSym M = this->M();
01413       write_matrix( M.gms(), M.uplo(), out );
01414     }
01415     else {
01416       const DMatrixSliceTri U = this->U();
01417       write_matrix( U.gms(), U.uplo(), out );
01418     }
01419   }
01420   // ToDo: You need to write both M and U if both are computed!
01421 }
01422 
01423 void MatrixSymPosDefCholFactor::unserialize( std::istream &in )
01424 {
01425   // Get the keywords for the matrix type
01426   std::string keywords;
01427   std::getline( in, keywords, '\n' );
01428   // For now make sure the types are exactly the same!
01429   const std::string this_keywords = build_serialization_string();
01430   TEST_FOR_EXCEPTION(
01431     this_keywords != keywords, std::logic_error
01432     ,"MatrixSymPosDefCholFactor::unserialize(...): Error, the matrix type being read in from file of "
01433     "\'"<<keywords<<"\' does not equal the type expected of \'"<<this_keywords<<"\'!"
01434     );
01435   // Read in the dimension of the matrix
01436   in >> M_size_;
01437   TEST_FOR_EXCEPTION(
01438     M_size_ < 0, std::logic_error
01439     ,"MatrixSymPosDefCholFactor::unserialize(...): Error, read in a size of M_size = "<<M_size_<<" < 0!"
01440     );
01441   allocate_storage(M_size_);
01442   // Read in the matrix into storage
01443   if(maintain_original_) {
01444     M_l_r_ = M_l_c_ = 1;
01445     DMatrixSliceSym M = this->M();
01446     read_matrix( in, M.uplo(), &M.gms() );
01447   }
01448   else {
01449     U_l_r_ = U_l_c_ = 1;
01450     DMatrixSliceTri U = this->U();
01451     read_matrix( in, U.uplo(), &U.gms() );
01452   }
01453 }
01454 
01455 // Private member functions
01456 
01457 void MatrixSymPosDefCholFactor::assert_storage() const
01458 {
01459   TEST_FOR_EXCEPT( !(  MU_store_.rows()  ) );
01460 }
01461 
01462 void MatrixSymPosDefCholFactor::allocate_storage(size_type max_size) const
01463 {
01464   namespace rcp = MemMngPack;
01465   if( allocates_storage_ && MU_store_.rows() < max_size + 1 ) {
01466     // We have the right to allocate storage so lets just do it.
01467     Teuchos::RCP<DMatrix>
01468       MU_store = Teuchos::rcp(new DMatrix( max_size + 1, max_size + 1 ));
01469     typedef MemMngPack::ReleaseResource_ref_count_ptr<DMatrix> ptr_t;
01470     const_cast<MatrixSymPosDefCholFactor*>(this)->release_resource_ptr_ = Teuchos::rcp(new ptr_t(MU_store));
01471     const_cast<MatrixSymPosDefCholFactor*>(this)->MU_store_.bind( (*MU_store)() );
01472     const_cast<MatrixSymPosDefCholFactor*>(this)->max_size_ = max_size;
01473   }
01474   else {
01475     TEST_FOR_EXCEPT( !(  MU_store_.rows()  >= max_size + 1  ) );
01476   }
01477 }
01478 
01479 void MatrixSymPosDefCholFactor::assert_initialized() const
01480 {
01481   TEST_FOR_EXCEPT( !(  M_size_  ) );
01482 }
01483 
01484 void MatrixSymPosDefCholFactor::resize_and_zero_off_diagonal(size_type n, value_type scale)
01485 {
01486   using DenseLinAlgPack::nonconst_tri_ele;
01487   using DenseLinAlgPack::assign;
01488 
01489   TEST_FOR_EXCEPT( !(  n <= my_min( MU_store_.rows(), MU_store_.cols() ) - 1  ) );
01490 
01491   // Resize the views
01492   set_view(
01493     n, scale, maintain_original_, 1, 1, maintain_factor_
01494     ,U_l_r_ ? 1 : 0, U_l_r_ ? 1 : 0 );
01495   if( maintain_original_ ) {
01496     if(!is_diagonal_ && n > 1)
01497       assign( &nonconst_tri_ele( M().gms()(2,n,1,n-1), BLAS_Cpp::lower ), 0.0 );
01498   }
01499   if( U_l_r_ ) {
01500     if(!is_diagonal_ && n > 1)
01501       assign( &nonconst_tri_ele( U().gms()(1,n-1,2,n), BLAS_Cpp::upper ), 0.0 );
01502   }
01503 }
01504 
01505 void MatrixSymPosDefCholFactor::update_factorization() const
01506 {
01507   if( factor_is_updated_ ) return; // The factor should already be updated.
01508   TEST_FOR_EXCEPTION(
01509     U_l_r_ == 0, std::logic_error
01510     ,"MatrixSymPosDefCholFactor::update_factorization() : "
01511     "Error, U_l_r == 0 was set in MatrixSymPosDefCholFactor::set_view(...) "
01512     "and therefore we can not update the factorization in the provided storage space." );
01513   // Update the factorization from scratch!
01514 #ifdef PROFILE_HACK_ENABLED
01515   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::update_factorization(...) ... update" );
01516 #endif
01517   MatrixSymPosDefCholFactor
01518     *nc_this = const_cast<MatrixSymPosDefCholFactor*>(this);
01519   DMatrixSliceTriEle U = DenseLinAlgPack::nonconst_tri_ele( nc_this->U().gms(), BLAS_Cpp::upper );
01520   DenseLinAlgPack::assign( &U, DenseLinAlgPack::tri_ele( M().gms(), BLAS_Cpp::lower ) );  // Copy in the original
01521   {
01522 #ifdef PROFILE_HACK_ENABLED
01523     ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::update_factorization(...) ... potrf" );
01524 #endif
01525     DenseLinAlgLAPack::potrf( &U );
01526   }
01527   nc_this->factor_is_updated_ = true;
01528 }
01529 
01530 std::string MatrixSymPosDefCholFactor::build_serialization_string() const
01531 {
01532   std::string str = "SYMMETRIC POS_DEF";
01533   if( !maintain_original_ )
01534     str.append(" CHOL_FACTOR");
01535   if( maintain_original_ )
01536     str.append(" LOWER");
01537   else
01538     str.append(" UPPER");
01539   return str;
01540 }
01541 
01542 void MatrixSymPosDefCholFactor::write_matrix( const DMatrixSlice &Q, BLAS_Cpp::Uplo Q_uplo, std::ostream &out )
01543 {
01544   const int Q_dim = Q.rows();
01545   if( Q_uplo == BLAS_Cpp::lower ) {
01546     for( int i = 1; i <= Q_dim; ++i ) {
01547       for( int j = 1; j <= i; ++j ) {
01548         out << " " << Q(i,j);
01549       }
01550       out << std::endl;
01551     }
01552   }
01553   else {
01554     for( int i = 1; i <= Q_dim; ++i ) {
01555       for( int j = i; j <= Q_dim; ++j ) {
01556         out << " " << Q(i,j);
01557       }
01558       out << std::endl;
01559     }
01560   }
01561 }
01562 
01563 void MatrixSymPosDefCholFactor::read_matrix(  std::istream &in, BLAS_Cpp::Uplo Q_uplo, DMatrixSlice *Q_out )
01564 {
01565   DMatrixSlice &Q = *Q_out;
01566   const int Q_dim = Q.rows();
01567   if( Q_uplo == BLAS_Cpp::lower ) {
01568     for( int i = 1; i <= Q_dim; ++i ) {
01569       for( int j = 1; j <= i; ++j ) {
01570 #ifdef TEUCHOS_DEBUG
01571         TEST_FOR_EXCEPTION(
01572           in.eof(), std::logic_error
01573           ,"MatrixSymPosDefCholFactor::read_matrix(in,lower,Q_out): Error, not finished reading in matrix yet (i="<<i<<",j="<<j<<")!"
01574           );
01575 #endif
01576         in >> Q(i,j);
01577       }
01578     }
01579   }
01580   else {
01581     for( int i = 1; i <= Q_dim; ++i ) {
01582       for( int j = i; j <= Q_dim; ++j ) {
01583 #ifdef TEUCHOS_DEBUG
01584         TEST_FOR_EXCEPTION(
01585           in.eof(), std::logic_error
01586           ,"MatrixSymPosDefCholFactor::read_matrix(in,upper,Q_out): Error, not finished reading in matrix yet (i="<<i<<",j="<<j<<")!"
01587           );
01588 #endif
01589         in >> Q(i,j);
01590       }
01591     }
01592   }
01593 }
01594 
01595 } // end namespace AbstractLinAlgPack

Generated on Wed May 12 21:52:27 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7