AbstractLinAlgPack_MatrixComposite.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_COMPOSITE_STD_H
#define MATRIX_COMPOSITE_STD_H

#include <deque>

#include "AbstractLinAlgPack_MatrixOp.hpp"
#include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
#include "AbstractLinAlgPack_VectorSpace.hpp"
#include "Teuchos_RCP.hpp"
#include "ReleaseResource.hpp"

namespace AbstractLinAlgPack {

class MatrixComposite : public MatrixOp {
public:

  // ///////////////////////////////////
  // Public types

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

  struct SubVectorEntry {
    typedef Teuchos::RCP<const GenPermMatrixSlice> GPMS_ptr_t;
    SubVectorEntry(
      size_type r_l, size_type c_l, value_type beta
      ,const Range1D& rng_G
      ,const GPMS_ptr_t& G, const release_resource_ptr_t& G_release, BLAS_Cpp::Transp G_trans
      ,const Vector* v
      ,const release_resource_ptr_t& v_release, BLAS_Cpp::Transp v_trans
      )
      :r_l_(r_l),c_l_(c_l),beta_(beta),rng_G_(rng_G),G_(G),G_release_(G_release),G_trans_(G_trans)
      ,v_(v),v_release_(v_release),v_trans_(v_trans)
      {}
    bool operator==(const SubVectorEntry v)
      {
        return 
          r_l_==v.r_l_ && c_l_==v.c_l_ && beta_==v.beta_
          && rng_G_==v.rng_G_ && G_.get()==v.G_.get() && G_release_.get()==v.G_release_.get() && G_trans_==v.G_trans_
          && v_==v.v_ && v_release_.get()==v.v_release_.get() && v_trans_==v.v_trans_;
      }
    size_type                  r_l_,   
                               c_l_;   

    value_type                 beta_;  

    Range1D                    rng_G_; 

    GPMS_ptr_t                 G_;     

    release_resource_ptr_t     G_release_;
    BLAS_Cpp::Transp           G_trans_;

    const Vector         *v_;     

    release_resource_ptr_t     v_release_;
    BLAS_Cpp::Transp           v_trans_;
  }; // end struct SubVectorEntry

  typedef std::deque<SubVectorEntry> vector_list_t;

  struct SubMatrixEntry {
    typedef Teuchos::RCP<const GenPermMatrixSlice> GPMS_ptr_t;
    SubMatrixEntry(
      size_type r_l, size_type r_u, size_type c_l, size_type c_u, value_type alpha
      ,const Range1D& rng_P
      ,const GPMS_ptr_t& P, const release_resource_ptr_t& P_release, BLAS_Cpp::Transp P_trans
      ,const MatrixOp* A, const release_resource_ptr_t& A_release, BLAS_Cpp::Transp A_trans
      ,const Range1D& rng_Q
      ,const GPMS_ptr_t& Q, const release_resource_ptr_t& Q_release, BLAS_Cpp::Transp Q_trans
      )
      :r_l_(r_l),r_u_(r_u),c_l_(c_l),c_u_(c_u),alpha_(alpha),rng_P_(rng_P),P_(P),P_release_(P_release),P_trans_(P_trans)
      ,A_(A),A_release_(A_release),A_trans_(A_trans),rng_Q_(rng_Q),Q_(Q),Q_release_(Q_release),Q_trans_(Q_trans)
      {}
    bool operator==(const SubMatrixEntry m)
      {
        return
          r_l_==m.r_l_ && r_u_==m.r_u_ && c_l_==m.c_l_ && c_u_==m.c_u_ && alpha_==m.alpha_
          && rng_P_==m.rng_P_ && P_.get()==m.P_.get() && P_release_.get()==m.P_release_.get() && P_trans_==m.P_trans_
          && A_==m.A_ && A_release_.get()==m.A_release_.get() && A_trans_==m.A_trans_
          && rng_Q_==m.rng_Q_ && Q_.get()==m.Q_.get() && Q_release_.get()==m.Q_release_.get() && Q_trans_==m.Q_trans_;
      }
    size_type                  r_l_, r_u_, c_l_, c_u_;
    value_type                 alpha_;
    Range1D                    rng_P_;  // rng_P_.size() > 0 => P_ is ignored, rng_P_.full_range() => all rows op(A)
    GPMS_ptr_t                 P_;
    release_resource_ptr_t     P_release_;
    BLAS_Cpp::Transp           P_trans_;
    const MatrixOp         *A_;
    release_resource_ptr_t     A_release_;
    BLAS_Cpp::Transp           A_trans_;
    Range1D                    rng_Q_; // rng_Q_.size() > 0 => Q_ is ignored, rng_Q_.full_range() => all columns op(A)
    GPMS_ptr_t                 Q_;
    release_resource_ptr_t     Q_release_;
    BLAS_Cpp::Transp           Q_trans_;
  }; // end struct SubMatrixEntry

  typedef std::deque<SubMatrixEntry> matrix_list_t;


  MatrixComposite( size_type rows = 0, size_type cols = 0 );

  void reinitialize( size_type rows = 0, size_type cols = 0 );

  void add_vector(
    size_type                      row_offset
    ,size_type                     col_offset
    ,value_type                    beta
    ,const GenPermMatrixSlice      *G
    ,const release_resource_ptr_t  &G_release
    ,BLAS_Cpp::Transp              G_trans
    ,const Vector            *v
    ,const release_resource_ptr_t  &v_release
    ,BLAS_Cpp::Transp              v_trans
    );

  void add_vector(
    size_type                      row_offset
    ,size_type                     col_offset
    ,value_type                    beta
    ,const Range1D                 &rng_G
    ,const Vector            *v
    ,const release_resource_ptr_t  &v_release
    ,BLAS_Cpp::Transp              v_trans
    );

  void add_vector(
    size_type                      row_offset
    ,size_type                     col_offset
    ,value_type                    beta
    ,const Vector            *v
    ,const release_resource_ptr_t  &v_release
    ,BLAS_Cpp::Transp              v_trans
    );

  void remove_vector( vector_list_t::iterator itr );

  void add_matrix(
    size_type                      row_offset
    ,size_type                     col_offset
    ,value_type                    alpha
    ,const GenPermMatrixSlice      *P
    ,const release_resource_ptr_t  &P_release
    ,BLAS_Cpp::Transp              P_trans
    ,const MatrixOp            *A
    ,const release_resource_ptr_t  &A_release
    ,BLAS_Cpp::Transp              A_trans
    ,const GenPermMatrixSlice      *Q
    ,const release_resource_ptr_t  &Q_release
    ,BLAS_Cpp::Transp              Q_trans
    );

  void add_matrix(
    size_type                      row_offset
    ,size_type                     col_offset
    ,value_type                    alpha
    ,const Range1D                 &rng_P
    ,const MatrixOp            *A
    ,const release_resource_ptr_t  &A_release
    ,BLAS_Cpp::Transp              A_trans
    ,const Range1D                 &rng_Q
    );

  void add_matrix(
    size_type                      row_offset
    ,size_type                     col_offset
    ,value_type                    alpha
    ,const Range1D                 &rng_P
    ,const MatrixOp            *A
    ,const release_resource_ptr_t  &A_release
    ,BLAS_Cpp::Transp              A_trans
    ,const GenPermMatrixSlice      *Q
    ,const release_resource_ptr_t  &Q_release
    ,BLAS_Cpp::Transp              Q_trans
    );

  void add_matrix(
    size_type                      row_offset
    ,size_type                     col_offset
    ,value_type                    alpha
    ,const GenPermMatrixSlice      *P
    ,const release_resource_ptr_t  &P_release
    ,BLAS_Cpp::Transp              P_trans
    ,const MatrixOp            *A
    ,const release_resource_ptr_t  &A_release
    ,BLAS_Cpp::Transp              A_trans
    ,const Range1D                 &rng_Q
    );

  void add_matrix(
    size_type                      row_offset
    ,size_type                     col_offset
    ,value_type                    alpha
    ,const MatrixOp            *A
    ,const release_resource_ptr_t  &A_release
    ,BLAS_Cpp::Transp              A_trans
    );

  void add_matrix(
    size_type                      row_offset
    ,size_type                     col_offset
    ,value_type                    alpha
    ,const GenPermMatrixSlice      *P
    ,const release_resource_ptr_t  &P_release
    ,BLAS_Cpp::Transp              P_trans
    );

  void remove_matrix( matrix_list_t::iterator itr );

  virtual void finish_construction(
    const VectorSpace::space_ptr_t&  space_cols
    ,const VectorSpace::space_ptr_t& space_rows
    );



  int                             num_vectors() const;
  vector_list_t::iterator         vectors_begin();
  vector_list_t::iterator         vectors_end();
  vector_list_t::const_iterator   vectors_begin() const;
  vector_list_t::const_iterator   vectors_end() const;
  int                             num_matrices() const;
  matrix_list_t::iterator         matrices_begin();
  matrix_list_t::iterator         matrices_end();
  matrix_list_t::const_iterator   matrices_begin() const;
  matrix_list_t::const_iterator   matrices_end() const;



  size_type rows() const;
  size_type cols() const;
  size_type nz() const;



  const VectorSpace& space_rows() const;
  const VectorSpace& space_cols() const;
  mat_ptr_t sub_view(const Range1D& row_rng, const Range1D& col_rng) const;
  void Vp_StMtV(VectorMutable* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
    , const Vector& v_rhs2, value_type beta) const;
  void Vp_StMtV(VectorMutable* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
    , const SpVectorSlice& sv_rhs2, value_type beta) const;
  void Vp_StPtMtV(VectorMutable* vs_lhs, value_type alpha
    , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
    , BLAS_Cpp::Transp M_rhs2_trans
    , const Vector& v_rhs3, value_type beta) const;
  void Vp_StPtMtV(VectorMutable* 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;


private:

  // ///////////////////////////////
  // private data members

  bool                       fully_constructed_;
  size_type                  rows_, cols_;
#ifdef DOXYGEN_COMPILE
  MatrixOp               *matrices;
    Vector               *vectors;
    VectorSpace                *space_cols;
    VectorSpace                *space_rows;
#else
  matrix_list_t              matrix_list_;
  vector_list_t              vector_list_;
  VectorSpace::space_ptr_t   space_cols_;
  VectorSpace::space_ptr_t   space_rows_;
#endif

  // ///////////////////////////////
  // private member functions

  void assert_fully_constructed() const;

}; // end class MatrixComposite

} // end namespace AbstractLinAlgPack

#endif // MATRIX_COMPOSITE_STD_H