AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_MatrixSymPosDefCholFactor.cpp
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 BLAS_Cpp::no_trans;
00392   using BLAS_Cpp::trans;
00393 #ifdef PROFILE_HACK_ENABLED
00394   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StMtV(...DVectorSlice...)" );
00395 #endif
00396   assert_initialized();
00397 
00398   DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), rows(), cols(), no_trans, x.dim() );
00399 
00400   if( maintain_original_ ) {
00401     //
00402     // M = symmetric
00403     //
00404     // y = b*y + a*M*x
00405     //
00406     DenseLinAlgPack::Vp_StMtV( y, a, M(), no_trans, x, b );
00407   }
00408   else {
00409     //
00410     // M = scale*U'*U
00411     //
00412     // y = b*y + a*op(M)*x
00413     //   = b*y = scale*a*U'*U*x
00414     //
00415     const DMatrixSliceTri
00416       U = this->U();
00417 
00418     if( b == 0.0 ) {
00419       // No temporary needed
00420       //
00421       // y = U*x
00422       DenseLinAlgPack::V_MtV( y, U, no_trans, x ); 
00423       // y = U'*y
00424       DenseLinAlgPack::V_MtV( y, U, trans, *y ); 
00425       // y *= scale*a
00426       if( a != 1.0 || scale_ != 1.0 )
00427         DenseLinAlgPack::Vt_S( y, scale_*a );
00428     }
00429     else {
00430       // We need a temporary
00431       DVector t;
00432       // t = U*x
00433       DenseLinAlgPack::V_MtV( &t, U, no_trans, x );
00434       // t = U'*t
00435       DenseLinAlgPack::V_MtV( &t(), U, trans, t() );
00436       // y *= b
00437       if(b != 1.0)
00438         DenseLinAlgPack::Vt_S( y, b );
00439       // y += scale*a*t
00440       DenseLinAlgPack::Vp_StV( y, scale_*a, t() );
00441     }
00442   }
00443 }
00444 
00445 void MatrixSymPosDefCholFactor::Vp_StMtV(
00446   DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00447   ,const SpVectorSlice& x, value_type b
00448   ) const
00449 {
00450   using BLAS_Cpp::no_trans;
00451   using BLAS_Cpp::trans;
00452 #ifdef PROFILE_HACK_ENABLED
00453   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StMtV(...SpVectorSlice...)" );
00454 #endif
00455   assert_initialized();
00456   if( maintain_original_ ) {
00457     const DMatrixSlice M = this->M().gms(); // This is lower triangular!
00458     const size_type n = M.rows();
00459     DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), n, n, no_trans, x.dim() );
00460     DenseLinAlgPack::Vt_S(y,b); // y = b*y
00461     //
00462     // Compute product column by column corresponding to x_itr->index() + x.offset()
00463     //
00464     // y += a * M * e(i) * x(i)
00465     //
00466     // [ y(1:i-1) ] += a * x(i) * [ ...  M(1:i-1,i) ... ]  stored as M(i,1:i-1) 
00467     // [ y(i:n)   ]               [ ...  M(i:n,i)   ... ]  stored as M(i:n,i)
00468     //
00469     for( SpVectorSlice::const_iterator x_itr = x.begin(); x_itr != x.end(); ++x_itr ) {
00470       const size_type i = x_itr->index() + x.offset();
00471       if( i > 1 )
00472         DenseLinAlgPack::Vp_StV( &(*y)(1,i-1), a * x_itr->value(), M.row(i)(1,i-1) );
00473       DenseLinAlgPack::Vp_StV( &(*y)(i,n), a * x_itr->value(), M.col(i)(i,n) );
00474     }
00475   }
00476   else {
00477     MatrixOpSerial::Vp_StMtV(y,a,M_trans,x,b); // ToDo: Specialize when needed!
00478   }
00479 }
00480 
00481 void MatrixSymPosDefCholFactor::Vp_StPtMtV(
00482   DVectorSlice* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00483   , BLAS_Cpp::Transp H_trans, const DVectorSlice& x, value_type b
00484   ) const
00485 {
00486 #ifdef PROFILE_HACK_ENABLED
00487   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StPtMtV(...DVectorSlice...)" );
00488 #endif
00489   assert_initialized();
00490   MatrixOpSerial::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // ToDo: Specialize when needed!
00491 }
00492 
00493 void MatrixSymPosDefCholFactor::Vp_StPtMtV(
00494   DVectorSlice* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00495   , BLAS_Cpp::Transp H_trans, const SpVectorSlice& x, value_type b
00496   ) const
00497 {
00498 #ifdef PROFILE_HACK_ENABLED
00499   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StPtMtV(...SpVectorSlice...)" );
00500 #endif
00501   assert_initialized();
00502   if( maintain_original_ ) {
00503     DenseLinAlgPack::Vt_S(y,b); // y = b*y
00504     const DMatrixSlice M = this->M().gms(); // This is lower triangular!
00505     // Compute product column by corresponding to x_itr->index() + x.offset()
00506     /*
00507     if( P.is_identity() ) {
00508       TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement
00509     }
00510     else {
00511       for( SpVectorSlice::const_iterator x_itr = x.begin(); x_itr != x.end(); ++x_itr ) {
00512         const size_type i = x_itr->index() + x.offset();
00513         
00514 
00515       }
00516     }
00517     */
00518     MatrixOpSerial::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // ToDo: Specialize when needed!
00519   }
00520   else {
00521     MatrixOpSerial::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // ToDo: Specialize when needed!
00522   }
00523 }
00524 
00525 // Overridden from MatrixSymOpSerial
00526 
00527 void MatrixSymPosDefCholFactor::Mp_StPtMtP(
00528   DMatrixSliceSym* S, value_type a
00529   , EMatRhsPlaceHolder dummy_place_holder
00530   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00531   , value_type b ) const
00532 {
00533   using BLAS_Cpp::no_trans;
00534   using BLAS_Cpp::trans;
00535 #ifdef PROFILE_HACK_ENABLED
00536   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Mp_StPtMtP(...)" );
00537 #endif
00538   assert_initialized();
00539   if( !maintain_original_ ) {
00540     MatrixSymOpSerial::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b);
00541   }
00542   else {
00543     MatrixSymOpSerial::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b); // ToDo: Override when needed!
00544   }
00545 }
00546 
00547 // Overridden from MatrixNonsingSerial
00548 
00549 void MatrixSymPosDefCholFactor::V_InvMtV(
00550   DVectorSlice* y, BLAS_Cpp::Transp M_trans, const DVectorSlice& x
00551   ) const
00552 {
00553   using BLAS_Cpp::no_trans;
00554   using BLAS_Cpp::trans;
00555 #ifdef PROFILE_HACK_ENABLED
00556   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::V_InvMtV(...DVectorSlice...)" );
00557 #endif
00558   assert_initialized();
00559 
00560   //
00561   // M = scale*U'*U
00562   //
00563   // y = inv(op(M))*x
00564   // =>
00565   // op(M)*y = x
00566   // =>
00567   // scale*U'*U*y = x
00568   // =>
00569   // y = (1/scale)*inv(U)*inv(U')*x
00570   //
00571   update_factorization();
00572   const DMatrixSliceTri
00573     U = this->U();
00574   DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), U.rows(), U.cols(), no_trans, x.dim() );
00575   // y = inv(U')*x
00576   DenseLinAlgPack::V_InvMtV( y, U, trans, x );
00577   // y = inv(U)*y
00578   DenseLinAlgPack::V_InvMtV( y, U, no_trans, *y );
00579   // y *= (1/scale)
00580   if( scale_ != 1.0 )
00581     DenseLinAlgPack::Vt_S( y, 1.0/scale_ );
00582 }
00583 
00584 void MatrixSymPosDefCholFactor::V_InvMtV(
00585   DVectorSlice* y, BLAS_Cpp::Transp M_trans, const SpVectorSlice& x
00586   ) const
00587 {
00588 #ifdef PROFILE_HACK_ENABLED
00589   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::V_InvMtV(...SpVectorSlice...)" );
00590 #endif
00591   assert_initialized();
00592   MatrixNonsingSerial::V_InvMtV(y,M_trans,x);
00593 }
00594 
00595 // Overridden from MatrixSymNonsingSerial
00596 
00597 void MatrixSymPosDefCholFactor::M_StMtInvMtM(
00598     DMatrixSliceSym* S, value_type a, const MatrixOpSerial& B
00599   , BLAS_Cpp::Transp B_trans, EMatrixDummyArg dummy_arg
00600   ) const
00601 {
00602 #ifdef PROFILE_HACK_ENABLED
00603   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::M_StMtInvMtM(...)" );
00604 #endif
00605 
00606 //  // Uncomment to use the defalut implementation (for debugging)
00607 //  MatrixSymFactorized::M_StMtInvMtM(S,a,B,B_trans,dummy_arg); return;
00608 
00609   using BLAS_Cpp::trans;
00610   using BLAS_Cpp::no_trans;
00611   using BLAS_Cpp::trans_not;
00612   using BLAS_Cpp::upper;
00613   using BLAS_Cpp::nonunit;
00614   using DenseLinAlgPack::tri;
00615   using DenseLinAlgPack::syrk;
00616   using DenseLinAlgPack::M_StInvMtM;
00617   using LinAlgOpPack::assign;
00618 
00619   assert_initialized();
00620   update_factorization();
00621 
00622   //
00623   // S = a * op(B) * inv(M) * op(B)'
00624   // 
00625   // M = scale*U'*U
00626   // =>
00627   // inv(M) = scale*inv(U'*U) = scale*inv(U)*inv(U')
00628   // =>
00629   // S = scale*a * op(B) * inv(U) * inv(U') * op(B)'
00630   // 
00631   // T = op(B)'
00632   // 
00633   // T = inv(U') * T (inplace with BLAS)
00634   // 
00635   // S = scale*a * T' * T
00636   // 
00637   DenseLinAlgPack::MtM_assert_sizes( 
00638     rows(), cols(), no_trans
00639     ,B.rows(), B.cols(), trans_not(B_trans)
00640     );
00641   DenseLinAlgPack::Mp_MtM_assert_sizes(
00642     S->rows(), S->cols(), no_trans
00643     ,B.rows(), B.cols(), B_trans
00644     ,B.rows(), B.cols(), trans_not(B_trans)
00645     );
00646   // T = op(B)'
00647   DMatrix T;
00648   assign( &T, B, trans_not(B_trans) );
00649   // T = inv(U') * T (inplace with BLAS)
00650   M_StInvMtM( &T(), 1.0, this->U(), trans, T(), no_trans );
00651   // S = scale*a * T' * T
00652   syrk( trans, scale_*a, T(), 0.0, S );
00653 }
00654 
00655 // Overridden from MatrixSymDenseInitialize
00656 
00657 void MatrixSymPosDefCholFactor::initialize( const DMatrixSliceSym& M )
00658 {
00659   // Initialize without knowing the inertia but is must be p.d.
00660   this->initialize(
00661     M, M.rows(), maintain_factor_
00662     , MatrixSymAddDelUpdateable::Inertia()
00663     , MatrixSymAddDelUpdateable::PivotTolerances()
00664     );
00665 }
00666 
00667 // Overridden from MatrixSymOpGetGMSSym
00668 
00669 const DenseLinAlgPack::DMatrixSliceSym MatrixSymPosDefCholFactor::get_sym_gms_view() const
00670 {
00671   TEUCHOS_TEST_FOR_EXCEPTION(
00672     !maintain_original_, std::logic_error
00673     ,"MatrixSymPosDefCholFactor::get_sym_gms_view(): Error, maintain_original must be "
00674     "true in order to call this method!" );
00675   return this->M();
00676 }
00677 
00678 void MatrixSymPosDefCholFactor::free_sym_gms_view(const DenseLinAlgPack::DMatrixSliceSym* sym_gms_view) const
00679 {
00680   TEUCHOS_TEST_FOR_EXCEPTION(
00681     !maintain_original_, std::logic_error
00682     ,"MatrixSymPosDefCholFactor::free_sym_gms_view(...): Error, maintain_original must be "
00683     "true in order to call this method!" );
00684   // Nothing todo
00685 }
00686 
00687 // Overridden from MatrixSymOpGetGMSSymMutable
00688 
00689 DenseLinAlgPack::DMatrixSliceSym MatrixSymPosDefCholFactor::get_sym_gms_view()
00690 {
00691   TEUCHOS_TEST_FOR_EXCEPTION(
00692     !maintain_original_, std::logic_error
00693     ,"MatrixSymPosDefCholFactor::get_sym_gms_view(): Error, maintain_original must be "
00694     "true in order to call this method!" );
00695   return this->M();
00696 }
00697 
00698 void MatrixSymPosDefCholFactor::commit_sym_gms_view(DenseLinAlgPack::DMatrixSliceSym* sym_gms_view)
00699 {
00700   TEUCHOS_TEST_FOR_EXCEPTION(
00701     !maintain_original_, std::logic_error
00702     ,"MatrixSymPosDefCholFactor::commit_sym_gms_view(...): Error, maintain_original must be "
00703     "true in order to call this method!" );
00704   this->initialize(*sym_gms_view);
00705 }
00706 
00707 // Overridden from MatrixExtractInvCholFactor
00708 
00709 void MatrixSymPosDefCholFactor::extract_inv_chol( DMatrixSliceTriEle* InvChol ) const
00710 {
00711   assert_initialized();
00712   update_factorization();
00713   //
00714   // The matrix is represented as the upper cholesky factor:
00715   //   M = scale * U' * U
00716   //
00717   //   inv(M) = inv(scale*U*U') = (1/sqrt(scale))*inv(U)*(1/sqrt(scale))*inv(U')
00718   //          = UInv * UInv'
00719   // =>
00720   //   UInv = (1/sqrt(scale))*inv(U)
00721   //
00722    // Here scale > 0 or an exception will be thrown.
00723   //
00724   // Compute the inverse cholesky factor as:
00725   //
00726   // Upper cholesky:
00727   //    sqrt(scale) * U * UInv = I => InvChol = UInv = (1/sqrt(scale))*inv(U) * I
00728   //    
00729   // Lower cholesky:
00730   //    sqrt(scale) * L * LInv = I => InvChol = LInv = (1/sqrt(scale))*inv(U) * inv(U') * I
00731   //
00732   TEUCHOS_TEST_FOR_EXCEPTION(
00733     scale_ < 0.0, std::logic_error
00734     ,"MatrixSymPosDefCholFactor::extract_inv_chol(...) : "
00735     "Error, we can not compute the inverse cholesky factor "
00736     "af a negative definite matrix." );
00737   DenseLinAlgPack::assign( &InvChol->gms(), 0.0 );  // Set InvChol to identity first.
00738   InvChol->gms().diag() = 1.0;
00739   DenseLinAlgPack::M_StInvMtM(                      // Comput InvChol using Level-3 BLAS
00740       &InvChol->gms(), 1.0 / std::sqrt(scale_), U()
00741     , InvChol->uplo() == BLAS_Cpp::upper ? BLAS_Cpp::no_trans : BLAS_Cpp::trans
00742     , InvChol->gms(), BLAS_Cpp::no_trans );
00743 }
00744 
00745 // Overridden from MatrixSymSecantUpdateble
00746 
00747 void MatrixSymPosDefCholFactor::init_identity( const VectorSpace& space_diag, value_type alpha )
00748 {
00749   const size_type n = space_diag.dim();
00750   allocate_storage( max_size_ ? max_size_ : n );
00751   //
00752   // B = alpha*I = = alpha*I*I = scale*U'*U
00753   // =>
00754   // U = sqrt(|alpha|) * I
00755   //
00756   // Here we will set scale = sign(alpha)
00757   //
00758   const value_type scale = alpha > 0.0 ? +1.0: -1.0; // Explicitly set the scale
00759   resize_and_zero_off_diagonal(n,scale);
00760   if( maintain_original_ ) {
00761     M().gms().diag()(1,n) = alpha;
00762   }
00763   if( maintain_factor_ ) {
00764     U().gms().diag()(1,n) = std::sqrt(std::fabs(alpha));
00765     factor_is_updated_ = true;
00766   }
00767   is_diagonal_ = true;
00768 }
00769 
00770 void MatrixSymPosDefCholFactor::init_diagonal( const Vector& diag_in )
00771 {
00772   VectorDenseEncap diag_encap(diag_in);
00773   const DVectorSlice diag = diag_encap(); // When diag_encap is destroyed, bye-bye view!
00774 
00775   allocate_storage( max_size_ ? max_size_ : diag.dim() );
00776   //
00777   // M = scale * U' * U = scale * (1/scale)*diag^(1/2) * (1/scale)*diag^(1/2)
00778   //
00779   // Here we will set scale = sign(diag(1)) and validate the rest
00780   //
00781   if( diag.dim() == 0 ) {
00782     M_size_ = 0;
00783     return; // We are unsizing this thing
00784   }
00785   const value_type scale = diag(1) > 0.0 ? +1.0: -1.0;
00786   resize_and_zero_off_diagonal(diag.dim(),scale);
00787   if( maintain_original_ ) {
00788     M().gms().diag() = diag;
00789     // ToDo: validate that scale*diag > 0
00790   }
00791   if( maintain_factor_ ) {
00792     DVectorSlice U_diag = U().gms().diag();
00793     U_diag = diag;
00794     if( scale_ != 1.0 )
00795       DenseLinAlgPack::Vt_S( &U_diag, 1.0/scale_ );
00796     DenseLinAlgPack::sqrt( &U_diag, U_diag );
00797     DenseLinAlgPack::assert_print_nan_inf( U_diag, "(1/scale)*diag", true, NULL );
00798     factor_is_updated_ = true;
00799   }
00800   is_diagonal_ = true;
00801 }
00802 
00803 void MatrixSymPosDefCholFactor::secant_update(
00804   VectorMutable* s_in, VectorMutable* y_in, VectorMutable* Bs_in
00805   )
00806 {
00807   using BLAS_Cpp::no_trans;
00808   using BLAS_Cpp::trans;
00809   using DenseLinAlgPack::dot;
00810   using DenseLinAlgPack::norm_2;
00811   using DenseLinAlgPack::norm_inf;
00812   namespace rcp = MemMngPack;
00813 
00814   assert_initialized();
00815 
00816   // Validate the input
00817   TEUCHOS_TEST_FOR_EXCEPT( !(  s_in && y_in  ) );
00818   DenseLinAlgPack::Vp_MtV_assert_sizes( y_in->dim(), M_size_, M_size_, no_trans, s_in->dim() );
00819 
00820   // Get the serial vectors
00821   VectorDenseMutableEncap s_encap(*s_in);
00822   VectorDenseMutableEncap y_encap(*y_in);
00823   VectorDenseMutableEncap Bs_encap( Bs_in ? *Bs_in : *y_in); // Must pass something on
00824   DVectorSlice
00825     *s  = &s_encap(),   // When s_encap, y_encap and Bs_encap are destroyed
00826     *y  = &y_encap(),   // these views go bye-bye!
00827     *Bs = ( Bs_in ? &Bs_encap() : NULL );
00828 
00829   // Check skipping the BFGS update
00830   const value_type
00831     sTy       = dot(*s,*y),
00832     sTy_scale = sTy*scale_;
00833   std::ostringstream omsg;
00834   if( !BFGS_sTy_suff_p_d(
00835       VectorMutableDense((*s)(),Teuchos::null)
00836       ,VectorMutableDense((*y)(),Teuchos::null)
00837       ,&sTy_scale,&omsg,"\nMatrixSymPosDefCholFactor::secant_update(...)"
00838       )
00839     )
00840   {
00841     throw UpdateSkippedException( omsg.str() ); 
00842   }
00843   // Compute Bs if it was not passed in
00844   DVector Bs_store;
00845   DVectorSlice Bs_view;
00846   if( !Bs ) {
00847     LinAlgOpPack::V_MtV( &Bs_store, *this, no_trans, *s );
00848     Bs_view.bind( Bs_store() );
00849     Bs = &Bs_view;
00850   }
00851   // Check that s'*Bs is positive and if not then throw exception
00852   const value_type sTBs = dot(*s,*Bs);
00853   TEUCHOS_TEST_FOR_EXCEPTION(
00854     scale_*sTBs <= 0.0 && scale_ > 0.0, UpdateFailedException
00855     ,"MatrixSymPosDefCholFactor::secant_update(...) : "
00856     "Error, B can't be positive definite if s'*Bs <= 0.0" );
00857   TEUCHOS_TEST_FOR_EXCEPTION(
00858     scale_*sTBs <= 0.0 && scale_ <= 0.0, UpdateFailedException
00859     ,"MatrixSymPosDefCholFactor::secant_update(...) : "
00860     "Error, B can't be negative definite if s'*Bs >= 0.0" );
00861   if( maintain_original_ ) {
00862     //
00863     // Compute the BFGS update of the original, nonfactored matrix
00864     //
00865     // Just preform two symmetric updates.
00866     //
00867     // B = B + (-1/s'*Bs) * Bs*Bs' + (1/s'*y) * y*y'
00868     //
00869     DMatrixSliceSym M = this->M();
00870     DenseLinAlgPack::syr( -1.0/sTBs, *Bs, &M );
00871     DenseLinAlgPack::syr( 1.0/sTy, *y, &M );
00872   }
00873   if( maintain_factor_ ) {
00874     //
00875     // Compute the BFGS update for the cholesky factor
00876     //
00877     // If this implementation is based on the one in Section 9.2, page 198-201 of:
00878     //
00879     // Dennis, J.E., R.B. Schnabel
00880     // "Numerical Methods for Unconstrained Optimization"
00881     //
00882     // Given that we have B = scale*U'*U and the BFGS update:
00883     //
00884     // B_new = B + y*y'/(y'*s) - (B*s)*(B*s)'/(s'*B*s)
00885     //
00886     // We can rewrite it in the following form:
00887     //
00888     // B_new = scale*(U' + a*u*v')*(U + a*v*u') = scale*U_new'*U_new
00889     //
00890     // where:
00891     //     v = sqrt(y'*s/(s'*B*s))*U*s
00892     //     u = y - U'*v
00893     //     a = 1/(v'*v)
00894     //
00895     DMatrixSliceTri U = this->U();
00896     // v = sqrt(y'*s/(s'*B*s))*U*s
00897     DVectorSlice v = *s; // Reuse s as storage for v
00898     DenseLinAlgPack::V_MtV( &v, U, no_trans, v ); // Direct call to xSYMV(...)
00899     DenseLinAlgPack::Vt_S( &v, std::sqrt( sTy / sTBs ) );
00900     // u = (y - U'*v)
00901     DVectorSlice u = *y;  // Reuse y as storage for u
00902     DenseLinAlgPack::Vp_StMtV( &u, -1.0, U, trans, v );
00903     // a = 1/(v'*v)
00904     const value_type  a = 1.0/dot(v,v);
00905     // Perform Givens rotations to make Q*(U' + a*u*v') -> U_new upper triangular:
00906     //
00907     // B_new = scale*(U' + a*u*v')*Q'*Q*(U + a*v*u') = scale*U_new'*U_new
00908     rank_2_chol_update(
00909       a, &v, u, v.dim() > 1 ? &U.gms().diag(-1) : NULL
00910       , &DMatrixSliceTriEle(U.gms(),BLAS_Cpp::upper), no_trans );
00911   }
00912   else {
00913     factor_is_updated_ = false;
00914   }
00915   is_diagonal_ = false;
00916 }
00917 
00918 // Overridden from MatrixSymAddDelUpdateble
00919 
00920 void MatrixSymPosDefCholFactor::initialize(
00921   value_type         alpha
00922   ,size_type         max_size
00923   )
00924 {
00925 #ifdef PROFILE_HACK_ENABLED
00926   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::initialize(alpha,max_size)" );
00927 #endif
00928   allocate_storage(max_size);
00929 
00930   if( alpha == 0.0 ) 
00931     throw SingularUpdateException(
00932       "MatrixSymPosDefCholFactor::initialize(...): "
00933       "Error, alpha == 0.0, matrix is singular.", 0.0 );
00934 
00935   // Resize the view
00936   if( maintain_original_ ) {
00937     M_l_r_ = 1;
00938     M_l_c_ = 1;
00939   }
00940   if( U_l_r_ ) {
00941     U_l_r_ = 1;
00942     U_l_c_ = 1;
00943   }
00944   M_size_ = 1;
00945   max_size_ = my_min( MU_store_.rows(), MU_store_.cols() ) - 1;
00946   scale_ = alpha > 0.0 ? +1.0 : -1.0;
00947   // Update the matrix
00948   if( maintain_original_ ) {
00949     M().gms()(1,1) = alpha;
00950   }
00951   if( U_l_r_ ) {
00952     U().gms()(1,1) = std::sqrt( scale_ * alpha );
00953     factor_is_updated_ = true;
00954   }
00955   is_diagonal_ = false;
00956 }
00957 
00958 void MatrixSymPosDefCholFactor::initialize(
00959   const DMatrixSliceSym      &A
00960   ,size_type         max_size
00961   ,bool              force_factorization
00962   ,Inertia           inertia
00963   ,PivotTolerances   pivot_tols
00964   )
00965 {
00966 #ifdef PROFILE_HACK_ENABLED
00967   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::initialize(A,max_size...)" );
00968 #endif
00969   typedef MatrixSymAddDelUpdateable::Inertia   Inertia;
00970 
00971   allocate_storage(max_size);
00972 
00973   const size_type
00974     n = A.rows();
00975 
00976   // Validate proper usage of inertia parameter
00977   TEUCHOS_TEST_FOR_EXCEPT( !(  inertia.zero_eigens == Inertia::UNKNOWN || inertia.zero_eigens == 0  ) );
00978   TEUCHOS_TEST_FOR_EXCEPT( !(  (inertia.neg_eigens == Inertia::UNKNOWN && inertia.pos_eigens == Inertia::UNKNOWN ) )
00979       || ( inertia.neg_eigens == n && inertia.pos_eigens == 0 )
00980       || ( inertia.neg_eigens == 0 && inertia.pos_eigens == n )
00981     );
00982 
00983   // We can infer if the matrix is p.d. or n.d. by the sign of the diagonal
00984   // elements.  If a matrix is s.p.d. (s.n.d) then A(i,i) > 0 (A(i,i) < 0)
00985   // for all i = 1...n so we will just check the first i.
00986   const value_type
00987     a_11 = A.gms()(1,1);
00988   const int sign_a_11 = ( a_11 == 0.0 ? 0 : ( a_11 > 0.0 ? +1 : -1 ) );
00989   if( sign_a_11 == 0.0 )
00990     std::logic_error(
00991       "MatrixSymPosDefCholFactor::initialize(...) : "
00992       "Error, A can not be positive definite or negative definete "
00993       "if A(1,1) == 0.0" );
00994   if( inertia.pos_eigens == n && sign_a_11 < 0 )
00995     std::logic_error(
00996       "MatrixSymPosDefCholFactor::initialize(...) : "
00997       "Error, A can not be positive definite "
00998       "if A(1,1) < 0.0" );
00999   if( inertia.neg_eigens == n && sign_a_11 > 0 )
01000     std::logic_error(
01001       "MatrixSymPosDefCholFactor::initialize(...) : "
01002       "Error, A can not be negative definite "
01003       "if A(1,1) > 0.0" );
01004   // Now we have got the sign
01005   const value_type
01006     scale = (value_type)sign_a_11;
01007   // Setup the view
01008   set_view(
01009     n,scale,maintain_original_,1,1
01010     ,maintain_factor_, U_l_r_ ? 1 : 0, U_l_r_ ? 1 : 0
01011     );
01012   // Now set the matrix and update the factors
01013   if( maintain_original_ ) {
01014     // Set M = S
01015     DenseLinAlgPack::assign( 
01016       &DMatrixSliceTriEle( M().gms(), BLAS_Cpp::lower )
01017       ,DMatrixSliceTriEle( A.gms(), A.uplo() )
01018       );
01019   }
01020   if( maintain_factor_ || force_factorization ) {
01021     // Copy S into U for an inplace factorization.
01022     DMatrixSliceTriEle U_ele = DMatrixSliceTriEle( U().gms(), BLAS_Cpp::upper );
01023     DenseLinAlgPack::assign( &U_ele, DMatrixSliceTriEle( A.gms(), A.uplo() ) );
01024     if( sign_a_11 < 0 )
01025       DenseLinAlgPack::Mt_S( &U_ele, -1.0 );
01026     try {
01027       DenseLinAlgLAPack::potrf( &U_ele );
01028       factor_is_updated_ = true;
01029     }
01030     catch( const DenseLinAlgLAPack::FactorizationException& excpt ) {
01031       M_size_ = 0;  // set unsized
01032       throw SingularUpdateException( excpt.what(), 0.0 );
01033     }
01034     catch(...) {
01035       M_size_ = 0;
01036       throw;
01037     }
01038     // Validate that the tolerances are met and throw appropriate
01039     // exceptions.  We already know that the matrix is technically
01040     // p.d. or n.d.  Now we must determine gamma = (min(|diag|)/max(|diag|))^2
01041     value_type
01042       min_diag = std::numeric_limits<value_type>::max(),
01043       max_diag = 0.0;
01044     DVectorSlice::iterator
01045       U_itr = U_ele.gms().diag().begin(),
01046       U_end = U_ele.gms().diag().end();
01047     while( U_itr != U_end ) {
01048       const value_type U_abs = std::abs(*U_itr++);
01049       if(U_abs < min_diag) min_diag = U_abs;
01050       if(U_abs > max_diag) max_diag = U_abs;
01051     }
01052     const value_type gamma = (min_diag*min_diag)/(max_diag*max_diag);
01053     // Validate gamma
01054     PivotTolerances use_pivot_tols = pivot_tols_;
01055     if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
01056       use_pivot_tols.warning_tol = pivot_tols.warning_tol;
01057     if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
01058       use_pivot_tols.singular_tol = pivot_tols.singular_tol;
01059     if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
01060       use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
01061     const bool throw_exception = (
01062       gamma == 0.0
01063       || gamma <= use_pivot_tols.warning_tol
01064       || gamma <= use_pivot_tols.singular_tol
01065       );
01066     // Create message and throw exceptions
01067     std::ostringstream onum_msg, omsg;
01068     if(throw_exception) {
01069       onum_msg
01070         << "gamma = (min(|diag(U)(k)|)/|max(|diag(U)(k)|))^2 = (" << min_diag <<"/"
01071         << max_diag << ")^2 = " << gamma;
01072       omsg
01073         << "MatrixSymPosDefCholFactor::initialize(...): ";
01074       if( gamma <= use_pivot_tols.singular_tol ) {
01075         M_size_ = 0;  // The initialization failed!
01076         omsg
01077           << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
01078           << use_pivot_tols.singular_tol;
01079         throw SingularUpdateException( omsg.str(), gamma );
01080       }
01081       else if( gamma <= use_pivot_tols.warning_tol ) {
01082         omsg
01083           << "Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
01084           << " < " << onum_msg.str() << " <= warning_tol = "
01085           << use_pivot_tols.warning_tol;
01086         throw WarnNearSingularUpdateException( omsg.str(), gamma ); // The initialization still succeeded through
01087       }
01088       else {
01089         TEUCHOS_TEST_FOR_EXCEPT(true); // Only local programming error?
01090       }
01091     }
01092   }
01093   else {
01094     factor_is_updated_ = false; // The factor is not updated!
01095   }
01096 }
01097 
01098 size_type MatrixSymPosDefCholFactor::max_size() const
01099 {
01100   return max_size_;
01101 }
01102 
01103 MatrixSymAddDelUpdateable::Inertia
01104 MatrixSymPosDefCholFactor::inertia() const
01105 {
01106   typedef MatrixSymAddDelUpdateable MSADU;
01107   typedef MSADU::Inertia Inertia;
01108   return ( M_size_
01109        ? ( scale_ > 0.0
01110          ? Inertia(0,0,M_size_)
01111          : Inertia(M_size_,0,0) )
01112        : Inertia(0,0,0) );
01113 }
01114 
01115 void MatrixSymPosDefCholFactor::set_uninitialized()
01116 {
01117   M_size_ = 0;
01118 }
01119 
01120 void MatrixSymPosDefCholFactor::augment_update(
01121   const DVectorSlice  *t
01122   ,value_type        alpha
01123   ,bool              force_refactorization
01124   ,EEigenValType     add_eigen_val
01125   ,PivotTolerances   pivot_tols
01126   )
01127 {
01128 #ifdef PROFILE_HACK_ENABLED
01129   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::augment_udpate(...)" );
01130 #endif
01131   using DenseLinAlgPack::dot;
01132   using DenseLinAlgPack::norm_inf;
01133   typedef MatrixSymAddDelUpdateable  MSADU;
01134 
01135   assert_initialized();
01136 
01137   // Validate the input
01138   TEUCHOS_TEST_FOR_EXCEPTION(
01139     rows() >= max_size(), MaxSizeExceededException
01140     ,"MatrixSymPosDefCholFactor::augment_update(...) : "
01141     "Error, the maximum size would be exceeded." );
01142   TEUCHOS_TEST_FOR_EXCEPTION(
01143     t && t->dim() != M_size_, std::length_error
01144     ,"MatrixSymPosDefCholFactor::augment_update(...): "
01145     "Error, t.dim() must be equal to this->rows()." );
01146   if( !(add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN || add_eigen_val != MSADU::EIGEN_VAL_ZERO ) )
01147     throw SingularUpdateException(
01148       "MatrixSymPosDefCholFactor::augment_update(...): "
01149       "Error, the client has specified a singular update in add_eigen_val.", -1.0 );
01150   if( alpha == 0.0 )
01151     throw SingularUpdateException(
01152       "MatrixSymPosDefCholFactor::augment_update(...): "
01153       "Error, alpha == 0.0 so the matrix is not positive definite "
01154       "or negative definite.", -1.0 );
01155   if( scale_ > 0.0 && alpha < 0.0 )
01156     throw WrongInertiaUpdateException(
01157       "MatrixSymPosDefCholFactor::augment_update(...): "
01158       "Error, alpha < 0.0 so the matrix is not postivie 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 negative definite ", -1.0 );
01163   if( scale_ > 0.0 && !( add_eigen_val == MSADU::EIGEN_VAL_POS || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
01164     throw WrongInertiaUpdateException(
01165       "MatrixSymPosDefCholFactor::augment_update(...): "
01166       "Error, alpha > 0.0 but the user specified a non-positive definite new matrix.", -1.0 );
01167   if( scale_ < 0.0 && !( add_eigen_val == MSADU::EIGEN_VAL_NEG || 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   // First try to augment the factor to verify that the matrix is still p.d. or n.d.
01172   bool throw_exception = false; // If true then throw exception
01173   std::ostringstream omsg;      // Will be set if an exception has to be thrown.
01174   value_type gamma;             // ...
01175   if( maintain_factor_ ) {
01176     //
01177     // The update is:
01178     //
01179     // B_new = [  B     t   ]
01180     //         [  t'  alpha ]
01181     //
01182     //       = scale * [     U'*U         (1/scale)*t   ]
01183     //                 [ (1/scale)*t'   (1/scale)*alpha ]
01184     //
01185     // We seek to find a new cholesky factor of the form:
01186     //
01187     // U_new = [ U11  u12  ]
01188     //         [  0   u22  ]
01189     //
01190     // B_new = scale*U_new'*U_new
01191     //
01192     //       = scale * [ U11'   0  ] * [ U11   u12 ]
01193     //                 [ u12'  u22 ]   [  0    u22 ]
01194     //
01195     //       = scale * [  U11'*U11          U11'*u12  ]
01196     //                 [  u12'*U11   u12'*u12 + u22^2 ]
01197     //
01198     // From the above we see that:
01199     // =>
01200     //   U11 = U
01201     //   u12 = inv(U') * (1/scale) * t
01202     //   u22 = sqrt( (1/scale)*alpha - u12'*u12 );
01203     //
01204     // We must compute gamma relative to the LU factorization
01205     // so we must square ||U11.diag()||inf.
01206 
01207     // Get references to the storage for the to-be-updated parts for the new factor.
01208     DVectorSlice u12 = MU_store_.col(U_l_c_+M_size_+1)(U_l_r_,U_l_r_+M_size_-1);
01209     value_type &u22 = MU_store_(U_l_r_+M_size_,U_l_c_+M_size_+1);
01210     // u12 = inv(U') * (1/scale) * t
01211     if(t) {
01212       DenseLinAlgPack::V_InvMtV( &u12, U(), BLAS_Cpp::trans, *t );
01213       if( scale_ != 1.0 ) DenseLinAlgPack::Vt_S( &u12, 1.0/scale_ );
01214     }
01215     else {
01216       u12 = 0.0;
01217     }
01218     // u22^2 = (1/scale)*alpha - u12'*u12;
01219     const value_type
01220       u22sqr          = (1/scale_) * alpha - ( t ? dot( u12, u12 ) : 0.0 ),
01221       u22sqrabs       = std::abs(u22sqr),
01222       nrm_U_diag      = norm_inf(U().gms().diag()),
01223       sqr_nrm_U_diag  = nrm_U_diag * nrm_U_diag;
01224     // Calculate gamma in proper context
01225     gamma = u22sqrabs / sqr_nrm_U_diag;
01226     // Check gamma
01227     const bool
01228       correct_inertia = u22sqr > 0.0;
01229     PivotTolerances use_pivot_tols = pivot_tols_;
01230     if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
01231       use_pivot_tols.warning_tol = pivot_tols.warning_tol;
01232     if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
01233       use_pivot_tols.singular_tol = pivot_tols.singular_tol;
01234     if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
01235       use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
01236     throw_exception = (
01237       gamma == 0.0
01238       || correct_inertia && gamma <= use_pivot_tols.warning_tol
01239       || correct_inertia && gamma <= use_pivot_tols.singular_tol
01240       || !correct_inertia
01241       );
01242     // Create message and throw exceptions
01243     std::ostringstream onum_msg;
01244     if(throw_exception) {
01245       onum_msg
01246         << "gamma = |u22^2|/(||diag(U11)||inf)^2 = |" << u22sqr <<"|/("
01247         << nrm_U_diag << ")^2 = " << gamma;
01248       omsg
01249         << "MatrixSymPosDefCholFactor::augment_update(...): ";
01250       if( correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
01251         omsg
01252           << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
01253           << use_pivot_tols.singular_tol;
01254         throw SingularUpdateException( omsg.str(), gamma );
01255       }
01256       else if( !correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
01257         omsg
01258           << "Singular update!\n" << onum_msg.str() << " <= wrong_inertia_tol = "
01259           << use_pivot_tols.wrong_inertia_tol;
01260         throw SingularUpdateException( omsg.str(), gamma );
01261       }
01262       
01263       else if( !correct_inertia ) {
01264         omsg
01265           << "Indefinite update!\n" << onum_msg.str() << " >= wrong_inertia_tol = "
01266           << use_pivot_tols.wrong_inertia_tol;
01267         throw WrongInertiaUpdateException( omsg.str(), gamma );
01268       }
01269       else if( correct_inertia && gamma <= use_pivot_tols.warning_tol ) {
01270         omsg
01271           << "Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
01272           << " < " << onum_msg.str() << " <= warning_tol = "
01273           << use_pivot_tols.warning_tol;
01274         // Don't throw the exception till the very end!
01275       }
01276       else {
01277         TEUCHOS_TEST_FOR_EXCEPT(true); // Only local programming error?
01278       }
01279     }
01280     // u22 = sqrt(u22^2)
01281     u22 = std::sqrt(u22sqrabs);
01282   }
01283   else {
01284     factor_is_updated_ = false;
01285   }
01286   // Now augment the original
01287   if( maintain_original_ ) {
01288     //
01289     // M_new = [  M    t    ]
01290     //         [  t'  alpha ]
01291     //
01292     DVectorSlice M12 = MU_store_.row(M_l_r_+M_size_+1)(M_l_c_,M_l_c_+M_size_-1);
01293     if(t)
01294       M12 = *t;
01295     else
01296       M12 = 0.0;
01297     MU_store_(M_l_r_+M_size_+1,M_l_c_+M_size_) = alpha;
01298   }
01299   ++M_size_; // Enlarge the matrix by one
01300   is_diagonal_ = false;
01301   if( throw_exception )
01302     throw WarnNearSingularUpdateException(omsg.str(),gamma);
01303 }
01304 
01305 void MatrixSymPosDefCholFactor::delete_update(
01306   size_type        jd
01307   ,bool            force_refactorization
01308   ,EEigenValType   drop_eigen_val
01309   ,PivotTolerances pivot_tols
01310   )
01311 {
01312 #ifdef PROFILE_HACK_ENABLED
01313   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::delete_udpate(...)" );
01314 #endif
01315   typedef MatrixSymAddDelUpdateable MSADU;
01316 
01317   TEUCHOS_TEST_FOR_EXCEPTION(
01318     jd < 1 || M_size_ < jd, std::out_of_range
01319     ,"MatrixSymPosDefCholFactor::delete_update(jd,...): "
01320     "Error, the indice jd must be 1 <= jd <= rows()" );
01321 
01322   TEUCHOS_TEST_FOR_EXCEPT( !( drop_eigen_val == MSADU::EIGEN_VAL_UNKNOWN
01323     || (scale_ > 0.0 && drop_eigen_val == MSADU::EIGEN_VAL_POS)
01324     || (scale_ < 0.0 && drop_eigen_val == MSADU::EIGEN_VAL_NEG)
01325     ) );
01326 
01327   if( maintain_original_ ) {
01328     //
01329     // Here we have the lower portion of M partitioned as:
01330     //
01331     //  1 |\
01332     //    |  \
01333     //    |    \
01334     //    | M11  \
01335     //    |________\ _
01336     // jd |_________|_|
01337     //    |         | |\
01338     //    |         | |  \
01339     //    |         | |    \
01340     //    |   M31   | | M33  \
01341     //  n |         | |        \
01342     //    ----------------------
01343     //    1         jd         n
01344     //
01345     // We want to move M31 up one row and M33 up and to the left.
01346     //
01347     DenseLinAlgPack::delete_row_col( jd, &DenseLinAlgPack::nonconst_tri_ele( M().gms(), BLAS_Cpp::lower ) );
01348   }
01349   if( maintain_factor_ ) {
01350     //
01351     // Here we have U partitioned as:
01352     //
01353     // 1              jd      n        
01354     // ------------------------- 1
01355     // \             | |       |
01356     //   \    U11    | |  U13  |
01357     //     \         | |       |
01358     //       \  u12->| |  u23' |
01359     //         \     | |   |   |
01360     //           \   | |  \./  |
01361     //             \ |_|_______|
01362     //        u22 -> |_|_______| jd  
01363     //                 \       |
01364     //                   \ U33 |
01365     //                     \   |
01366     //                       \ | n
01367     //
01368     // To preform the update we need to first update U33 as:
01369     //   U33'*U33 + u23'*u23 ==> U33'*U33
01370     // Then we need to move U12 (one column at a time) to the
01371     // left one column and overwrite u12.
01372     // Then we will move the updated U33 (one column at a time)
01373     // up and to the left one position to overwrite u22 and
01374     // the left part of u23'.  We then decrease the dimension
01375     // of U by one and we are finished updating the factorization.
01376     //
01377     // See RAB notes from 1/21/99 and 1/26/99 for details on this update.
01378     //
01379     const size_type n = M_size_;
01380     // Resize workspace if it has not been done yet.
01381     size_type work_size = 3 * max_size_;
01382     if(work_.dim() < work_size)
01383       work_.resize(work_size);
01384     // Update the factors
01385     {
01386       DMatrixSlice U = this->U().gms();
01387       // Update U33 where it sits.
01388       if(jd < n) {
01389         size_type _size = n-jd; // Set storage for u23, c and s
01390         Range1D rng(1,_size);
01391         DVectorSlice
01392           u23 = work_(rng),
01393           c   = work_(rng+_size),
01394           s   = work_(rng+2*_size);
01395         Range1D U_rng(jd+1,n);  // Set U33 and u23
01396         DMatrixSlice U33 = U(U_rng,U_rng);
01397         u23 = U.row(jd)(U_rng);
01398         // Update U33
01399         value_type dummy;
01400         ALAP_DCHUD_CALL (
01401           U33.col_ptr(1), U33.max_rows()
01402           ,U_rng.size(), u23.start_ptr(), &dummy, 1, 0, &dummy
01403           ,&dummy, c.start_ptr(), s.start_ptr() );
01404       }
01405       // Move U13 and U33 to delete row and column jd
01406       DenseLinAlgPack::delete_row_col( jd, &DenseLinAlgPack::nonconst_tri_ele( U, BLAS_Cpp::upper ) );
01407     }
01408   }
01409   else {
01410     factor_is_updated_ = false;
01411   }
01412   // Strink the size of M and U
01413   --M_size_;
01414 }
01415 
01416 // Overridden from Serializable
01417 
01418 // ToDo: Refactor this code and create an external utility matrix
01419 // serialization class that will convert from matrix type to matrix
01420 // type.
01421 
01422 void MatrixSymPosDefCholFactor::serialize( std::ostream &out ) const
01423 {
01424   // Write key words on top line
01425   out << build_serialization_string() << std::endl;
01426   // Set the precision (very important!)
01427   out.precision(std::numeric_limits<value_type>::digits10+4);
01428   // Write the dimmension
01429   out << M_size_ << std::endl;
01430   if(M_size_) {
01431     // Write the matrix values
01432     if( maintain_original_ ) {
01433       const DMatrixSliceSym M = this->M();
01434       write_matrix( M.gms(), M.uplo(), out );
01435     }
01436     else {
01437       const DMatrixSliceTri U = this->U();
01438       write_matrix( U.gms(), U.uplo(), out );
01439     }
01440   }
01441   // ToDo: You need to write both M and U if both are computed!
01442 }
01443 
01444 void MatrixSymPosDefCholFactor::unserialize( std::istream &in )
01445 {
01446   // Get the keywords for the matrix type
01447   std::string keywords;
01448   std::getline( in, keywords, '\n' );
01449   // For now make sure the types are exactly the same!
01450   const std::string this_keywords = build_serialization_string();
01451   TEUCHOS_TEST_FOR_EXCEPTION(
01452     this_keywords != keywords, std::logic_error
01453     ,"MatrixSymPosDefCholFactor::unserialize(...): Error, the matrix type being read in from file of "
01454     "\'"<<keywords<<"\' does not equal the type expected of \'"<<this_keywords<<"\'!"
01455     );
01456   // Read in the dimension of the matrix
01457   in >> M_size_;
01458   TEUCHOS_TEST_FOR_EXCEPTION(
01459     M_size_ < 0, std::logic_error
01460     ,"MatrixSymPosDefCholFactor::unserialize(...): Error, read in a size of M_size = "<<M_size_<<" < 0!"
01461     );
01462   allocate_storage(M_size_);
01463   // Read in the matrix into storage
01464   if(maintain_original_) {
01465     M_l_r_ = M_l_c_ = 1;
01466     DMatrixSliceSym M = this->M();
01467     read_matrix( in, M.uplo(), &M.gms() );
01468   }
01469   else {
01470     U_l_r_ = U_l_c_ = 1;
01471     DMatrixSliceTri U = this->U();
01472     read_matrix( in, U.uplo(), &U.gms() );
01473   }
01474 }
01475 
01476 // Private member functions
01477 
01478 void MatrixSymPosDefCholFactor::assert_storage() const
01479 {
01480   TEUCHOS_TEST_FOR_EXCEPT( !(  MU_store_.rows()  ) );
01481 }
01482 
01483 void MatrixSymPosDefCholFactor::allocate_storage(size_type max_size) const
01484 {
01485   namespace rcp = MemMngPack;
01486   if( allocates_storage_ && MU_store_.rows() < max_size + 1 ) {
01487     // We have the right to allocate storage so lets just do it.
01488     Teuchos::RCP<DMatrix>
01489       MU_store = Teuchos::rcp(new DMatrix( max_size + 1, max_size + 1 ));
01490     typedef MemMngPack::ReleaseResource_ref_count_ptr<DMatrix> ptr_t;
01491     const_cast<MatrixSymPosDefCholFactor*>(this)->release_resource_ptr_ = Teuchos::rcp(new ptr_t(MU_store));
01492     const_cast<MatrixSymPosDefCholFactor*>(this)->MU_store_.bind( (*MU_store)() );
01493     const_cast<MatrixSymPosDefCholFactor*>(this)->max_size_ = max_size;
01494   }
01495   else {
01496     TEUCHOS_TEST_FOR_EXCEPT( !(  MU_store_.rows()  >= max_size + 1  ) );
01497   }
01498 }
01499 
01500 void MatrixSymPosDefCholFactor::assert_initialized() const
01501 {
01502   TEUCHOS_TEST_FOR_EXCEPT( !(  M_size_  ) );
01503 }
01504 
01505 void MatrixSymPosDefCholFactor::resize_and_zero_off_diagonal(size_type n, value_type scale)
01506 {
01507   using DenseLinAlgPack::nonconst_tri_ele;
01508   using DenseLinAlgPack::assign;
01509 
01510   TEUCHOS_TEST_FOR_EXCEPT( !(  n <= my_min( MU_store_.rows(), MU_store_.cols() ) - 1  ) );
01511 
01512   // Resize the views
01513   set_view(
01514     n, scale, maintain_original_, 1, 1, maintain_factor_
01515     ,U_l_r_ ? 1 : 0, U_l_r_ ? 1 : 0 );
01516   if( maintain_original_ ) {
01517     if(!is_diagonal_ && n > 1)
01518       assign( &nonconst_tri_ele( M().gms()(2,n,1,n-1), BLAS_Cpp::lower ), 0.0 );
01519   }
01520   if( U_l_r_ ) {
01521     if(!is_diagonal_ && n > 1)
01522       assign( &nonconst_tri_ele( U().gms()(1,n-1,2,n), BLAS_Cpp::upper ), 0.0 );
01523   }
01524 }
01525 
01526 void MatrixSymPosDefCholFactor::update_factorization() const
01527 {
01528   if( factor_is_updated_ ) return; // The factor should already be updated.
01529   TEUCHOS_TEST_FOR_EXCEPTION(
01530     U_l_r_ == 0, std::logic_error
01531     ,"MatrixSymPosDefCholFactor::update_factorization() : "
01532     "Error, U_l_r == 0 was set in MatrixSymPosDefCholFactor::set_view(...) "
01533     "and therefore we can not update the factorization in the provided storage space." );
01534   // Update the factorization from scratch!
01535 #ifdef PROFILE_HACK_ENABLED
01536   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::update_factorization(...) ... update" );
01537 #endif
01538   MatrixSymPosDefCholFactor
01539     *nc_this = const_cast<MatrixSymPosDefCholFactor*>(this);
01540   DMatrixSliceTriEle U = DenseLinAlgPack::nonconst_tri_ele( nc_this->U().gms(), BLAS_Cpp::upper );
01541   DenseLinAlgPack::assign( &U, DenseLinAlgPack::tri_ele( M().gms(), BLAS_Cpp::lower ) );  // Copy in the original
01542   {
01543 #ifdef PROFILE_HACK_ENABLED
01544     ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::update_factorization(...) ... potrf" );
01545 #endif
01546     DenseLinAlgLAPack::potrf( &U );
01547   }
01548   nc_this->factor_is_updated_ = true;
01549 }
01550 
01551 std::string MatrixSymPosDefCholFactor::build_serialization_string() const
01552 {
01553   std::string str = "SYMMETRIC POS_DEF";
01554   if( !maintain_original_ )
01555     str.append(" CHOL_FACTOR");
01556   if( maintain_original_ )
01557     str.append(" LOWER");
01558   else
01559     str.append(" UPPER");
01560   return str;
01561 }
01562 
01563 void MatrixSymPosDefCholFactor::write_matrix( const DMatrixSlice &Q, BLAS_Cpp::Uplo Q_uplo, std::ostream &out )
01564 {
01565   const int Q_dim = Q.rows();
01566   if( Q_uplo == BLAS_Cpp::lower ) {
01567     for( int i = 1; i <= Q_dim; ++i ) {
01568       for( int j = 1; j <= i; ++j ) {
01569         out << " " << Q(i,j);
01570       }
01571       out << std::endl;
01572     }
01573   }
01574   else {
01575     for( int i = 1; i <= Q_dim; ++i ) {
01576       for( int j = i; j <= Q_dim; ++j ) {
01577         out << " " << Q(i,j);
01578       }
01579       out << std::endl;
01580     }
01581   }
01582 }
01583 
01584 void MatrixSymPosDefCholFactor::read_matrix(  std::istream &in, BLAS_Cpp::Uplo Q_uplo, DMatrixSlice *Q_out )
01585 {
01586   DMatrixSlice &Q = *Q_out;
01587   const int Q_dim = Q.rows();
01588   if( Q_uplo == BLAS_Cpp::lower ) {
01589     for( int i = 1; i <= Q_dim; ++i ) {
01590       for( int j = 1; j <= i; ++j ) {
01591 #ifdef TEUCHOS_DEBUG
01592         TEUCHOS_TEST_FOR_EXCEPTION(
01593           in.eof(), std::logic_error
01594           ,"MatrixSymPosDefCholFactor::read_matrix(in,lower,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   else {
01602     for( int i = 1; i <= Q_dim; ++i ) {
01603       for( int j = i; j <= Q_dim; ++j ) {
01604 #ifdef TEUCHOS_DEBUG
01605         TEUCHOS_TEST_FOR_EXCEPTION(
01606           in.eof(), std::logic_error
01607           ,"MatrixSymPosDefCholFactor::read_matrix(in,upper,Q_out): Error, not finished reading in matrix yet (i="<<i<<",j="<<j<<")!"
01608           );
01609 #endif
01610         in >> Q(i,j);
01611       }
01612     }
01613   }
01614 }
01615 
01616 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends