RTOpPack_Types.hpp

// @HEADER
// *****************************************************************************
// RTOp: Interfaces and Support Software for Vector Reduction Transformation
//       Operations
//
// Copyright 2006 NTESS and the RTOp contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER


#ifndef RTOPPACK_TYPES_HPP
#define RTOPPACK_TYPES_HPP


#include "RTOp_ConfigDefs.hpp"
#include "Teuchos_Array.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_ArrayRCP.hpp"
#include "Teuchos_Range1D.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include "Teuchos_TypeNameTraits.hpp"
#include "Teuchos_Assert.hpp"
#include "Teuchos_implicit_cast.hpp"
#include "Teuchos_FancyOStream.hpp"


namespace RTOpPack {


//
// Basic types
//

typedef Teuchos_Ordinal Ordinal;
using Teuchos::Ptr;
using Teuchos::RCP;
using Teuchos::ArrayRCP;
using Teuchos::ArrayView;
using Teuchos::Array;
using Teuchos::Range1D;
using Teuchos::ScalarTraits;
using Teuchos::TypeNameTraits;
using Teuchos::FancyOStream;

typedef Teuchos_Ordinal index_type;
typedef char  char_type;


//
// Exceptions
//


class UnknownError : public std::logic_error
{public: UnknownError(const std::string& what_arg) : std::logic_error(what_arg) {}};
class InvalidUsage : public std::logic_error
{public: InvalidUsage(const std::string& what_arg) : std::logic_error(what_arg) {}};
class InvalidNumVecs : public std::logic_error
{public: InvalidNumVecs(const std::string& what_arg) : std::logic_error(what_arg) {}};
class InvalidNumTargVecs : public std::logic_error
{public: InvalidNumTargVecs(const std::string& what_arg) : std::logic_error(what_arg) {}};
class IncompatibleVecs : public std::logic_error
{public: IncompatibleVecs(const std::string& what_arg) : std::logic_error(what_arg) {}};
class IncompatibleReductObj : public std::logic_error
{public: IncompatibleReductObj(const std::string& what_arg) : std::logic_error(what_arg) {}};


//
// VectorBase subviews
//


template<class Scalar>
class ConstSubVectorView {
public:
  ConstSubVectorView() : globalOffset_(0), subDim_(0), stride_(0) {}
  ConstSubVectorView(const ArrayRCP<const Scalar> &values_in)
    :globalOffset_(0), subDim_(0), stride_(0)
    { initialize(0, values_in.size(), values_in, 1); }
  ConstSubVectorView(Ordinal globalOffset_in, Ordinal subDim_in,
    const ArrayRCP<const Scalar> &values_in, ptrdiff_t stride_in)
    :globalOffset_(0), subDim_(0), stride_(0)
    { initialize(globalOffset_in, subDim_in, values_in, stride_in); }
  ConstSubVectorView( const ConstSubVectorView<Scalar>& sv )
    :globalOffset_(sv.globalOffset()), subDim_(sv.subDim()),
     values_(sv.values()), stride_(sv.stride()) 
    {}
  void initialize(Ordinal globalOffset_in, Ordinal subDim_in,
    const ArrayRCP<const Scalar> &values_in, ptrdiff_t stride_in)
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_ASSERT(globalOffset_in >= 0);
      if (!is_null(values_in)) {
        TEUCHOS_ASSERT(subDim_in >= 0);
        TEUCHOS_ASSERT(stride_in != 0);
        TEUCHOS_ASSERT(
          subDim_in*std::abs(Teuchos::as<int>(stride_in)) - 1 <= values_in.upperOffset());
        TEUCHOS_ASSERT(values_in.lowerOffset() <= 0);
      }
      else {
        TEUCHOS_ASSERT(subDim_in==0);
      }
#endif
      globalOffset_=globalOffset_in;
      subDim_=subDim_in;
      values_=values_in;
      stride_=stride_in;
    }
  void uninitialize()
    { globalOffset_ = 0; subDim_=0; values_ = Teuchos::null; stride_ = 0; }
  void setGlobalOffset(Ordinal globalOffset_in)
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_ASSERT(globalOffset_in >= 0);
#endif
      globalOffset_ = globalOffset_in;
    } 
  Ordinal globalOffset() const { return globalOffset_; }
  Ordinal subDim() const { return subDim_; }
  const ArrayRCP<const Scalar>  values() const { return values_; }
  ptrdiff_t stride() const { return stride_; }
  const Scalar& operator[](Ordinal i) const
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(i, 0, subDim_);
#endif
      return valuesBegin()[stride_*i];
    }
  const Scalar& operator()(Ordinal i) const { return (*this)[i]; }
private:
  Ordinal globalOffset_;
  Ordinal subDim_;
  ArrayRCP<const Scalar> values_;
  ptrdiff_t stride_;
  const typename ArrayRCP<const Scalar>::iterator valuesBegin() const
    {
      if (stride_ > 0)
        return values_.begin();
      return values_.begin() + (subDim_*std::abs(Teuchos::as<int>(stride_)) - 1);
    } 
public:
};


template<class Scalar>
class SubVectorView : public ConstSubVectorView<Scalar> {
public:
  SubVectorView() {}
  SubVectorView(const ArrayRCP<Scalar> &values_in)
    :ConstSubVectorView<Scalar>(values_in)
    {}
  SubVectorView(Ordinal globalOffset_in, Ordinal subDim_in,
    const ArrayRCP<Scalar> &values_in, ptrdiff_t stride_in)
    :ConstSubVectorView<Scalar>(globalOffset_in, subDim_in, values_in, stride_in)
    {}
  SubVectorView(Ordinal subDim_in)
    :ConstSubVectorView<Scalar>(0, subDim_in, Teuchos::arcp<Scalar>(subDim_in), 1)
    {}
  SubVectorView(const SubVectorView<Scalar> & sv)
    :ConstSubVectorView<Scalar>(sv)
    {}
  void initialize(Ordinal globalOffset_in, Ordinal subDim_in,
    const ArrayRCP<Scalar> &values_in, ptrdiff_t stride_in)
    { ConstSubVectorView<Scalar>::initialize(globalOffset_in, subDim_in, values_in, stride_in); }
  const ArrayRCP<Scalar> values() const
    { return Teuchos::arcp_const_cast<Scalar>(ConstSubVectorView<Scalar>::values());  }
  Scalar& operator[](Ordinal i) const
    { return const_cast<Scalar&>(ConstSubVectorView<Scalar>::operator[](i)); }
  Scalar& operator()(Ordinal i) const { return (*this)[i]; }
public:
};


template<class Scalar>
void assign_entries( const Ptr<const SubVectorView<Scalar> > &msv,
  const ConstSubVectorView<Scalar> &sv )
{
#ifdef TEUCHOS_DEBUG
  TEUCHOS_ASSERT_EQUALITY(msv->subDim(), sv.subDim());
#endif
  for( int i = 0; i < sv.subDim(); ++i ) {
    (*msv)(i) = sv(i);
  }
}


template<class Scalar>
std::ostream& operator<<(std::ostream &out, const ConstSubVectorView<Scalar> &sv)
{
  out
    << "{"
    << "globalOffset="<<sv.globalOffset()
    << ",subDim="<<sv.subDim()
    << ",values="<<sv.values()
    << ",stride="<<sv.stride()
    << "}";
  return out;
}


//
// MultiVectorBase subviews
//


template<class Scalar>
class ConstSubMultiVectorView {
public:
  ConstSubMultiVectorView()
    :globalOffset_(0), subDim_(0), colOffset_(0), numSubCols_(0),
     leadingDim_(0)
    {}
  ConstSubMultiVectorView(
    Ordinal globalOffset_in, Ordinal subDim_in,
    Ordinal colOffset_in, Ordinal numSubCols_in,
    const ArrayRCP<const Scalar> &values_in, Ordinal leadingDim_in
    )
    :globalOffset_(0), subDim_(0), colOffset_(0), numSubCols_(0),
     leadingDim_(0)
    {
      initialize(globalOffset_in, subDim_in, colOffset_in, numSubCols_in, values_in,
        leadingDim_in);
    }
  ConstSubMultiVectorView( const ConstSubMultiVectorView<Scalar>& smv )
    :globalOffset_(smv.globalOffset()), subDim_(smv.subDim()),
     colOffset_(smv.colOffset()), numSubCols_(smv.numSubCols()),
     values_(smv.values()), leadingDim_(smv.leadingDim())
    {}
  void initialize(
    Ordinal globalOffset_in, Ordinal subDim_in,
    Ordinal colOffset_in, Ordinal numSubCols_in,
    const ArrayRCP<const Scalar> &values_in, Ordinal leadingDim_in
    )
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_ASSERT(globalOffset_in >= 0);
      TEUCHOS_ASSERT(colOffset_in >= 0);
      if (!is_null(values_in)) {
        TEUCHOS_ASSERT(subDim_in >= 0);
        TEUCHOS_ASSERT(leadingDim_in >= subDim_in);
        TEUCHOS_ASSERT(numSubCols_in*leadingDim_in - 1 <= values_in.upperOffset());
        TEUCHOS_ASSERT(values_in.lowerOffset() <= 0);
      }
      else {
        TEUCHOS_ASSERT(subDim_in == 0);
      }
#endif
      globalOffset_=globalOffset_in;
      subDim_=subDim_in;
      colOffset_=colOffset_in;
      numSubCols_=numSubCols_in;
      values_=values_in;
      leadingDim_=leadingDim_in;
    }
  void uninitialize()
    {
      globalOffset_ = 0; subDim_=0; colOffset_=0, numSubCols_=0;
      values_=Teuchos::null; leadingDim_=0;
    }
  void setGlobalOffset(Ordinal globalOffset_in)
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_ASSERT(globalOffset_in >= 0);
#endif
      globalOffset_ = globalOffset_in;
    } 
  Ordinal globalOffset() const { return globalOffset_; }
  Ordinal subDim() const { return subDim_; }
  Ordinal colOffset() const { return colOffset_; }
  Ordinal numSubCols() const { return numSubCols_; }
  const ArrayRCP<const Scalar> values() const { return values_; }
  Ordinal leadingDim() const { return leadingDim_; }
  const Scalar& operator()(Ordinal i, Ordinal j) const
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(i, 0, subDim_);
      TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(j, 0, numSubCols_ );
#endif
      return values_[ i + leadingDim_*j ];
    }
  ConstSubVectorView<Scalar> col( const Ordinal j ) const
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(j, 0, numSubCols_ );
#endif
      return ConstSubVectorView<Scalar>(
        globalOffset(), subDim(), values().persistingView(j*leadingDim(),subDim()), 1 );
    }
private:
  Ordinal globalOffset_;
  Ordinal subDim_;
  Ordinal colOffset_;
  Ordinal numSubCols_;
  ArrayRCP<const Scalar> values_;
  Ordinal leadingDim_;
public:
};


template<class Scalar>
class SubMultiVectorView : public ConstSubMultiVectorView<Scalar> {
public:
  SubMultiVectorView() {}
  SubMultiVectorView(
    Ordinal numRows_in, Ordinal numCols_in
    )
    :ConstSubMultiVectorView<Scalar>(0, numRows_in, 0, numCols_in,
      Teuchos::arcp<Scalar>(numRows_in*numCols_in), numRows_in)
    {}
  SubMultiVectorView(
    Ordinal globalOffset_in, Ordinal subDim_in,
    Ordinal colOffset_in, Ordinal numSubCols_in,
    const ArrayRCP<Scalar> &values_in, Ordinal leadingDim_in
    )
    :ConstSubMultiVectorView<Scalar>(globalOffset_in, subDim_in,
      colOffset_in, numSubCols_in, values_in, leadingDim_in)
    {}
  SubMultiVectorView( const SubMultiVectorView<Scalar> & smv)
    :ConstSubMultiVectorView<Scalar>(smv)
    {}
 void initialize(
   Ordinal globalOffset_in, Ordinal subDim_in,
   Ordinal colOffset_in, Ordinal numSubCols_in,
   const ArrayRCP<Scalar> &values_in, Ordinal leadingDim_in
   )
   {
     ConstSubMultiVectorView<Scalar>::initialize(globalOffset_in,
       subDim_in, colOffset_in, numSubCols_in, values_in, leadingDim_in);
   }
  const ArrayRCP<Scalar> values() const
    {
      return Teuchos::arcp_const_cast<Scalar>(
        ConstSubMultiVectorView<Scalar>::values());
    }
  Scalar& operator()(Ordinal i, Ordinal j) const
    { return const_cast<Scalar&>(ConstSubMultiVectorView<Scalar>::operator()(i,j)); }
  SubVectorView<Scalar> col( const Ordinal j ) const
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(j, 0, this->numSubCols());
#endif
      return SubVectorView<Scalar>(this->globalOffset(), this->subDim(),
        values().persistingView(j*this->leadingDim(),this->subDim()), 1);
    }
public:
};


template<class Scalar>
void assign_entries( const Ptr<const SubMultiVectorView<Scalar> > &msmv,
  const ConstSubMultiVectorView<Scalar> &smv )
{
#ifdef TEUCHOS_DEBUG
  TEUCHOS_ASSERT_EQUALITY(msmv->subDim(), smv.subDim());
  TEUCHOS_ASSERT_EQUALITY(msmv->numSubCols(), smv.numSubCols());
#endif
  for( Ordinal j = 0; j < smv.numSubCols(); ++j ) {
    for( Ordinal i = 0; i < smv.subDim(); ++i ) {
      (*msmv)(i,j) = smv(i,j);
    }
  }
}


//
// Primitive Type Traits
//


template <class Scalar, class ConcreteObj>
class PrimitiveTypeTraits {
public:
  typedef Scalar primitiveType;
  static int numPrimitiveObjs()
    { return Scalar::this_type_is_missing_a_specialization(); }
  static int numIndexObjs()
    { return Scalar::this_type_is_missing_a_specialization(); }
  static int numCharObjs()
    { return Scalar::this_type_is_missing_a_specialization(); }
  static void extractPrimitiveObjs(
    const Scalar &obj,
    const ArrayView<primitiveType> &primitiveObjs,
    const ArrayView<index_type> &indexObjs,
    const ArrayView<char> &charObjs
    )
    {
      Scalar::this_type_is_missing_a_specialization(obj);
    }
  static void loadPrimitiveObjs(
    const ArrayView<const primitiveType> &primitiveObjs,
    const ArrayView<const index_type> &indexObjs,
    const ArrayView<const char> &charObjs,
    const Ptr<Scalar> &obj
    )
    {
      *obj = Scalar::this_type_is_missing_a_specialization();
    }
};



template <class Scalar>
class PrimitiveTypeTraits<Scalar, Scalar> {
public:
  typedef Scalar primitiveType;
  static int numPrimitiveObjs() { return 1; }
  static int numIndexObjs() { return 0; }
  static int numCharObjs() { return 0; }
  static void extractPrimitiveObjs(
    const Scalar &obj,
    const ArrayView<primitiveType> &primitiveObjs,
    const ArrayView<index_type> &indexObjs,
    const ArrayView<char> &charObjs
    )
    {
      assertInput(primitiveObjs, indexObjs, charObjs);
      primitiveObjs[0] = obj;
    }
  static void loadPrimitiveObjs(
    const ArrayView<const primitiveType> &primitiveObjs,
    const ArrayView<const index_type> &indexObjs,
    const ArrayView<const char> &charObjs,
    const Ptr<Scalar> &obj
    )
    {
      assertInput(primitiveObjs, indexObjs, charObjs);
      *obj = primitiveObjs[0];
    }
private:
  static void assertInput(
    const ArrayView<const primitiveType> &primitiveObjs,
    const ArrayView<const index_type> &indexObjs,
    const ArrayView<const char> &charObjs
    )
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT( primitiveObjs.size()!=1 || indexObjs.size()!=0
        || charObjs.size()!=0 );
#else
      (void)primitiveObjs;
      (void)indexObjs;
      (void)charObjs;
#endif
    }
};


template <class Scalar>
class PrimitiveTypeTraits<Scalar, index_type> {
public:
  typedef PrimitiveTypeTraits<Scalar,Scalar> ScalarPrimitiveTypeTraits;
  typedef typename ScalarPrimitiveTypeTraits::primitiveType primitiveType;
  static int numPrimitiveObjs() { return 0; }
  static int numIndexObjs() { return 1; }
  static int numCharObjs() { return 0; }
  static void extractPrimitiveObjs(
    const index_type &obj,
    const ArrayView<primitiveType> &primitiveObjs,
    const ArrayView<index_type> &indexObjs,
    const ArrayView<char> &charObjs
    )
    {
      assertInput(primitiveObjs, indexObjs, charObjs);
      indexObjs[0] = obj;
    }
  static void loadPrimitiveObjs(
    const ArrayView<const primitiveType> &primitiveObjs,
    const ArrayView<const index_type> &indexObjs,
    const ArrayView<const char> &charObjs,
    const Ptr<index_type> &obj
    )
    {
      assertInput(primitiveObjs, indexObjs, charObjs);
      *obj = indexObjs[0];
    }
private:
  static void assertInput(
    const ArrayView<const primitiveType> &primitiveObjs,
    const ArrayView<const index_type> &indexObjs,
    const ArrayView<const char> &charObjs
    )
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT( primitiveObjs.size()!=0 || indexObjs.size()!=1
        || charObjs.size()!=0 );
#endif
    }
};


#if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)


template <class Scalar>
class PrimitiveTypeTraits<std::complex<Scalar>, std::complex<Scalar> > {
public:
  typedef PrimitiveTypeTraits<Scalar,Scalar> ScalarPrimitiveTypeTraits;
  typedef typename ScalarPrimitiveTypeTraits::primitiveType primitiveType;
  static int numPrimitiveObjs()
    { return 2*ScalarPrimitiveTypeTraits::numPrimitiveObjs(); }
  static int numIndexObjs() { return 0; }
  static int numCharObjs() { return 0; }
  static void extractPrimitiveObjs(
    const std::complex<Scalar> &obj,
    const ArrayView<primitiveType> &primitiveObjs,
    const ArrayView<index_type> &indexObjs,
    const ArrayView<char> &charObjs
    )
    {
      using Teuchos::null;
      const int numScalarPrimitiveObjs =
        ScalarPrimitiveTypeTraits::numPrimitiveObjs();
      assertInput(primitiveObjs, indexObjs, charObjs);
      ScalarPrimitiveTypeTraits::extractPrimitiveObjs(
        obj.real(), primitiveObjs(0,numScalarPrimitiveObjs), null, null );
      ScalarPrimitiveTypeTraits::extractPrimitiveObjs(
        obj.imag(), primitiveObjs(numScalarPrimitiveObjs,numScalarPrimitiveObjs), null, null );
    }
  static void loadPrimitiveObjs(
    const ArrayView<const primitiveType> &primitiveObjs,
    const ArrayView<const index_type> &indexObjs,
    const ArrayView<const char> &charObjs,
    const Ptr<std::complex<Scalar> > &obj
    )
    {
      using Teuchos::null;
      using Teuchos::outArg;
      assertInput(primitiveObjs, indexObjs, charObjs);
      const int numScalarPrimitiveObjs =
        ScalarPrimitiveTypeTraits::numPrimitiveObjs();
      Scalar real, imag;
      ScalarPrimitiveTypeTraits::loadPrimitiveObjs(
        primitiveObjs(0,numScalarPrimitiveObjs), null, null,
        outArg(real) );
      ScalarPrimitiveTypeTraits::loadPrimitiveObjs(
        primitiveObjs(numScalarPrimitiveObjs,numScalarPrimitiveObjs), null, null,
        outArg(imag) );
      *obj = std::complex<Scalar>( real, imag );
    }
private:
  static void assertInput(
    const ArrayView<const primitiveType> &primitiveObjs,
    const ArrayView<const index_type> &indexObjs,
    const ArrayView<const char> &charObjs
    )
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT(
        primitiveObjs.size()!=2*ScalarPrimitiveTypeTraits::numPrimitiveObjs()
        || indexObjs.size()!=0
        || charObjs.size()!=0 );
#endif
    }
};


template <class Scalar>
class PrimitiveTypeTraits<std::complex<Scalar>, Scalar> {
public:
  typedef PrimitiveTypeTraits<Scalar,Scalar> ScalarPrimitiveTypeTraits;
  typedef typename ScalarPrimitiveTypeTraits::primitiveType primitiveType;
  static int numPrimitiveObjs()
    { return ScalarPrimitiveTypeTraits::numPrimitiveObjs(); }
  static int numIndexObjs() { return 0; }
  static int numCharObjs() { return 0; }
  static void extractPrimitiveObjs(
    const Scalar &obj,
    const ArrayView<primitiveType> &primitiveObjs,
    const ArrayView<index_type> &indexObjs,
    const ArrayView<char> &charObjs
    )
    {
      using Teuchos::null;
      assertInput(primitiveObjs, indexObjs, charObjs);
      ScalarPrimitiveTypeTraits::extractPrimitiveObjs(
        obj, primitiveObjs, null, null );
    }
  static void loadPrimitiveObjs(
    const ArrayView<const primitiveType> &primitiveObjs,
    const ArrayView<const index_type> &indexObjs,
    const ArrayView<const char> &charObjs,
    const Ptr<Scalar > &obj
    )
    {
      using Teuchos::null;
      assertInput(primitiveObjs, indexObjs, charObjs);
      ScalarPrimitiveTypeTraits::loadPrimitiveObjs(
        primitiveObjs, null, null, obj );
    }
private:
  static void assertInput(
    const ArrayView<const primitiveType> &primitiveObjs,
    const ArrayView<const index_type> &indexObjs,
    const ArrayView<const char> &charObjs
    )
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT(
        primitiveObjs.size()!=ScalarPrimitiveTypeTraits::numPrimitiveObjs()
        || indexObjs.size()!=0
        || charObjs.size()!=0 );
#endif
    }
};


#endif // defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)



//
// Forward declaration for templated types
//


template<class Scalar>  class RTOpT;


} // namespace RTOpPack


#endif // RTOPPACK_TYPES_HPP