BlockVectorView.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2018-2020 TotalEnergies
 * Copyright (c) 2019-     GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_LINEARALGEBRA_UTILITIES_BLOCKVECTORVIEW_HPP_
#define GEOS_LINEARALGEBRA_UTILITIES_BLOCKVECTORVIEW_HPP_
 
#include "common/common.hpp"
#include "linearAlgebra/common/common.hpp"
 
namespace geos
{
 
template< typename VECTOR >
class BlockVectorView
{
public:
 
  using Vector = VECTOR;
 
 
  BlockVectorView & operator=( BlockVectorView const & rhs ) = delete;
 
  BlockVectorView & operator=( BlockVectorView && rhs ) noexcept = delete;
 
  virtual ~BlockVectorView() = default;
 
 
 
  void copy( BlockVectorView const & src );
 
 
 
  void scale( real64 const factor );
 
  void zero();
 
  void rand( unsigned const seed );
 
  real64 dot( BlockVectorView const & x ) const;
 
  real64 norm2() const;
 
  real64 normInf() const;
 
  void axpy( real64 const alpha,
             BlockVectorView const & x );
 
  void axpby( real64 const alpha,
              BlockVectorView const & x,
              real64 const beta );
 
 
 
  localIndex blockSize() const;
 
  globalIndex globalSize() const;
 
  localIndex localSize() const;
 
  void print( std::ostream & os = std::cout ) const;
 
  VECTOR const & block( localIndex const blockIndex ) const
  {
    GEOS_LAI_ASSERT( m_vectors[blockIndex] != nullptr );
    return *m_vectors[blockIndex];
  }
 
  VECTOR & block( localIndex const blockIndex )
  {
    GEOS_LAI_ASSERT( m_vectors[blockIndex] != nullptr );
    return *m_vectors[blockIndex];
  }
 
  VECTOR const & operator()( localIndex const blockIndex ) const
  {
    return block( blockIndex );
  }
 
  VECTOR & operator()( localIndex const blockIndex )
  {
    return block( blockIndex );
  }
 
 
protected:
 
 
  explicit BlockVectorView( localIndex const nBlocks )
    : m_vectors( nBlocks )
  { }
 
  BlockVectorView( BlockVectorView const & ) = default;
 
  BlockVectorView( BlockVectorView && ) = default;
 
 
  void resize( localIndex const size )
  {
    m_vectors.resize( size );
  }
 
  void setPointer( localIndex i, VECTOR * vec )
  {
    m_vectors( i ) = vec;
  }
 
private:
 
  array1d< VECTOR * > m_vectors;
};
 
template< typename VECTOR >
void BlockVectorView< VECTOR >::copy( BlockVectorView const & src )
{
  GEOS_LAI_ASSERT_EQ( blockSize(), src.blockSize() );
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    block( i ).copy( src.block( i ) );
  }
}
 
template< typename VECTOR >
void BlockVectorView< VECTOR >::scale( real64 const factor )
{
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    block( i ).scale( factor );
  }
}
 
template< typename VECTOR >
void BlockVectorView< VECTOR >::zero()
{
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    block( i ).zero();
  }
}
 
template< typename VECTOR >
void BlockVectorView< VECTOR >::rand( unsigned const seed )
{
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    block( i ).rand( seed );
  }
}
 
template< typename VECTOR >
real64 BlockVectorView< VECTOR >::dot( BlockVectorView const & src ) const
{
  GEOS_LAI_ASSERT_EQ( blockSize(), src.blockSize() );
  real64 accum = 0;
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    accum += block( i ).dot( src.block( i ) );
  }
  return accum;
}
 
template< typename VECTOR >
real64 BlockVectorView< VECTOR >::norm2() const
{
  real64 accum = 0;
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    real64 const temp = block( i ).norm2();
    accum = accum + temp * temp;
  }
  return std::sqrt( accum );
}
 
template< typename VECTOR >
real64 BlockVectorView< VECTOR >::normInf() const
{
  real64 result = 0.0;
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    result = std::fmax( result, block( i ).normInf() );
  }
  return result;
}
 
template< typename VECTOR >
void BlockVectorView< VECTOR >::axpy( real64 const alpha,
                                      BlockVectorView const & x )
{
  GEOS_LAI_ASSERT_EQ( blockSize(), x.blockSize() );
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    block( i ).axpy( alpha, x.block( i ) );
  }
}
 
template< typename VECTOR >
void BlockVectorView< VECTOR >::axpby( real64 const alpha,
                                       BlockVectorView const & x,
                                       real64 const beta )
{
  GEOS_LAI_ASSERT_EQ( blockSize(), x.blockSize() );
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    block( i ).axpby( alpha, x.block( i ), beta );
  }
}
 
template< typename VECTOR >
localIndex BlockVectorView< VECTOR >::blockSize() const
{
  return m_vectors.size();
}
 
template< typename VECTOR >
globalIndex BlockVectorView< VECTOR >::globalSize() const
{
  globalIndex size = 0;
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    size += block( i ).globalSize();
  }
  return size;
}
 
template< typename VECTOR >
localIndex BlockVectorView< VECTOR >::localSize() const
{
  localIndex size = 0;
  for( localIndex i = 0; i < blockSize(); i++ )
  {
    size += block( i ).localSize();
  }
  return size;
}
 
template< typename VECTOR >
void BlockVectorView< VECTOR >::print( std::ostream & os ) const
{
  for( localIndex i = 0; i < m_vectors.size(); i++ )
  {
    os << "Block " << i << " of " << blockSize() << ":" << std::endl;
    os << "=============" << std::endl;
    block( i ).print( os );
  }
}
 
template< typename VECTOR >
std::ostream & operator<<( std::ostream & os,
                           BlockVectorView< VECTOR > const & vec )
{
  vec.print( os );
  return os;
}
 
} // end geosx namespace
 
 
#endif /* GEOS_LINEARALGEBRA_UTILITIES_BLOCKVECTORVIEW_HPP_ */