BlockOperatorView.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 Total, S.A
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_LINEARALGEBRA_UTILITIES_BLOCKOPERATORVIEW_HPP_
#define GEOS_LINEARALGEBRA_UTILITIES_BLOCKOPERATORVIEW_HPP_
 
#include "codingUtilities/SFINAE_Macros.hpp"
#include "linearAlgebra/common/LinearOperator.hpp"
#include "linearAlgebra/utilities/BlockVector.hpp"
#include "common/traits.hpp"
#include "linearAlgebra/common/common.hpp"
 
namespace geos
{
 
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;
 
 
 
  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< traits::VectorBasedTraits< Vector >::template 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(); } );
  }
 
  virtual localIndex numLocalRows() const override
  {
    return computeRowSize( []( OPERATOR const & block ) { return block.numLocalRows(); } );
  }
 
  virtual localIndex numLocalCols() const override
  {
    return computeColSize( []( OPERATOR const & block ) { return block.numLocalCols(); } );
  }
 
  virtual MPI_Comm comm() const override
  {
    return m_operators( 0, 0 )->comm();
  }
 
 
 
  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
  {
    GEOS_LAI_ASSERT( m_operators( blockRowIndex, blockColIndex ) != nullptr );
    return *m_operators( blockRowIndex, blockColIndex );
  }
 
  OPERATOR & block( localIndex const blockRowIndex, localIndex const blockColIndex )
  {
    GEOS_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 )
  {
    GEOS_LAI_ASSERT_GT( nRows, 0 );
    GEOS_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 )
  {
    GEOS_LAI_ASSERT_GE( blockRowIndex, 0 );
    GEOS_LAI_ASSERT_GT( numBlockRows(), blockRowIndex );
    GEOS_LAI_ASSERT_GE( blockColIndex, 0 );
    GEOS_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 geos namespace
 
 
#endif /*GEOS_LINEARALGEBRA_UTILITIES_BLOCKOPERATORVIEW_HPP_*/