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