MatVecFields.h

Go to the documentation of this file.

/*
 * Copyright (c) 2014-2017 The University of Utah
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to
 * deal in the Software without restriction, including without limitation the
 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
 * sell copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 * IN THE SOFTWARE.
 * ----------------------------------------------------------------------------
*/

#ifndef SpatialOps_MatVecFields_h
#define SpatialOps_MatVecFields_h
#include <cmath>
#include <stdexcept>
#include <sstream>
#include <spatialops/SpatialOpsConfigure.h>
#include <spatialops/structured/SpatialField.h>
#include <spatialops/structured/SpatialFieldStore.h>

namespace SpatialOps
{

template< typename OpT >
class VecScalarOp;  
template< typename OpT >
class MatScalarOp;  
template< typename OpT >
class VecVecOp;  
template< typename OpT >
class MatVecOp;  
template< typename OpT >
class MatMatOp;  
template< typename OpT >
class MatUnaryOp;  
template< typename FieldT >
class EigenDecomposition;  
template< typename FieldT >
class LinearSolve;  
template< typename FieldT >
class MatScalDiagOffset;  

template< typename FieldT >
class FieldVector
{
private:
  std::vector<SpatFldPtr<FieldT> > vals_; 
  typedef typename FieldT::value_type ValT; 
public:

  FieldVector( const int elements,
               const MemoryWindow window,
               const BoundaryCellInfo bc,
               const GhostData ghosts,
               const short int devIdx = CPU_INDEX )
    : vals_( )
  {
    int i = elements;
    vals_.reserve ( i );
    while( i-- ){
      vals_.push_back ( SpatialFieldStore::get_from_window<FieldT>( window, bc, ghosts, devIdx ) );
    }
  }

  FieldVector( const std::vector<SpatFldPtr<FieldT> >& vec )
    : vals_( vec )
  {
  }

  template< typename OpT >
  void operator = ( const VecScalarOp< OpT >& rhsOp );

  template< typename OpT >
  void operator = ( const VecVecOp< OpT >& rhsOp );

  template< typename OpT >
  void operator = ( const MatVecOp< OpT >& rhsOp );

  template< typename OpT >
  void operator = ( const MatUnaryOp< OpT >& rhsOp );

  void operator = ( const ValT scalar );

  FieldT&       operator() ( const int j )   
  { return *(vals_[j]); }
  const FieldT& operator() ( const int j ) const   
  { return *(vals_[j]); }
  FieldT&       at ( const int j )   
  { return *(vals_[j]); }
  const FieldT& at ( const int j ) const   
  { return *(vals_[j]); }

  std::vector<SpatFldPtr<FieldT> >& raw_vector(void) 
  { return vals_; }
  const std::vector<SpatFldPtr<FieldT> >& raw_vector(void) const 
  { return vals_; }

  int elements() const    
  { return vals_.size(); }

  void add_device(short int deviceIndex)
  {
    int i = vals_.size();
    while( i-- ){
      vals_[i]->add_device(deviceIndex);
    }
  }
  void validate_device(short int deviceIndex)
  {
    int i = vals_.size();
    while( i-- ){
      vals_[i]->validate_device(deviceIndex);
    }
  }
  inline void set_device_as_active( const short int deviceIndex )
  {
    int i = vals_.size();
    while( i-- ){
      vals_[i]->set_device_as_active(deviceIndex);
    }
  }

#ifdef __CUDACC__

  inline bool is_valid(int const deviceIndex) const
  {
    int i = vals_.size();
    bool valid = true;
    while( i-- ){
      valid &= vals_[i]->is_valid(deviceIndex);
    }
    return valid;
  }
  inline int active_device_index(void) const
  {
    /*Note: Assume the first field has the same active field as all others*/
    return vals_[0]->active_device_index();
  }

  inline void set_stream( const cudaStream_t& stream )
  {
    int i = vals_.size();
    while( i-- ){
      vals_[i]->set_stream(stream);
    }
  }

  inline void set_stream( const size_t i, const cudaStream_t& stream )
  {
    vals_[i]->set_stream(stream);
  }

  inline cudaStream_t const & get_stream( const size_t i ) const
  {
    return vals_[i]->get_stream();
  }

  inline cudaEvent_t const & get_last_event( const size_t i ) const
  {
    return vals_[i]->get_last_event();
  }
#endif /* __CUDACC__ */

};

template< typename FieldT >
class FieldMatrix
{
private:
  std::vector<SpatFldPtr<FieldT> > vals_;  
  typedef typename FieldT::value_type ValT; 
public:

  FieldMatrix( const int elements,
               const MemoryWindow window,
               const BoundaryCellInfo bc,
               const GhostData ghosts,
               const short int devIdx = CPU_INDEX )
    : vals_(),
      elements_( elements )
  {
    int ij = elements * elements;
    vals_.reserve ( ij );
    while( ij-- ){
      vals_.push_back ( SpatialFieldStore::get_from_window<FieldT>( window, bc, ghosts, devIdx ) );
    }
  }

  FieldMatrix( const std::vector<SpatFldPtr<FieldT> >& vec )
    : vals_( vec ),
      elements_( std::sqrt( vec.size() ) )
  {
#ifndef NDEBUG
    double elements = std::sqrt( vec.size() );
    assert( elements == std::floor( elements ) ); //Ensure that the vector coming in could be square.
#endif
  }

  void add_to_diagonal(const typename FieldT::value_type c);

  MatUnaryOp< EigenDecomposition<FieldT> > real_eigenvalues();
  void eigenvalues(FieldVector<FieldT>& realVector,
                   FieldVector<FieldT>& complexVector);

  MatVecOp< LinearSolve< FieldT > > solve( const FieldVector< FieldT >& vec ); 

  template< typename OpT >
  void operator = ( const MatScalarOp< OpT >& rhsOp );

  template< typename OpT >
  void operator = ( const MatMatOp< OpT >& rhsOp );

  void operator = ( const ValT scalar );

  FieldT& operator() ( const int i, const int j )   
  { return *vals_[i * elements_ + j]; }

  FieldT& operator() ( const int ij )   
  { return *(vals_[ij]); }

  FieldT& at ( const int ij )   
  { return *(vals_[ij]); }

  const FieldT& operator() ( const int i, const int j ) const   
  { return *vals_[i * elements_ + j]; }

  const FieldT& operator() ( const int ij ) const   
  { return *(vals_[ij]); }

  const FieldT& at ( const int ij ) const   
  { return *(vals_[ij]); }

  FieldVector<FieldT> row (const int i) const   
  {
    std::vector<SpatFldPtr<FieldT> > rowVector(vals_.begin() + (i * elements_), vals_.begin() + ((i * elements_) + elements_));

    return FieldVector<FieldT>(rowVector);
  }

  std::vector<SpatFldPtr<FieldT> >& raw_vector(void) 
  { return vals_; }
  const std::vector<SpatFldPtr<FieldT> >& raw_vector(void) const 
  { return vals_; }

  int elements() const   
  { return elements_; }

  void add_device(short int deviceIndex)
  {
    int i = vals_.size();
    while( i-- ){
      vals_[i]->add_device(deviceIndex);
    }
  }
  void validate_device(short int deviceIndex)
  {
    int i = vals_.size();
    while( i-- ){
      vals_[i]->validate_device(deviceIndex);
    }
  }
  inline void set_device_as_active( const short int deviceIndex )
  {
    int i = vals_.size();
    while( i-- ){
      vals_[i]->set_device_as_active(deviceIndex);
    }
  }

#ifdef __CUDACC__

  inline bool is_valid(int const deviceIndex) const
  {
    int i = vals_.size();
    bool valid = true;
    while( i-- ){
      valid &= vals_[i]->is_valid(deviceIndex);
    }
    return valid;
  }
  inline int active_device_index(void) const
  {
    /*Note: Assume the first field has the same active field as all others*/
    return vals_[0]->active_device_index();
  }

  inline void set_stream( const cudaStream_t& stream )
  {
    int i = vals_.size();
    while( i-- ){
      vals_[i]->set_stream(stream);
    }
  }

  inline void set_stream( const size_t i, const cudaStream_t& stream )
  {
    vals_[i]->set_stream(stream);
  }

  inline cudaStream_t const & get_stream( const size_t i ) const
  {
    return vals_[i]->get_stream();
  }

  inline cudaEvent_t const & get_last_event( const size_t i ) const
  {
    return vals_[i]->get_last_event();
  }
#endif /* __CUDACC__ */

private:
  const int elements_;
};

template< typename FieldT >
void FieldVector<FieldT>::operator = ( const ValT scalar )
{
  int i = vals_.size();
  while( i-- ){
    *(vals_[i]) <<= scalar;
  }
}

template< typename FieldT >
void FieldMatrix<FieldT>::operator = ( const ValT scalar )
{
  int i = vals_.size();
  while( i-- ){
    *(vals_[i]) <<= scalar;
  }
}

template< typename FieldT >
template< typename OpT >
void FieldVector< FieldT >::operator = ( const VecScalarOp< OpT >& rhsOp ){
  rhsOp.eval( this );
}

template< typename FieldT >
template< typename OpT >
void FieldMatrix< FieldT >::operator = ( const MatScalarOp< OpT >& rhsOp ){
  rhsOp.eval( this );
}

template< typename FieldT >
template< typename OpT >
void FieldVector< FieldT >::operator = ( const VecVecOp< OpT >& rhsOp ){
  rhsOp.eval( this );
}

template< typename FieldT >
template< typename OpT >
void FieldVector< FieldT >::operator = ( const MatVecOp< OpT >& rhsOp ){
  rhsOp.eval( this );
}

template< typename FieldT >
template< typename OpT >
void FieldMatrix< FieldT >::operator = ( const MatMatOp< OpT >& rhsOp ){
  rhsOp.eval( this );
}

template< typename FieldT >
template< typename OpT >
void FieldVector< FieldT>::operator = ( const MatUnaryOp< OpT >& rhsOp ){
  rhsOp.eval( this );
}

} // namespace SpatialOps

#endif // SpatialOps_MatVecField_h