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 // 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 #ifndef MATRIX_SYM_POS_DEF_CHOL_FACTOR_H
00043 #define MATRIX_SYM_POS_DEF_CHOL_FACTOR_H
00044 
00045 #include "AbstractLinAlgPack_Types.hpp"
00046 #include "AbstractLinAlgPack_MatrixExtractInvCholFactor.hpp"
00047 #include "AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
00048 #include "AbstractLinAlgPack_MatrixSymOpNonsingSerial.hpp"
00049 #include "AbstractLinAlgPack_MatrixSymDenseInitialize.hpp"
00050 #include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
00051 #include "AbstractLinAlgPack_MatrixSymSecant.hpp"
00052 #include "DenseLinAlgPack_DMatrixClass.hpp"
00053 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00054 #include "SerializationPack_Serializable.hpp"
00055 #include "Teuchos_RCP.hpp"
00056 #include "ReleaseResource.hpp"
00057 
00058 namespace AbstractLinAlgPack {
00119 class MatrixSymPosDefCholFactor
00120   : virtual public AbstractLinAlgPack::MatrixSymOpNonsingSerial     // doxygen needs full name
00121   , virtual public AbstractLinAlgPack::MatrixSymDenseInitialize     // ""
00122   , virtual public AbstractLinAlgPack::MatrixSymOpGetGMSSymMutable  // ""
00123   , virtual public MatrixExtractInvCholFactor
00124   , virtual public MatrixSymSecant
00125   , virtual public MatrixSymAddDelUpdateable
00126   , virtual public SerializationPack::Serializable
00127 {
00128 public:
00129   
00132 
00134   typedef Teuchos::RCP<
00135     MemMngPack::ReleaseResource>  release_resource_ptr_t;
00136 
00139   class PostMod {
00140   public:
00142     PostMod(
00143       bool   maintain_original = true
00144       ,bool  maintain_factor   = false
00145       ,bool  allow_factor      = true
00146       )
00147       :maintain_original_(maintain_original)
00148       ,maintain_factor_(maintain_factor)
00149       ,allow_factor_(allow_factor)
00150     {}
00152     void initialize(MatrixSymPosDefCholFactor* p) const
00153     {
00154       p->init_setup(
00155         NULL,Teuchos::null,0 // Allocate own storage
00156         ,maintain_original_,maintain_factor_,allow_factor_
00157         );
00158     }
00159   private:
00160     bool  maintain_original_;
00161     bool  maintain_factor_;
00162     bool  allow_factor_;
00163   }; // end PostMod
00164 
00166 
00169 
00176   MatrixSymPosDefCholFactor();
00177 
00182   MatrixSymPosDefCholFactor(
00183     DMatrixSlice                    *MU_store
00184     ,const release_resource_ptr_t&    release_resource_ptr = Teuchos::null
00185     ,size_type                        max_size             = 0
00186     ,bool                             maintain_original    = true
00187     ,bool                             maintain_factor      = false
00188     ,bool                             allow_factor         = true
00189     ,bool                             set_full_view        = true
00190     ,value_type                       scale                = 1.0
00191     );
00192 
00246   void init_setup(
00247     DMatrixSlice                    *MU_store
00248     ,const release_resource_ptr_t&    release_resource_ptr = Teuchos::null
00249     ,size_type                        max_size             = 0
00250     ,bool                             maintain_original    = true
00251     ,bool                             maintain_factor      = false
00252     ,bool                             allow_factor         = true
00253     ,bool                             set_full_view        = true
00254     ,value_type                       scale                = 1.0
00255     );
00299   void set_view(
00300     size_t            M_size
00301     ,value_type       scale
00302     ,bool             maintain_original
00303     ,size_t           M_l_r
00304     ,size_t           M_l_c
00305     ,bool             maintain_factor
00306     ,size_t           U_l_r
00307     ,size_t           U_l_c
00308     );
00311   void pivot_tols( PivotTolerances pivot_tols );
00313   PivotTolerances pivot_tols() const;
00314 
00316 
00319 
00322   void get_view_setup(
00323     size_t            *M_size
00324     ,value_type       *scale
00325     ,bool             *maintain_original
00326     ,size_t           *M_l_r
00327     ,size_t           *M_l_c
00328     ,bool             *maintain_factor
00329     ,size_t           *U_l_r
00330     ,size_t           *U_l_c
00331     ) const;
00334   bool allocates_storage() const;
00335 
00338   DMatrixSlice& MU_store();
00340   const DMatrixSlice& MU_store() const;
00343   DMatrixSliceTri U();
00345   const DMatrixSliceTri U() const;
00348   DMatrixSliceSym M();
00350   const DMatrixSliceSym M() const;
00351 
00353 
00356 
00358   size_type rows() const;
00359 
00361 
00364 
00366   void zero_out();
00368   std::ostream& output(std::ostream& out) const;
00370   bool Mp_StM(
00371     MatrixOp* m_lhs, value_type alpha
00372     ,BLAS_Cpp::Transp trans_rhs
00373     ) const;
00375   bool Mp_StM(
00376     value_type alpha,const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs
00377     );
00379   bool syrk(
00380     const MatrixOp      &mwo_rhs
00381     ,BLAS_Cpp::Transp   M_trans
00382     ,value_type         alpha
00383     ,value_type         beta
00384     );
00385 
00387 
00390 
00392   void Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00393     , const DVectorSlice& vs_rhs2, value_type beta) const;
00395   void Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00396     , const SpVectorSlice& vs_rhs2, value_type beta) const;
00398   void Vp_StPtMtV(DVectorSlice* vs_lhs, value_type alpha
00399     , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00400     , BLAS_Cpp::Transp M_rhs2_trans
00401     , const DVectorSlice& vs_rhs3, value_type beta) const;
00403   void Vp_StPtMtV(DVectorSlice* vs_lhs, value_type alpha
00404     , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00405     , BLAS_Cpp::Transp M_rhs2_trans
00406     , const SpVectorSlice& sv_rhs3, value_type beta) const;
00407 
00409 
00412 
00413   void Mp_StPtMtP( DMatrixSliceSym* sym_lhs, value_type alpha
00414     , EMatRhsPlaceHolder dummy_place_holder
00415     , const GenPermMatrixSlice& gpms_rhs, BLAS_Cpp::Transp gpms_rhs_trans
00416     , value_type beta ) const;
00417 
00419 
00422 
00424   void V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00425     , const DVectorSlice& vs_rhs2) const;
00427   void V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00428     , const SpVectorSlice& sv_rhs2) const;
00429 
00431 
00434 
00436   void M_StMtInvMtM(
00437     DMatrixSliceSym* sym_gms_lhs, value_type alpha
00438     , const MatrixOpSerial& mwo, BLAS_Cpp::Transp mwo_trans, EMatrixDummyArg
00439     ) const;
00440 
00442 
00445 
00447   void initialize( const DMatrixSliceSym& M );
00448 
00450 
00453 
00455   const DenseLinAlgPack::DMatrixSliceSym get_sym_gms_view() const;
00457   void free_sym_gms_view(const DenseLinAlgPack::DMatrixSliceSym* sym_gms_view) const;
00458 
00460 
00463 
00465   DenseLinAlgPack::DMatrixSliceSym get_sym_gms_view();
00467   void commit_sym_gms_view(DenseLinAlgPack::DMatrixSliceSym* sym_gms_view);
00468 
00470 
00473 
00475   void extract_inv_chol( DMatrixSliceTriEle* InvChol ) const;
00476   
00478 
00481 
00483   void init_identity( const VectorSpace& space_diag, value_type alpha );
00485   void init_diagonal( const Vector& diag );
00487   void secant_update(VectorMutable* s, VectorMutable* y, VectorMutable* Bs);
00488 
00490 
00493 
00495   void initialize(
00496     value_type         alpha
00497     ,size_type         max_size
00498     );
00500   void initialize(
00501     const DMatrixSliceSym      &A
00502     ,size_type         max_size
00503     ,bool              force_factorization
00504     ,Inertia           inertia
00505     ,PivotTolerances   pivot_tols
00506     );
00508   size_type max_size() const;
00510   Inertia inertia() const;
00512   void set_uninitialized();
00514   void augment_update(
00515     const DVectorSlice  *t
00516     ,value_type        alpha
00517     ,bool              force_refactorization
00518     ,EEigenValType     add_eigen_val
00519     ,PivotTolerances   pivot_tols
00520     );
00522   void delete_update(
00523     size_type          jd
00524     ,bool              force_refactorization
00525     ,EEigenValType     drop_eigen_val
00526     ,PivotTolerances   pivot_tols
00527     );
00528 
00530 
00531 
00534 
00536   void serialize( std::ostream &out ) const;
00538   void unserialize( std::istream &in );
00539 
00541 
00542 private:
00543   
00544   // /////////////////////////////
00545   // Private data members
00546 
00547   bool                    maintain_original_;
00548   bool                    maintain_factor_;
00549   bool                    factor_is_updated_;    // Is set to true if maintain_factor_=true
00550   bool                    allocates_storage_;    // If true then this object allocates the storage
00551   release_resource_ptr_t  release_resource_ptr_;
00552   DMatrixSlice            MU_store_;
00553   size_t                  max_size_;
00554   size_t                  M_size_,               // M_size == 0 is flag that we are completely uninitialized
00555                           M_l_r_,
00556                           M_l_c_,
00557                           U_l_r_,
00558                           U_l_c_;
00559   value_type              scale_;
00560   bool                    is_diagonal_;
00561   PivotTolerances         pivot_tols_;
00562   DVector                 work_; // workspace.
00563 
00564   // /////////////////////////////
00565   // Private member functions
00566 
00567   void assert_storage() const;
00568   void allocate_storage(size_type max_size) const;
00569   void assert_initialized() const;
00570   void resize_and_zero_off_diagonal(size_type n, value_type scale);
00571   void update_factorization() const;
00572   std::string build_serialization_string() const;
00573   static void write_matrix( const DMatrixSlice &Q, BLAS_Cpp::Uplo Q_uplo, std::ostream &out );
00574   static void read_matrix( std::istream &in, BLAS_Cpp::Uplo Q_uplo, DMatrixSlice *Q );
00575 
00576 }; // end class MatrixSymPosDefCholFactor
00577 
00578 // ///////////////////////////////////////////////////////
00579 // Inline members for MatrixSymPosDefCholFactor
00580 
00581 inline
00582 bool MatrixSymPosDefCholFactor::allocates_storage() const
00583 {
00584   return allocates_storage_;
00585 }
00586 
00587 inline
00588 DMatrixSlice& MatrixSymPosDefCholFactor::MU_store()
00589 {
00590   return MU_store_;
00591 }
00592 
00593 inline
00594 const DMatrixSlice& MatrixSymPosDefCholFactor::MU_store() const
00595 {
00596   return MU_store_;
00597 }
00598 
00599 inline
00600 void MatrixSymPosDefCholFactor::get_view_setup(
00601   size_t            *M_size
00602   ,value_type       *scale
00603   ,bool             *maintain_original
00604   ,size_t           *M_l_r
00605   ,size_t           *M_l_c
00606   ,bool             *maintain_factor
00607   ,size_t           *U_l_r
00608   ,size_t           *U_l_c
00609   ) const
00610 {
00611   *M_size               = M_size_;
00612   *scale                = scale_;
00613   *maintain_original    = maintain_original_;
00614   *M_l_r                = maintain_original_ ? M_l_r_ : 0;
00615   *M_l_c                = maintain_original_ ? M_l_c_ : 0;
00616   *maintain_factor      = maintain_factor_;
00617   *U_l_r                = maintain_factor_ ? U_l_r_ : 0;
00618   *U_l_c                = maintain_factor_ ? U_l_c_ : 0;
00619 }
00620 
00621 inline
00622 DMatrixSliceTri MatrixSymPosDefCholFactor::U()
00623 {
00624   return DenseLinAlgPack::nonconst_tri(
00625     MU_store_(U_l_r_,U_l_r_+M_size_-1,U_l_c_+1,U_l_c_+M_size_)
00626     , BLAS_Cpp::upper, BLAS_Cpp::nonunit
00627     );
00628 }
00629 
00630 inline
00631 const DMatrixSliceTri MatrixSymPosDefCholFactor::U() const
00632 {
00633   return DenseLinAlgPack::tri(
00634     MU_store_(U_l_r_,U_l_r_+M_size_-1,U_l_c_+1,U_l_c_+M_size_)
00635     , BLAS_Cpp::upper, BLAS_Cpp::nonunit
00636     );
00637 }
00638 
00639 inline
00640 DMatrixSliceSym MatrixSymPosDefCholFactor::M()
00641 {
00642   return DenseLinAlgPack::nonconst_sym(
00643     MU_store_(M_l_r_+1,M_l_r_+M_size_,M_l_c_,M_l_c_+M_size_-1)
00644     , BLAS_Cpp::lower
00645     );
00646 }
00647 
00648 inline
00649 const DMatrixSliceSym MatrixSymPosDefCholFactor::M() const
00650 {
00651   return DenseLinAlgPack::sym(
00652     MU_store_(M_l_r_+1,M_l_r_+M_size_,M_l_c_,M_l_c_+M_size_-1)
00653     , BLAS_Cpp::lower
00654     );
00655 }
00656 
00657 } // end namespace AbstractLinAlgPack
00658 
00659 #endif // MATRIX_SYM_POS_DEF_CHOL_FACTOR_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends