AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp

// @HEADER
// ***********************************************************************
// 
// Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
//                  Copyright (2003) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#ifndef MATRIX_SYM_POS_DEF_CHOL_FACTOR_H
#define MATRIX_SYM_POS_DEF_CHOL_FACTOR_H

#include "AbstractLinAlgPack_Types.hpp"
#include "AbstractLinAlgPack_MatrixExtractInvCholFactor.hpp"
#include "AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
#include "AbstractLinAlgPack_MatrixSymOpNonsingSerial.hpp"
#include "AbstractLinAlgPack_MatrixSymDenseInitialize.hpp"
#include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
#include "AbstractLinAlgPack_MatrixSymSecant.hpp"
#include "DenseLinAlgPack_DMatrixClass.hpp"
#include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
#include "SerializationPack_Serializable.hpp"
#include "Teuchos_RCP.hpp"
#include "ReleaseResource.hpp"

namespace AbstractLinAlgPack {
class MatrixSymPosDefCholFactor
  : virtual public AbstractLinAlgPack::MatrixSymOpNonsingSerial     // doxygen needs full name
  , virtual public AbstractLinAlgPack::MatrixSymDenseInitialize     // ""
  , virtual public AbstractLinAlgPack::MatrixSymOpGetGMSSymMutable  // ""
  , virtual public MatrixExtractInvCholFactor
  , virtual public MatrixSymSecant
  , virtual public MatrixSymAddDelUpdateable
  , virtual public SerializationPack::Serializable
{
public:
  

  typedef Teuchos::RCP<
    MemMngPack::ReleaseResource>  release_resource_ptr_t;

  class PostMod {
  public:
    PostMod(
      bool   maintain_original = true
      ,bool  maintain_factor   = false
      ,bool  allow_factor      = true
      )
      :maintain_original_(maintain_original)
      ,maintain_factor_(maintain_factor)
      ,allow_factor_(allow_factor)
    {}
    void initialize(MatrixSymPosDefCholFactor* p) const
    {
      p->init_setup(
        NULL,Teuchos::null,0 // Allocate own storage
        ,maintain_original_,maintain_factor_,allow_factor_
        );
    }
  private:
    bool  maintain_original_;
    bool  maintain_factor_;
    bool  allow_factor_;
  }; // end PostMod



  MatrixSymPosDefCholFactor();

  MatrixSymPosDefCholFactor(
    DMatrixSlice                    *MU_store
    ,const release_resource_ptr_t&    release_resource_ptr = Teuchos::null
    ,size_type                        max_size             = 0
    ,bool                             maintain_original    = true
    ,bool                             maintain_factor      = false
    ,bool                             allow_factor         = true
    ,bool                             set_full_view        = true
    ,value_type                       scale                = 1.0
    );

  void init_setup(
    DMatrixSlice                    *MU_store
    ,const release_resource_ptr_t&    release_resource_ptr = Teuchos::null
    ,size_type                        max_size             = 0
    ,bool                             maintain_original    = true
    ,bool                             maintain_factor      = false
    ,bool                             allow_factor         = true
    ,bool                             set_full_view        = true
    ,value_type                       scale                = 1.0
    );
  void set_view(
    size_t            M_size
    ,value_type       scale
    ,bool             maintain_original
    ,size_t           M_l_r
    ,size_t           M_l_c
    ,bool             maintain_factor
    ,size_t           U_l_r
    ,size_t           U_l_c
    );
  void pivot_tols( PivotTolerances pivot_tols );
  PivotTolerances pivot_tols() const;



  void get_view_setup(
    size_t            *M_size
    ,value_type       *scale
    ,bool             *maintain_original
    ,size_t           *M_l_r
    ,size_t           *M_l_c
    ,bool             *maintain_factor
    ,size_t           *U_l_r
    ,size_t           *U_l_c
    ) const;
  bool allocates_storage() const;

  DMatrixSlice& MU_store();
  const DMatrixSlice& MU_store() const;
  DMatrixSliceTri U();
  const DMatrixSliceTri U() const;
  DMatrixSliceSym M();
  const DMatrixSliceSym M() const;



  size_type rows() const;



  void zero_out();
  std::ostream& output(std::ostream& out) const;
  bool Mp_StM(
    MatrixOp* m_lhs, value_type alpha
    ,BLAS_Cpp::Transp trans_rhs
    ) const;
  bool Mp_StM(
    value_type alpha,const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs
    );
  bool syrk(
    const MatrixOp      &mwo_rhs
    ,BLAS_Cpp::Transp   M_trans
    ,value_type         alpha
    ,value_type         beta
    );



  void Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
    , const DVectorSlice& vs_rhs2, value_type beta) const;
  void Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
    , const SpVectorSlice& vs_rhs2, value_type beta) const;
  void Vp_StPtMtV(DVectorSlice* vs_lhs, value_type alpha
    , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
    , BLAS_Cpp::Transp M_rhs2_trans
    , const DVectorSlice& vs_rhs3, value_type beta) const;
  void Vp_StPtMtV(DVectorSlice* vs_lhs, value_type alpha
    , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
    , BLAS_Cpp::Transp M_rhs2_trans
    , const SpVectorSlice& sv_rhs3, value_type beta) const;



  void Mp_StPtMtP( DMatrixSliceSym* sym_lhs, value_type alpha
    , EMatRhsPlaceHolder dummy_place_holder
    , const GenPermMatrixSlice& gpms_rhs, BLAS_Cpp::Transp gpms_rhs_trans
    , value_type beta ) const;



  void V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
    , const DVectorSlice& vs_rhs2) const;
  void V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
    , const SpVectorSlice& sv_rhs2) const;



  void M_StMtInvMtM(
    DMatrixSliceSym* sym_gms_lhs, value_type alpha
    , const MatrixOpSerial& mwo, BLAS_Cpp::Transp mwo_trans, EMatrixDummyArg
    ) const;



  void initialize( const DMatrixSliceSym& M );



  const DenseLinAlgPack::DMatrixSliceSym get_sym_gms_view() const;
  void free_sym_gms_view(const DenseLinAlgPack::DMatrixSliceSym* sym_gms_view) const;



  DenseLinAlgPack::DMatrixSliceSym get_sym_gms_view();
  void commit_sym_gms_view(DenseLinAlgPack::DMatrixSliceSym* sym_gms_view);



  void extract_inv_chol( DMatrixSliceTriEle* InvChol ) const;
  


  void init_identity( const VectorSpace& space_diag, value_type alpha );
  void init_diagonal( const Vector& diag );
  void secant_update(VectorMutable* s, VectorMutable* y, VectorMutable* Bs);



  void initialize(
    value_type         alpha
    ,size_type         max_size
    );
  void initialize(
    const DMatrixSliceSym      &A
    ,size_type         max_size
    ,bool              force_factorization
    ,Inertia           inertia
    ,PivotTolerances   pivot_tols
    );
  size_type max_size() const;
  Inertia inertia() const;
  void set_uninitialized();
  void augment_update(
    const DVectorSlice  *t
    ,value_type        alpha
    ,bool              force_refactorization
    ,EEigenValType     add_eigen_val
    ,PivotTolerances   pivot_tols
    );
  void delete_update(
    size_type          jd
    ,bool              force_refactorization
    ,EEigenValType     drop_eigen_val
    ,PivotTolerances   pivot_tols
    );




  void serialize( std::ostream &out ) const;
  void unserialize( std::istream &in );


private:
  
  // /////////////////////////////
  // Private data members

  bool                    maintain_original_;
  bool                    maintain_factor_;
  bool                    factor_is_updated_;    // Is set to true if maintain_factor_=true
  bool                    allocates_storage_;    // If true then this object allocates the storage
  release_resource_ptr_t  release_resource_ptr_;
  DMatrixSlice            MU_store_;
  size_t                  max_size_;
  size_t                  M_size_,               // M_size == 0 is flag that we are completely uninitialized
                          M_l_r_,
                          M_l_c_,
                          U_l_r_,
                          U_l_c_;
  value_type              scale_;
  bool                    is_diagonal_;
  PivotTolerances         pivot_tols_;
  DVector                 work_; // workspace.

  // /////////////////////////////
  // Private member functions

  void assert_storage() const;
  void allocate_storage(size_type max_size) const;
  void assert_initialized() const;
  void resize_and_zero_off_diagonal(size_type n, value_type scale);
  void update_factorization() const;
  std::string build_serialization_string() const;
  static void write_matrix( const DMatrixSlice &Q, BLAS_Cpp::Uplo Q_uplo, std::ostream &out );
  static void read_matrix( std::istream &in, BLAS_Cpp::Uplo Q_uplo, DMatrixSlice *Q );

}; // end class MatrixSymPosDefCholFactor

// ///////////////////////////////////////////////////////
// Inline members for MatrixSymPosDefCholFactor

inline
bool MatrixSymPosDefCholFactor::allocates_storage() const
{
  return allocates_storage_;
}

inline
DMatrixSlice& MatrixSymPosDefCholFactor::MU_store()
{
  return MU_store_;
}

inline
const DMatrixSlice& MatrixSymPosDefCholFactor::MU_store() const
{
  return MU_store_;
}

inline
void MatrixSymPosDefCholFactor::get_view_setup(
  size_t            *M_size
  ,value_type       *scale
  ,bool             *maintain_original
  ,size_t           *M_l_r
  ,size_t           *M_l_c
  ,bool             *maintain_factor
  ,size_t           *U_l_r
  ,size_t           *U_l_c
  ) const
{
  *M_size               = M_size_;
  *scale                = scale_;
  *maintain_original    = maintain_original_;
  *M_l_r                = maintain_original_ ? M_l_r_ : 0;
  *M_l_c                = maintain_original_ ? M_l_c_ : 0;
  *maintain_factor      = maintain_factor_;
  *U_l_r                = maintain_factor_ ? U_l_r_ : 0;
  *U_l_c                = maintain_factor_ ? U_l_c_ : 0;
}

inline
DMatrixSliceTri MatrixSymPosDefCholFactor::U()
{
  return DenseLinAlgPack::nonconst_tri(
    MU_store_(U_l_r_,U_l_r_+M_size_-1,U_l_c_+1,U_l_c_+M_size_)
    , BLAS_Cpp::upper, BLAS_Cpp::nonunit
    );
}

inline
const DMatrixSliceTri MatrixSymPosDefCholFactor::U() const
{
  return DenseLinAlgPack::tri(
    MU_store_(U_l_r_,U_l_r_+M_size_-1,U_l_c_+1,U_l_c_+M_size_)
    , BLAS_Cpp::upper, BLAS_Cpp::nonunit
    );
}

inline
DMatrixSliceSym MatrixSymPosDefCholFactor::M()
{
  return DenseLinAlgPack::nonconst_sym(
    MU_store_(M_l_r_+1,M_l_r_+M_size_,M_l_c_,M_l_c_+M_size_-1)
    , BLAS_Cpp::lower
    );
}

inline
const DMatrixSliceSym MatrixSymPosDefCholFactor::M() const
{
  return DenseLinAlgPack::sym(
    MU_store_(M_l_r_+1,M_l_r_+M_size_,M_l_c_,M_l_c_+M_size_-1)
    , BLAS_Cpp::lower
    );
}

} // end namespace AbstractLinAlgPack

#endif // MATRIX_SYM_POS_DEF_CHOL_FACTOR_H