AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef MATRIX_SYM_POS_DEF_CHOL_FACTOR_H
00030 #define MATRIX_SYM_POS_DEF_CHOL_FACTOR_H
00031 
00032 #include "AbstractLinAlgPack_Types.hpp"
00033 #include "AbstractLinAlgPack_MatrixExtractInvCholFactor.hpp"
00034 #include "AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
00035 #include "AbstractLinAlgPack_MatrixSymOpNonsingSerial.hpp"
00036 #include "AbstractLinAlgPack_MatrixSymDenseInitialize.hpp"
00037 #include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
00038 #include "AbstractLinAlgPack_MatrixSymSecant.hpp"
00039 #include "DenseLinAlgPack_DMatrixClass.hpp"
00040 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00041 #include "SerializationPack_Serializable.hpp"
00042 #include "Teuchos_RCP.hpp"
00043 #include "ReleaseResource.hpp"
00044 
00045 namespace AbstractLinAlgPack {
00106 class MatrixSymPosDefCholFactor
00107   : virtual public AbstractLinAlgPack::MatrixSymOpNonsingSerial     // doxygen needs full name
00108   , virtual public AbstractLinAlgPack::MatrixSymDenseInitialize     // ""
00109   , virtual public AbstractLinAlgPack::MatrixSymOpGetGMSSymMutable  // ""
00110   , virtual public MatrixExtractInvCholFactor
00111   , virtual public MatrixSymSecant
00112   , virtual public MatrixSymAddDelUpdateable
00113   , virtual public SerializationPack::Serializable
00114 {
00115 public:
00116   
00119 
00121   typedef Teuchos::RCP<
00122     MemMngPack::ReleaseResource>  release_resource_ptr_t;
00123 
00126   class PostMod {
00127   public:
00129     PostMod(
00130       bool   maintain_original = true
00131       ,bool  maintain_factor   = false
00132       ,bool  allow_factor      = true
00133       )
00134       :maintain_original_(maintain_original)
00135       ,maintain_factor_(maintain_factor)
00136       ,allow_factor_(allow_factor)
00137     {}
00139     void initialize(MatrixSymPosDefCholFactor* p) const
00140     {
00141       p->init_setup(
00142         NULL,Teuchos::null,0 // Allocate own storage
00143         ,maintain_original_,maintain_factor_,allow_factor_
00144         );
00145     }
00146   private:
00147     bool  maintain_original_;
00148     bool  maintain_factor_;
00149     bool  allow_factor_;
00150   }; // end PostMod
00151 
00153 
00156 
00163   MatrixSymPosDefCholFactor();
00164 
00169   MatrixSymPosDefCholFactor(
00170     DMatrixSlice                    *MU_store
00171     ,const release_resource_ptr_t&    release_resource_ptr = Teuchos::null
00172     ,size_type                        max_size             = 0
00173     ,bool                             maintain_original    = true
00174     ,bool                             maintain_factor      = false
00175     ,bool                             allow_factor         = true
00176     ,bool                             set_full_view        = true
00177     ,value_type                       scale                = 1.0
00178     );
00179 
00233   void init_setup(
00234     DMatrixSlice                    *MU_store
00235     ,const release_resource_ptr_t&    release_resource_ptr = Teuchos::null
00236     ,size_type                        max_size             = 0
00237     ,bool                             maintain_original    = true
00238     ,bool                             maintain_factor      = false
00239     ,bool                             allow_factor         = true
00240     ,bool                             set_full_view        = true
00241     ,value_type                       scale                = 1.0
00242     );
00286   void set_view(
00287     size_t            M_size
00288     ,value_type       scale
00289     ,bool             maintain_original
00290     ,size_t           M_l_r
00291     ,size_t           M_l_c
00292     ,bool             maintain_factor
00293     ,size_t           U_l_r
00294     ,size_t           U_l_c
00295     );
00298   void pivot_tols( PivotTolerances pivot_tols );
00300   PivotTolerances pivot_tols() const;
00301 
00303 
00306 
00309   void get_view_setup(
00310     size_t            *M_size
00311     ,value_type       *scale
00312     ,bool             *maintain_original
00313     ,size_t           *M_l_r
00314     ,size_t           *M_l_c
00315     ,bool             *maintain_factor
00316     ,size_t           *U_l_r
00317     ,size_t           *U_l_c
00318     ) const;
00321   bool allocates_storage() const;
00322 
00325   DMatrixSlice& MU_store();
00327   const DMatrixSlice& MU_store() const;
00330   DMatrixSliceTri U();
00332   const DMatrixSliceTri U() const;
00335   DMatrixSliceSym M();
00337   const DMatrixSliceSym M() const;
00338 
00340 
00343 
00345   size_type rows() const;
00346 
00348 
00351 
00353   void zero_out();
00355   std::ostream& output(std::ostream& out) const;
00357   bool Mp_StM(
00358     MatrixOp* m_lhs, value_type alpha
00359     ,BLAS_Cpp::Transp trans_rhs
00360     ) const;
00362   bool Mp_StM(
00363     value_type alpha,const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs
00364     );
00366   bool syrk(
00367     const MatrixOp      &mwo_rhs
00368     ,BLAS_Cpp::Transp   M_trans
00369     ,value_type         alpha
00370     ,value_type         beta
00371     );
00372 
00374 
00377 
00379   void Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00380     , const DVectorSlice& vs_rhs2, value_type beta) const;
00382   void Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00383     , const SpVectorSlice& vs_rhs2, value_type beta) const;
00385   void Vp_StPtMtV(DVectorSlice* vs_lhs, value_type alpha
00386     , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00387     , BLAS_Cpp::Transp M_rhs2_trans
00388     , const DVectorSlice& vs_rhs3, value_type beta) const;
00390   void Vp_StPtMtV(DVectorSlice* vs_lhs, value_type alpha
00391     , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00392     , BLAS_Cpp::Transp M_rhs2_trans
00393     , const SpVectorSlice& sv_rhs3, value_type beta) const;
00394 
00396 
00399 
00400   void Mp_StPtMtP( DMatrixSliceSym* sym_lhs, value_type alpha
00401     , EMatRhsPlaceHolder dummy_place_holder
00402     , const GenPermMatrixSlice& gpms_rhs, BLAS_Cpp::Transp gpms_rhs_trans
00403     , value_type beta ) const;
00404 
00406 
00409 
00411   void V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00412     , const DVectorSlice& vs_rhs2) const;
00414   void V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00415     , const SpVectorSlice& sv_rhs2) const;
00416 
00418 
00421 
00423   void M_StMtInvMtM(
00424     DMatrixSliceSym* sym_gms_lhs, value_type alpha
00425     , const MatrixOpSerial& mwo, BLAS_Cpp::Transp mwo_trans, EMatrixDummyArg
00426     ) const;
00427 
00429 
00432 
00434   void initialize( const DMatrixSliceSym& M );
00435 
00437 
00440 
00442   const DenseLinAlgPack::DMatrixSliceSym get_sym_gms_view() const;
00444   void free_sym_gms_view(const DenseLinAlgPack::DMatrixSliceSym* sym_gms_view) const;
00445 
00447 
00450 
00452   DenseLinAlgPack::DMatrixSliceSym get_sym_gms_view();
00454   void commit_sym_gms_view(DenseLinAlgPack::DMatrixSliceSym* sym_gms_view);
00455 
00457 
00460 
00462   void extract_inv_chol( DMatrixSliceTriEle* InvChol ) const;
00463   
00465 
00468 
00470   void init_identity( const VectorSpace& space_diag, value_type alpha );
00472   void init_diagonal( const Vector& diag );
00474   void secant_update(VectorMutable* s, VectorMutable* y, VectorMutable* Bs);
00475 
00477 
00480 
00482   void initialize(
00483     value_type         alpha
00484     ,size_type         max_size
00485     );
00487   void initialize(
00488     const DMatrixSliceSym      &A
00489     ,size_type         max_size
00490     ,bool              force_factorization
00491     ,Inertia           inertia
00492     ,PivotTolerances   pivot_tols
00493     );
00495   size_type max_size() const;
00497   Inertia inertia() const;
00499   void set_uninitialized();
00501   void augment_update(
00502     const DVectorSlice  *t
00503     ,value_type        alpha
00504     ,bool              force_refactorization
00505     ,EEigenValType     add_eigen_val
00506     ,PivotTolerances   pivot_tols
00507     );
00509   void delete_update(
00510     size_type          jd
00511     ,bool              force_refactorization
00512     ,EEigenValType     drop_eigen_val
00513     ,PivotTolerances   pivot_tols
00514     );
00515 
00517 
00518 
00521 
00523   void serialize( std::ostream &out ) const;
00525   void unserialize( std::istream &in );
00526 
00528 
00529 private:
00530   
00531   // /////////////////////////////
00532   // Private data members
00533 
00534   bool                    maintain_original_;
00535   bool                    maintain_factor_;
00536   bool                    factor_is_updated_;    // Is set to true if maintain_factor_=true
00537   bool                    allocates_storage_;    // If true then this object allocates the storage
00538   release_resource_ptr_t  release_resource_ptr_;
00539   DMatrixSlice            MU_store_;
00540   size_t                  max_size_;
00541   size_t                  M_size_,               // M_size == 0 is flag that we are completely uninitialized
00542                           M_l_r_,
00543                           M_l_c_,
00544                           U_l_r_,
00545                           U_l_c_;
00546   value_type              scale_;
00547   bool                    is_diagonal_;
00548   PivotTolerances         pivot_tols_;
00549   DVector                 work_; // workspace.
00550 
00551   // /////////////////////////////
00552   // Private member functions
00553 
00554   void assert_storage() const;
00555   void allocate_storage(size_type max_size) const;
00556   void assert_initialized() const;
00557   void resize_and_zero_off_diagonal(size_type n, value_type scale);
00558   void update_factorization() const;
00559   std::string build_serialization_string() const;
00560   static void write_matrix( const DMatrixSlice &Q, BLAS_Cpp::Uplo Q_uplo, std::ostream &out );
00561   static void read_matrix( std::istream &in, BLAS_Cpp::Uplo Q_uplo, DMatrixSlice *Q );
00562 
00563 }; // end class MatrixSymPosDefCholFactor
00564 
00565 // ///////////////////////////////////////////////////////
00566 // Inline members for MatrixSymPosDefCholFactor
00567 
00568 inline
00569 bool MatrixSymPosDefCholFactor::allocates_storage() const
00570 {
00571   return allocates_storage_;
00572 }
00573 
00574 inline
00575 DMatrixSlice& MatrixSymPosDefCholFactor::MU_store()
00576 {
00577   return MU_store_;
00578 }
00579 
00580 inline
00581 const DMatrixSlice& MatrixSymPosDefCholFactor::MU_store() const
00582 {
00583   return MU_store_;
00584 }
00585 
00586 inline
00587 void MatrixSymPosDefCholFactor::get_view_setup(
00588   size_t            *M_size
00589   ,value_type       *scale
00590   ,bool             *maintain_original
00591   ,size_t           *M_l_r
00592   ,size_t           *M_l_c
00593   ,bool             *maintain_factor
00594   ,size_t           *U_l_r
00595   ,size_t           *U_l_c
00596   ) const
00597 {
00598   *M_size               = M_size_;
00599   *scale                = scale_;
00600   *maintain_original    = maintain_original_;
00601   *M_l_r                = maintain_original_ ? M_l_r_ : 0;
00602   *M_l_c                = maintain_original_ ? M_l_c_ : 0;
00603   *maintain_factor      = maintain_factor_;
00604   *U_l_r                = maintain_factor_ ? U_l_r_ : 0;
00605   *U_l_c                = maintain_factor_ ? U_l_c_ : 0;
00606 }
00607 
00608 inline
00609 DMatrixSliceTri MatrixSymPosDefCholFactor::U()
00610 {
00611   return DenseLinAlgPack::nonconst_tri(
00612     MU_store_(U_l_r_,U_l_r_+M_size_-1,U_l_c_+1,U_l_c_+M_size_)
00613     , BLAS_Cpp::upper, BLAS_Cpp::nonunit
00614     );
00615 }
00616 
00617 inline
00618 const DMatrixSliceTri MatrixSymPosDefCholFactor::U() const
00619 {
00620   return DenseLinAlgPack::tri(
00621     MU_store_(U_l_r_,U_l_r_+M_size_-1,U_l_c_+1,U_l_c_+M_size_)
00622     , BLAS_Cpp::upper, BLAS_Cpp::nonunit
00623     );
00624 }
00625 
00626 inline
00627 DMatrixSliceSym MatrixSymPosDefCholFactor::M()
00628 {
00629   return DenseLinAlgPack::nonconst_sym(
00630     MU_store_(M_l_r_+1,M_l_r_+M_size_,M_l_c_,M_l_c_+M_size_-1)
00631     , BLAS_Cpp::lower
00632     );
00633 }
00634 
00635 inline
00636 const DMatrixSliceSym MatrixSymPosDefCholFactor::M() const
00637 {
00638   return DenseLinAlgPack::sym(
00639     MU_store_(M_l_r_+1,M_l_r_+M_size_,M_l_c_,M_l_c_+M_size_-1)
00640     , BLAS_Cpp::lower
00641     );
00642 }
00643 
00644 } // end namespace AbstractLinAlgPack
00645 
00646 #endif // MATRIX_SYM_POS_DEF_CHOL_FACTOR_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends