AbstractLinAlgPack_MultiVectorMutableDense.cpp

// @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.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#include <assert.h>

#include "AbstractLinAlgPack_MultiVectorMutableDense.hpp"
#include "AbstractLinAlgPack_VectorMutableDense.hpp"
#include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
#include "AbstractLinAlgPack_SpVectorOp.hpp"
#include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
#include "DenseLinAlgPack_DMatrixOut.hpp"
#include "ReleaseResource_ref_count_ptr.hpp"
#include "Teuchos_Workspace.hpp"
#include "Teuchos_TestForException.hpp"

namespace AbstractLinAlgPack {

// Constructors / initializers

MultiVectorMutableDense::MultiVectorMutableDense(
  const size_type                    rows
  ,const size_type                   cols
  )
{
  this->initialize(rows,cols);
}

MultiVectorMutableDense::MultiVectorMutableDense(
  DMatrixSlice                     gms
  ,BLAS_Cpp::Transp                  gms_trans
  ,const release_resource_ptr_t&     gms_release
  )
{
  this->initialize(gms,gms_trans,gms_release);
}

void MultiVectorMutableDense::initialize(
  const size_type                    rows
  ,const size_type                   cols
  )
{
  namespace rcp = MemMngPack;
  namespace rmp = MemMngPack;
  typedef Teuchos::RCP<DMatrix> vec_ptr_t;
  vec_ptr_t gms_ptr = Teuchos::rcp(new DMatrix(rows,cols));
  this->initialize(
    (*gms_ptr)()
    ,BLAS_Cpp::no_trans
    ,Teuchos::rcp(
      new rmp::ReleaseResource_ref_count_ptr<DMatrix>(
        gms_ptr
        )
      )
    );
}

void MultiVectorMutableDense::initialize(
  DMatrixSlice                     gms
  ,BLAS_Cpp::Transp                  gms_trans
  ,const release_resource_ptr_t&     gms_release
  )
{
  gms_.bind(gms);
  gms_trans_   = gms_trans;
  gms_release_ = gms_release;
}

// Overridden from MatrixOpGetGMS

const DMatrixSlice MultiVectorMutableDense::get_gms_view() const
{
  if(gms_trans_ == BLAS_Cpp::trans) {
    TEST_FOR_EXCEPT(true); // ToDo: We need to create a copy and transpose it!
  }
  return get_gms(); // No memory to allocate!
}

void MultiVectorMutableDense::free_gms_view(const DMatrixSlice* gms_view) const
{
  if(gms_trans_ == BLAS_Cpp::trans) {
    TEST_FOR_EXCEPT(true); // ToDo: We need to free the copy that we created in get_gms_view()
  }
  else {
    // Nothing to free!
  }
}

// Overridden from MatrixOpGetGMSMutable

DMatrixSlice MultiVectorMutableDense::get_gms_view()
{
  if(gms_trans_ == BLAS_Cpp::trans) {
    TEST_FOR_EXCEPT(true); // ToDo: We need to create a copy and transpose it!
  }
  return set_gms(); // No memory to allocate!
}

void MultiVectorMutableDense::commit_gms_view(DMatrixSlice* gms_view)
{
  if(gms_trans_ == BLAS_Cpp::trans) {
    TEST_FOR_EXCEPT(true); // ToDo: We need to free the copy that we created in get_gms_view()
  }
  else {
    // Nothing to free!
  }
}

// Overridden from MultiVector

MultiVectorMutableDense::access_by_t
MultiVectorMutableDense::access_by() const
{
  return ROW_ACCESS | COL_ACCESS | DIAG_ACCESS; // row, column and diagonal access is available!
}

// Overridden from MultiVectorMutable

MultiVectorMutableDense::vec_mut_ptr_t
MultiVectorMutableDense::col(index_type j)
{
  namespace rcp = MemMngPack;
  return Teuchos::rcp(
    new VectorMutableDense(
      DenseLinAlgPack::col( set_gms(), gms_trans(), j )
      ,Teuchos::null ) );
}

MultiVectorMutableDense::vec_mut_ptr_t
MultiVectorMutableDense::row(index_type i)
{
  namespace rcp = MemMngPack;
  return Teuchos::rcp(
    new VectorMutableDense(
      DenseLinAlgPack::row( set_gms(), gms_trans(), i )
      ,Teuchos::null ) );
}

MultiVectorMutableDense::vec_mut_ptr_t
MultiVectorMutableDense::diag(int k)
{
  namespace rcp = MemMngPack;
  return Teuchos::rcp(
    new VectorMutableDense(
      gms_.diag( gms_trans() == BLAS_Cpp::no_trans ? k : -k )
      ,Teuchos::null ) );
}

MultiVectorMutableDense::multi_vec_mut_ptr_t
MultiVectorMutableDense::mv_sub_view(const Range1D& row_rng, const Range1D& col_rng)
{
  namespace rcp = MemMngPack;
  return Teuchos::rcp(
    new MultiVectorMutableDense(
      gms_(
        gms_trans() == BLAS_Cpp::no_trans   ? row_rng : col_rng
        ,gms_trans() == BLAS_Cpp::no_trans  ? col_rng : row_rng )
      ,gms_trans()
      ,Teuchos::null ) );
}

// Overridden from MatrixBase

size_type MultiVectorMutableDense::rows() const
{
  return BLAS_Cpp::rows( get_gms().rows(), get_gms().cols(), gms_trans() );
}

size_type MultiVectorMutableDense::cols() const
{
  return BLAS_Cpp::cols( get_gms().rows(), get_gms().cols(), gms_trans() );
}

// Overridden from MatrixOp

void MultiVectorMutableDense::zero_out()
{
  gms_ = 0.0;
}

void MultiVectorMutableDense::Mt_S( value_type alpha )
{
  DenseLinAlgPack::Mt_S(&gms_,alpha);
}

MatrixOp& MultiVectorMutableDense::operator=(const MatrixOp& mwo_rhs)
{
  DenseLinAlgPack::assign( &set_gms(), MatrixDenseEncap(mwo_rhs)(), gms_trans() );
  return *this;
}

std::ostream& MultiVectorMutableDense::output(std::ostream& out) const
{
  if(gms_trans() == BLAS_Cpp::no_trans)
    return out << gms_;
  return MatrixOpSerial::output(out);
}

bool MultiVectorMutableDense::Mp_StM(
  MatrixOp* mwo_lhs, value_type alpha
  ,BLAS_Cpp::Transp trans_rhs
  ) const
{
  return MultiVectorMutable::Mp_StM(mwo_lhs,alpha,trans_rhs);
}

bool MultiVectorMutableDense::Mp_StM(
  value_type alpha,const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs
  )
{
  return MultiVectorMutable::Mp_StM(alpha,M_rhs,trans_rhs);
}

bool MultiVectorMutableDense::syrk(
  BLAS_Cpp::Transp M_trans, value_type alpha
  ,value_type beta, MatrixSymOp* sym_lhs
  ) const
{
#ifdef TEUCHOS_DEBUG
  TEST_FOR_EXCEPTION(
    sym_lhs == NULL, std::invalid_argument
    ,"MultiVectorMutableDense::syrk(...) : Error!" );
#endif
  MatrixSymOpGetGMSSymMutable
    *sym_get_lhs = dynamic_cast<MatrixSymOpGetGMSSymMutable*>(sym_lhs);
  if(!sym_get_lhs)
    return false;
  MatrixDenseSymMutableEncap  sym_gms_lhs(sym_get_lhs);
  DenseLinAlgPack::syrk( M_trans, alpha, get_gms(), beta, &sym_gms_lhs() );
  return true;
}

bool MultiVectorMutableDense::Mp_StMtM(
  MatrixOp* mwo_lhs, value_type alpha
  ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
  ,BLAS_Cpp::Transp trans_rhs2
  ,value_type beta ) const
{
  if(MultiVector::Mp_StMtM(mwo_lhs,alpha,mwo_rhs1,trans_rhs1,trans_rhs2,beta))
    return true;
  return MatrixOpSerial::Mp_StMtM(mwo_lhs,alpha,mwo_rhs1,trans_rhs1,trans_rhs2,beta);
}

bool MultiVectorMutableDense::Mp_StMtM(
  MatrixOp* mwo_lhs, value_type alpha
  ,BLAS_Cpp::Transp trans_rhs1
  ,const MatrixOp& mwo_rhs2, BLAS_Cpp::Transp trans_rhs2
  ,value_type beta ) const
{
  if(MultiVector::Mp_StMtM(mwo_lhs,alpha,trans_rhs1,mwo_rhs2,trans_rhs2,beta))
    return true;
  return MatrixOpSerial::Mp_StMtM(mwo_lhs,alpha,trans_rhs1,mwo_rhs2,trans_rhs2,beta);
}

// Overridden from MatrixOpSerial

void MultiVectorMutableDense::Vp_StMtV(
  DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
  , const DVectorSlice& vs_rhs2, value_type beta) const
{
  DenseLinAlgPack::Vp_StMtV(
    vs_lhs,alpha,gms_,BLAS_Cpp::trans_trans(gms_trans(),trans_rhs1),vs_rhs2,beta);
}

void MultiVectorMutableDense::Vp_StMtV(
  DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
  , const SpVectorSlice& sv_rhs2, value_type beta) const
{
  AbstractLinAlgPack::Vp_StMtV(
    vs_lhs,alpha,gms_,BLAS_Cpp::trans_trans(gms_trans(),trans_rhs1),sv_rhs2,beta);
}

// protected

// Overridden from MultiVector

void MultiVectorMutableDense::apply_op(
  EApplyBy apply_by, const RTOpPack::RTOp& primary_op
  ,const size_t num_multi_vecs,      const MultiVector**   multi_vecs
  ,const size_t num_targ_multi_vecs, MultiVectorMutable**  targ_multi_vecs
  ,RTOpPack::ReductTarget* reduct_objs[]
  ,const index_type primary_first_ele   , const index_type primary_sub_dim   , const index_type primary_global_offset
  ,const index_type secondary_first_ele , const index_type secondary_sub_dim 
  ) const
{
  MultiVector::apply_op(
    apply_by, primary_op, num_multi_vecs, multi_vecs, num_targ_multi_vecs, targ_multi_vecs
    ,reduct_objs
    ,primary_first_ele, primary_sub_dim, primary_global_offset, secondary_first_ele, secondary_sub_dim
    ); // ToDo: Provide Specialized implementation if needed?
}

void MultiVectorMutableDense::apply_op(
  EApplyBy apply_by, const RTOpPack::RTOp& primary_op, const RTOpPack::RTOp& secondary_op
  ,const size_t num_multi_vecs,      const MultiVector**   multi_vecs
  ,const size_t num_targ_multi_vecs, MultiVectorMutable**  targ_multi_vecs
  ,RTOpPack::ReductTarget *reduct_obj
  ,const index_type primary_first_ele   , const index_type primary_sub_dim   , const index_type primary_global_offset
  ,const index_type secondary_first_ele , const index_type secondary_sub_dim 
  ) const
{
  MultiVector::apply_op(
    apply_by, primary_op, secondary_op, num_multi_vecs, multi_vecs, num_targ_multi_vecs, targ_multi_vecs
    ,reduct_obj
    ,primary_first_ele, primary_sub_dim, primary_global_offset, secondary_first_ele, secondary_sub_dim
    ); // ToDo: Provide Specialized implementation if needed?
}

} // end namespace AbstractLinAlgPack