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 Total, S.A
 * Copyright (c) 2019-     GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */

#ifndef GEOSX_LINEARALGEBRA_UTILITIES_BLOCKVECTORVIEW_HPP_
#define GEOSX_LINEARALGEBRA_UTILITIES_BLOCKVECTORVIEW_HPP_

#include "linearAlgebra/common.hpp"

namespace geosx
{

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 = 1984 );

  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
  {
    GEOSX_LAI_ASSERT( m_vectors[blockIndex] != nullptr );
    return *m_vectors[blockIndex];
  }

  VECTOR & block( localIndex const blockIndex )
  {
    GEOSX_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 )
{
  GEOSX_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
{
  GEOSX_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 )
{
  GEOSX_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 )
{
  GEOSX_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 /* GEOSX_LINEARALGEBRA_UTILITIES_BLOCKVECTORVIEW_HPP_ */