BlockOperatorView.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_BLOCKOPERATORVIEW_HPP_
#define GEOSX_LINEARALGEBRA_UTILITIES_BLOCKOPERATORVIEW_HPP_

#include "codingUtilities/SFINAE_Macros.hpp"
#include "linearAlgebra/interfaces/LinearOperator.hpp"
#include "linearAlgebra/utilities/BlockVector.hpp"

namespace geosx
{

template< typename VECTOR, typename OPERATOR = LinearOperator< VECTOR > >
class BlockOperatorView : public LinearOperator< BlockVectorView< VECTOR > >
{

public:

  using Base = LinearOperator< BlockVectorView< VECTOR > >;

  using Vector = typename Base::Vector;


  virtual ~BlockOperatorView() override = default;

  BlockOperatorView & operator=( BlockOperatorView const & ) = delete;

  BlockOperatorView & operator=( BlockOperatorView && ) = delete;


private:

  HAS_MEMBER_FUNCTION( numLocalRows, localIndex, );
  HAS_MEMBER_FUNCTION( numLocalCols, localIndex, );
  HAS_MEMBER_FUNCTION_NO_RTYPE( applyTranspose, std::declval< Vector const & >(), std::declval< Vector & >() );

public:


  virtual void apply( Vector const & src,
                      Vector & dst ) const override
  {
    for( localIndex i = 0; i < numBlockRows(); i++ )
    {
      dst.block( i ).zero();
      VECTOR temp( dst.block( i ) );
      for( localIndex j = 0; j < numBlockCols(); j++ )
      {
        if( m_operators( i, j ) != nullptr )
        {
          m_operators( i, j )->apply( src.block( j ), temp );
          dst.block( i ).axpy( 1.0, temp );
        }
      }
    }
  }

  template< typename OP = OPERATOR >
  std::enable_if_t< HasMemberFunction_applyTranspose< OP > >
  applyTranspose( BlockVectorView< VECTOR > const & src,
                  BlockVectorView< VECTOR > & dst ) const
  {
    for( localIndex j = 0; j < numBlockCols(); j++ )
    {
      dst.block( j ).zero();
      VECTOR temp( dst.block( j ) );
      for( localIndex i = 0; i < numBlockRows(); i++ )
      {
        if( m_operators( j, i ) != nullptr )
        {
          m_operators( j, i )->applyTranspose( src.block( i ), temp );
          dst.block( j ).axpy( 1.0, temp );
        }
      }
    }
  }

  virtual globalIndex numGlobalRows() const override
  {
    return computeRowSize( []( OPERATOR const & block ) { return block.numGlobalRows(); } );
  }

  virtual globalIndex numGlobalCols() const override
  {
    return computeColSize( []( OPERATOR const & block ) { return block.numGlobalCols(); } );
  }

  template< typename OP = OPERATOR >
  std::enable_if_t< HasMemberFunction_numLocalRows< OP >, localIndex >
  numLocalRows() const
  {
    return computeRowSize( []( OPERATOR const & block ) { return block.numLocalRows(); } );
  }

  template< typename OP = OPERATOR >
  std::enable_if_t< HasMemberFunction_numLocalCols< OP >, localIndex >
  numLocalCols() const
  {
    return computeColSize( []( OPERATOR const & block ) { return block.numLocalCols(); } );
  }



  localIndex numBlockRows() const
  {
    return m_operators.size( 0 );
  }

  localIndex numBlockCols() const
  {
    return m_operators.size( 1 );
  }

  OPERATOR const & block( localIndex const blockRowIndex, localIndex const blockColIndex ) const
  {
    GEOSX_LAI_ASSERT( m_operators( blockRowIndex, blockColIndex ) != nullptr );
    return *m_operators( blockRowIndex, blockColIndex );
  }

  OPERATOR & block( localIndex const blockRowIndex, localIndex const blockColIndex )
  {
    GEOSX_LAI_ASSERT( m_operators( blockRowIndex, blockColIndex ) != nullptr );
    return *m_operators( blockRowIndex, blockColIndex );
  }

  OPERATOR const & operator()( localIndex const blockRowIndex, localIndex const blockColIndex = 0 ) const
  {
    return block( blockRowIndex, blockColIndex );
  }

  OPERATOR & operator()( localIndex const blockRowIndex, localIndex const blockColIndex = 0 )
  {
    return block( blockRowIndex, blockColIndex );
  }


protected:

  BlockOperatorView( localIndex const nRows, localIndex const nCols )
    : m_operators( nRows, nCols )
  {
    GEOSX_LAI_ASSERT_GT( nRows, 0 );
    GEOSX_LAI_ASSERT_GT( nCols, 0 );
  }

  BlockOperatorView( BlockOperatorView< VECTOR, OPERATOR > const & x ) = default;

  BlockOperatorView( BlockOperatorView< VECTOR, OPERATOR > && x ) = default;

  void setPointer( localIndex const blockRowIndex, localIndex const blockColIndex, OPERATOR * op )
  {
    GEOSX_LAI_ASSERT_GE( blockRowIndex, 0 );
    GEOSX_LAI_ASSERT_GT( numBlockRows(), blockRowIndex );
    GEOSX_LAI_ASSERT_GE( blockColIndex, 0 );
    GEOSX_LAI_ASSERT_GT( numBlockCols(), blockColIndex );
    m_operators( blockRowIndex, blockColIndex ) = op;
  }

private:

  template< typename FUNC >
  auto computeRowSize( FUNC func ) const -> decltype( func( std::declval< OPERATOR const >() ) );

  template< typename FUNC >
  auto computeColSize( FUNC func ) const -> decltype( func( std::declval< OPERATOR const >() ) );

  array2d< OPERATOR * > m_operators;
};

template< typename VECTOR, typename OPERATOR >
template< typename FUNC >
auto BlockOperatorView< VECTOR, OPERATOR >::computeRowSize( FUNC func ) const -> decltype( func( std::declval< OPERATOR const >() ) )
{
  using sizeType = decltype( func( std::declval< OPERATOR const >() ) );
  sizeType rowSize = 0;
  for( localIndex i = 0; i < numBlockRows(); ++i )
  {
    for( localIndex j = 0; j < numBlockCols(); ++j )
    {
      if( m_operators( i, j ) != nullptr )
      {
        rowSize += func( block( i, j ) );
        break;
      }
    }
  }
  return rowSize;
}

template< typename VECTOR, typename OPERATOR >
template< typename FUNC >
auto BlockOperatorView< VECTOR, OPERATOR >::computeColSize( FUNC func ) const -> decltype( func( std::declval< OPERATOR const >() ) )
{
  using sizeType = decltype( func( std::declval< OPERATOR const >() ) );
  sizeType colSize = 0;
  for( localIndex j = 0; j < numBlockCols(); j++ )
  {
    for( localIndex i = 0; i < numBlockRows(); ++i )
    {
      if( m_operators( i, j ) != nullptr )
      {
        colSize += func( block( i, j ) );
        break;
      }
    }
  }
  return colSize;
}

}// end geosx namespace


#endif /*GEOSX_LINEARALGEBRA_UTILITIES_BLOCKOPERATORVIEW_HPP_*/