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