RTOpPack_Types.hpp

// @HEADER
// ***********************************************************************
// 
//      Thyra: Interfaces and Support Code for the Interoperability of Abstract Numerical Algorithms 
//                 Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#ifndef RTOPPACK_TYPES_HPP
#define RTOPPACK_TYPES_HPP

#include "RTOp_config.h"
#include "Teuchos_TestForException.hpp"

namespace RTOpPack {

//
// Basic types
//

typedef RTOp_value_type value_type;
typedef RTOp_index_type index_type;
typedef RTOp_char_type  char_type;
typedef RTOp_index_type Index;

//
// 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) {}};

//
// VectorBase subviews
//

template<class Scalar>
class SubVectorT {
public:
  SubVectorT() : globalOffset_(0), subDim_(0), values_(NULL), stride_(0) {}
  SubVectorT(RTOp_index_type globalOffset, RTOp_index_type subDim, const Scalar *values, ptrdiff_t stride)
    :globalOffset_(globalOffset), subDim_(subDim), values_(values), stride_(stride) 
    {}
  SubVectorT( const SubVectorT<Scalar>& sv )
    :globalOffset_(sv.globalOffset()), subDim_(sv.subDim()), values_(sv.values()), stride_(sv.stride()) 
    {}
  void initialize(RTOp_index_type globalOffset, RTOp_index_type subDim, const Scalar *values, ptrdiff_t stride)
    { globalOffset_=globalOffset; subDim_=subDim; values_=values; stride_=stride;  }
  void set_uninitialized()
    { globalOffset_ = 0; subDim_=0; values_=NULL; stride_ = 0; }
  void setGlobalOffset(RTOp_index_type globalOffset) { globalOffset_ = globalOffset; } 
  RTOp_index_type   globalOffset() const { return globalOffset_; }
  RTOp_index_type   subDim()       const { return subDim_;  }
  const Scalar*     values()       const { return values_;  }
  ptrdiff_t         stride()       const { return stride_;  }
  const Scalar& operator[](RTOp_index_type i) const { return values_[ stride_*i ]; }
  const Scalar& operator()(RTOp_index_type i) const { return values_[ stride_*(i-1) ]; }
private:
  RTOp_index_type     globalOffset_;
  RTOp_index_type     subDim_;
  const Scalar        *values_;
  ptrdiff_t           stride_;
};

template<class Scalar>
class MutableSubVectorT : public SubVectorT<Scalar> {
public:
  MutableSubVectorT() {}
  MutableSubVectorT(RTOp_index_type globalOffset, RTOp_index_type subDim, Scalar *values, ptrdiff_t stride)
    :SubVectorT<Scalar>(globalOffset, subDim, values, stride)
    {}
  MutableSubVectorT( const MutableSubVectorT<Scalar> & sv)
    :SubVectorT<Scalar>(sv)
    {}
  void initialize(RTOp_index_type globalOffset, RTOp_index_type subDim, Scalar *values, ptrdiff_t stride)
    { SubVectorT<Scalar>::initialize(globalOffset, subDim, values, stride); }
  void set_uninitialized()
    { SubVectorT<Scalar>::set_uninitialized(); }
  Scalar* values() const { return const_cast<Scalar*>(SubVectorT<Scalar>::values());  }
  Scalar& operator[](RTOp_index_type i) const { return const_cast<Scalar&>(SubVectorT<Scalar>::operator[](i)); }
  Scalar& operator()(RTOp_index_type i) const { return const_cast<Scalar&>(SubVectorT<Scalar>::operator()(i)); }
};

template<class Scalar>
class SparseSubVectorT {
public:
  SparseSubVectorT()
    :globalOffset_(0),subDim_(0),subNz_(0)
    ,values_(NULL),valuesStride_(0),indices_(NULL)
    ,indicesStride_(0),localOffset_(0),isSorted_(0)
    {}
  SparseSubVectorT(
    RTOp_index_type globalOffset, RTOp_index_type subDim, RTOp_index_type subNz
    ,const Scalar values[], ptrdiff_t valuesStride
    ,const RTOp_index_type indices[], ptrdiff_t indicesStride
    ,ptrdiff_t localOffset, int isSorted
    )
    :globalOffset_(globalOffset),subDim_(subDim),subNz_(subNz)
    ,values_(values),valuesStride_(valuesStride),indices_(indices)
    ,indicesStride_(indicesStride),localOffset_(localOffset),isSorted_(isSorted)
    {}
  SparseSubVectorT(
    RTOp_index_type globalOffset, RTOp_index_type subDim
    ,const Scalar values[], ptrdiff_t valuesStride
    )
    :globalOffset_(globalOffset),subDim_(subDim),subNz_(subDim)
    ,values_(values),valuesStride_(valuesStride),indices_(NULL)
    ,indicesStride_(0),localOffset_(0),isSorted_(0)
    {}
  SparseSubVectorT( const SubVectorT<Scalar>& sv )
    :globalOffset_(sv.globalOffset()), subDim_(sv.subDim()), subNz_(sv.subDim()), values_(sv.values())
    ,valuesStride_(sv.stride()), indices_(NULL),indicesStride_(0),localOffset_(0),isSorted_(0)
    {}
  void initialize(
    RTOp_index_type globalOffset, RTOp_index_type subDim, RTOp_index_type subNz
    ,const Scalar values[], ptrdiff_t valuesStride
    ,const RTOp_index_type indices[], ptrdiff_t indicesStride
    ,ptrdiff_t localOffset, int isSorted
    )
    {
      globalOffset_ = globalOffset; subDim_ = subDim; subNz_ = subNz;
      values_ = values; valuesStride_ = valuesStride; indices_ = indices;
      indicesStride_ = indicesStride; localOffset_ = localOffset; isSorted_ = isSorted;
    }
  void initialize(
    RTOp_index_type globalOffset, RTOp_index_type subDim
    ,const Scalar values[], ptrdiff_t valuesStride
    )
    {
      globalOffset_ = globalOffset; subDim_ = subDim; subNz_ = subDim;
      values_ = values; valuesStride_ = valuesStride; indices_ = NULL;
      indicesStride_ = 0; localOffset_ = 0; isSorted_ = 1;
    }
  void set_uninitialized()
    {
      globalOffset_ = 0; subDim_ = 0; subNz_ = 0;
      values_ = NULL; valuesStride_ = 0; indices_ = NULL;
      indicesStride_ = 0; localOffset_ = 0; isSorted_ = 1;
    }
  void setGlobalOffset(RTOp_index_type globalOffset) { globalOffset_ = globalOffset; } 
  RTOp_index_type                  globalOffset() const { return globalOffset_; } 
  RTOp_index_type                  subDim() const { return subDim_; }
  RTOp_index_type                  subNz() const { return subNz_; }
  const Scalar*                    values() const { return values_; }
  ptrdiff_t                        valuesStride() const { return valuesStride_; }
  const RTOp_index_type*           indices() const { return indices_; }
  ptrdiff_t                        indicesStride() const { return indicesStride_; }
  ptrdiff_t                        localOffset() const { return localOffset_; }
  int                              isSorted() const { return isSorted_; }
private:
  RTOp_index_type                  globalOffset_;
  RTOp_index_type                  subDim_;
  RTOp_index_type                  subNz_;
  const Scalar                     *values_;
  ptrdiff_t                        valuesStride_;
  const RTOp_index_type            *indices_;
  ptrdiff_t                        indicesStride_;
  ptrdiff_t                        localOffset_;
  int                              isSorted_;
};


template<class Scalar>
void assign_entries( const MutableSubVectorT<Scalar> *msv, const SubVectorT<Scalar> &sv )
{
#ifdef _DEBUG
  TEST_FOR_EXCEPT(msv==NULL);
  TEST_FOR_EXCEPT(msv->subDim() != sv.subDim());
#endif
  for( int i = 1; i <= sv.subDim(); ++i ) {
    (*msv)(i) = sv(i);
  }
}

//
// MultiVectorBase subviews
//

template<class Scalar>
class SubMultiVectorT {
public:
  SubMultiVectorT()
    :globalOffset_(0), subDim_(0), colOffset_(0), numSubCols_(0)
    ,values_(NULL), leadingDim_(0)
    {}
  SubMultiVectorT(
    RTOp_index_type globalOffset, RTOp_index_type subDim
    ,RTOp_index_type colOffset, RTOp_index_type numSubCols
    ,const Scalar *values, RTOp_index_type leadingDim
    )
    :globalOffset_(globalOffset), subDim_(subDim)
    ,colOffset_(colOffset), numSubCols_(numSubCols)
    ,values_(values), leadingDim_(leadingDim)
    {}
  SubMultiVectorT( const SubMultiVectorT<Scalar>& smv )
    :globalOffset_(smv.globalOffset()), subDim_(smv.subDim())
    ,colOffset_(smv.colOffset()), numSubCols_(smv.numSubCols())
    ,values_(smv.values()), leadingDim_(smv.leadingDim())
    {}
  void initialize(
    RTOp_index_type globalOffset, RTOp_index_type subDim
    ,RTOp_index_type colOffset, RTOp_index_type numSubCols
    ,const Scalar *values, RTOp_index_type leadingDim
    )
    { globalOffset_=globalOffset; subDim_=subDim; colOffset_=colOffset; numSubCols_=numSubCols;
    values_=values; leadingDim_=leadingDim; }
  void set_uninitialized()
    { globalOffset_ = 0; subDim_=0; colOffset_=0, numSubCols_=0; values_=NULL; leadingDim_=0; }
  void setGlobalOffset(RTOp_index_type globalOffset) { globalOffset_ = globalOffset; } 
  RTOp_index_type   globalOffset()   const { return globalOffset_; }
  RTOp_index_type   subDim()         const { return subDim_; }
  RTOp_index_type   colOffset()      const { return colOffset_; }
  RTOp_index_type   numSubCols()     const { return numSubCols_; }
  const Scalar*     values()         const { return values_; }
  RTOp_index_type   leadingDim()     const { return leadingDim_;  }
  const Scalar& operator()(RTOp_index_type i, RTOp_index_type j) const
    { return values_[ (i-1) + leadingDim_*(j-1) ]; }
  SubVectorT<Scalar> col( const RTOp_index_type j ) const
    { return SubVectorT<Scalar>(globalOffset(),subDim(),values()+(j-1)*leadingDim(),1); }
private:
  RTOp_index_type     globalOffset_;
  RTOp_index_type     subDim_;
  RTOp_index_type     colOffset_;
  RTOp_index_type     numSubCols_;
  const Scalar        *values_;
  RTOp_index_type     leadingDim_;
};

template<class Scalar>
class MutableSubMultiVectorT : public SubMultiVectorT<Scalar> {
public:
  MutableSubMultiVectorT() {}
  MutableSubMultiVectorT(
    RTOp_index_type globalOffset, RTOp_index_type subDim
    ,RTOp_index_type colOffset, RTOp_index_type numSubCols
    ,const Scalar *values, RTOp_index_type leadingDim
    )
    :SubMultiVectorT<Scalar>(globalOffset,subDim,colOffset,numSubCols,values,leadingDim)
    {}
  MutableSubMultiVectorT( const MutableSubMultiVectorT<Scalar> & smv)
    :SubMultiVectorT<Scalar>(smv)
    {}
  void initialize(
    RTOp_index_type globalOffset, RTOp_index_type subDim
    ,RTOp_index_type colOffset, RTOp_index_type numSubCols
    ,const Scalar *values, RTOp_index_type leadingDim
    )
    { SubMultiVectorT<Scalar>::initialize(globalOffset,subDim,colOffset,numSubCols,values,leadingDim); }
  void set_uninitialized()
    { SubMultiVectorT<Scalar>::set_uninitialized(); }
  Scalar* values() const { return const_cast<Scalar*>(SubMultiVectorT<Scalar>::values());  }
  Scalar& operator()(RTOp_index_type i, RTOp_index_type j) const
    { return const_cast<Scalar&>(SubMultiVectorT<Scalar>::operator()(i,j)); }
  MutableSubVectorT<Scalar> col( const RTOp_index_type j ) const
    { return MutableSubVectorT<Scalar>(SubMultiVectorT<Scalar>::globalOffset(),SubMultiVectorT<Scalar>::subDim(),values()+(j-1)*SubMultiVectorT<Scalar>::leadingDim(),1); }
};

template<class Scalar>
void assign_entries( const MutableSubMultiVectorT<Scalar> *msmv, const SubMultiVectorT<Scalar> &smv )
{
#ifdef _DEBUG
  TEST_FOR_EXCEPT(msmv==NULL);
  TEST_FOR_EXCEPT(msmv->subDim() != smv.subDim());
  TEST_FOR_EXCEPT(msmv->numSubCols() != smv.numSubCols());
#endif
  for( RTOp_index_type j = 1; j <= smv.numSubCols(); ++j ) {
    for( RTOp_index_type i = 1; i <= smv.subDim(); ++i ) {
      (*msmv)(i,j) = smv(i,j);
    }
  }
}

//
// Templated types
//

template<class Scalar>  class RTOpT;

//
// Typedefs
//

typedef SubVectorT<RTOp_value_type>              SubVector;
typedef MutableSubVectorT<RTOp_value_type>       MutableSubVector;
typedef SparseSubVectorT<RTOp_value_type>        SparseSubVector;
typedef SubMultiVectorT<RTOp_value_type>         SubMultiVector;
typedef MutableSubMultiVectorT<RTOp_value_type>  MutableSubMultiVector;
typedef RTOpT<RTOp_value_type>                   RTOp;

} // namespace RTOpPack

#endif // RTOPPACK_TYPES_HPP

---