GEOS
BlockOperatorView.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 TotalEnergies
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
20 #ifndef GEOS_LINEARALGEBRA_UTILITIES_BLOCKOPERATORVIEW_HPP_
21 #define GEOS_LINEARALGEBRA_UTILITIES_BLOCKOPERATORVIEW_HPP_
22 
23 #include "codingUtilities/SFINAE_Macros.hpp"
26 #include "common/traits.hpp"
28 
29 namespace geos
30 {
31 
41 template< typename VECTOR, typename OPERATOR = LinearOperator< VECTOR > >
42 class BlockOperatorView : public LinearOperator< BlockVectorView< VECTOR > >
43 {
44 
45 public:
46 
49 
51  using Vector = typename Base::Vector;
52 
57 
63 
69 
71 
76 
84  virtual void apply( Vector const & src,
85  Vector & dst ) const override
86  {
87  for( localIndex i = 0; i < numBlockRows(); i++ )
88  {
89  dst.block( i ).zero();
90  VECTOR temp( dst.block( i ) );
91  for( localIndex j = 0; j < numBlockCols(); j++ )
92  {
93  if( m_operators( i, j ) != nullptr )
94  {
95  m_operators( i, j )->apply( src.block( j ), temp );
96  dst.block( i ).axpy( 1.0, temp );
97  }
98  }
99  }
100  }
101 
111  template< typename OP = OPERATOR >
112  std::enable_if_t< traits::VectorBasedTraits< Vector >::template HasMemberFunction_applyTranspose< OP > >
114  BlockVectorView< VECTOR > & dst ) const
115  {
116  for( localIndex j = 0; j < numBlockCols(); j++ )
117  {
118  dst.block( j ).zero();
119  VECTOR temp( dst.block( j ) );
120  for( localIndex i = 0; i < numBlockRows(); i++ )
121  {
122  if( m_operators( j, i ) != nullptr )
123  {
124  m_operators( j, i )->applyTranspose( src.block( i ), temp );
125  dst.block( j ).axpy( 1.0, temp );
126  }
127  }
128  }
129  }
130 
135  virtual globalIndex numGlobalRows() const override
136  {
137  return computeRowSize( []( OPERATOR const & block ) { return block.numGlobalRows(); } );
138  }
139 
144  virtual globalIndex numGlobalCols() const override
145  {
146  return computeColSize( []( OPERATOR const & block ) { return block.numGlobalCols(); } );
147  }
148 
154  virtual localIndex numLocalRows() const override
155  {
156  return computeRowSize( []( OPERATOR const & block ) { return block.numLocalRows(); } );
157  }
158 
164  virtual localIndex numLocalCols() const override
165  {
166  return computeColSize( []( OPERATOR const & block ) { return block.numLocalCols(); } );
167  }
168 
176  virtual MPI_Comm comm() const override
177  {
178  return m_operators( 0, 0 )->comm();
179  }
180 
182 
187 
193  {
194  return m_operators.size( 0 );
195  }
196 
202  {
203  return m_operators.size( 1 );
204  }
205 
212  OPERATOR const & block( localIndex const blockRowIndex, localIndex const blockColIndex ) const
213  {
214  GEOS_LAI_ASSERT( m_operators( blockRowIndex, blockColIndex ) != nullptr );
215  return *m_operators( blockRowIndex, blockColIndex );
216  }
217 
221  OPERATOR & block( localIndex const blockRowIndex, localIndex const blockColIndex )
222  {
223  GEOS_LAI_ASSERT( m_operators( blockRowIndex, blockColIndex ) != nullptr );
224  return *m_operators( blockRowIndex, blockColIndex );
225  }
226 
230  OPERATOR const & operator()( localIndex const blockRowIndex, localIndex const blockColIndex = 0 ) const
231  {
232  return block( blockRowIndex, blockColIndex );
233  }
234 
238  OPERATOR & operator()( localIndex const blockRowIndex, localIndex const blockColIndex = 0 )
239  {
240  return block( blockRowIndex, blockColIndex );
241  }
242 
244 
245 protected:
246 
252  BlockOperatorView( localIndex const nRows, localIndex const nCols )
253  : m_operators( nRows, nCols )
254  {
255  GEOS_LAI_ASSERT_GT( nRows, 0 );
256  GEOS_LAI_ASSERT_GT( nCols, 0 );
257  }
258 
264 
270 
277  void setPointer( localIndex const blockRowIndex, localIndex const blockColIndex, OPERATOR * op )
278  {
279  GEOS_LAI_ASSERT_GE( blockRowIndex, 0 );
280  GEOS_LAI_ASSERT_GT( numBlockRows(), blockRowIndex );
281  GEOS_LAI_ASSERT_GE( blockColIndex, 0 );
282  GEOS_LAI_ASSERT_GT( numBlockCols(), blockColIndex );
283  m_operators( blockRowIndex, blockColIndex ) = op;
284  }
285 
286 private:
287 
288  template< typename FUNC >
289  auto computeRowSize( FUNC func ) const -> decltype( func( std::declval< OPERATOR const >() ) );
290 
291  template< typename FUNC >
292  auto computeColSize( FUNC func ) const -> decltype( func( std::declval< OPERATOR const >() ) );
293 
295  array2d< OPERATOR * > m_operators;
296 };
297 
298 template< typename VECTOR, typename OPERATOR >
299 template< typename FUNC >
300 auto BlockOperatorView< VECTOR, OPERATOR >::computeRowSize( FUNC func ) const -> decltype( func( std::declval< OPERATOR const >() ) )
301 {
302  using sizeType = decltype( func( std::declval< OPERATOR const >() ) );
303  sizeType rowSize = 0;
304  for( localIndex i = 0; i < numBlockRows(); ++i )
305  {
306  for( localIndex j = 0; j < numBlockCols(); ++j )
307  {
308  if( m_operators( i, j ) != nullptr )
309  {
310  rowSize += func( block( i, j ) );
311  break;
312  }
313  }
314  }
315  return rowSize;
316 }
317 
318 template< typename VECTOR, typename OPERATOR >
319 template< typename FUNC >
320 auto BlockOperatorView< VECTOR, OPERATOR >::computeColSize( FUNC func ) const -> decltype( func( std::declval< OPERATOR const >() ) )
321 {
322  using sizeType = decltype( func( std::declval< OPERATOR const >() ) );
323  sizeType colSize = 0;
324  for( localIndex j = 0; j < numBlockCols(); j++ )
325  {
326  for( localIndex i = 0; i < numBlockRows(); ++i )
327  {
328  if( m_operators( i, j ) != nullptr )
329  {
330  colSize += func( block( i, j ) );
331  break;
332  }
333  }
334  }
335  return colSize;
336 }
337 
338 }// end geos namespace
339 
340 
341 #endif /*GEOS_LINEARALGEBRA_UTILITIES_BLOCKOPERATORVIEW_HPP_*/
Abstract view of a block operator.
BlockOperatorView & operator=(BlockOperatorView &&)=delete
Deleted move assignment.
virtual globalIndex numGlobalCols() const override
Get the number of global columns.
virtual globalIndex numGlobalRows() const override
Get the number of global rows.
virtual MPI_Comm comm() const override
Get the MPI communicator the matrix was created with.
OPERATOR const & operator()(localIndex const blockRowIndex, localIndex const blockColIndex=0) const
Get the operator corresponding to a sub-block.
BlockOperatorView(BlockOperatorView< VECTOR, OPERATOR > &&x)=default
Move constructor.
localIndex numBlockCols() const
Get number of block columns.
virtual void apply(Vector const &src, Vector &dst) const override
Apply operator to a vector.
virtual localIndex numLocalCols() const override
Get the number of local columns.
localIndex numBlockRows() const
Get number of block rows.
typename Base::Vector Vector
Alias for vector type.
BlockOperatorView(localIndex const nRows, localIndex const nCols)
Create an operator with given number of block rows and columns.
OPERATOR & operator()(localIndex const blockRowIndex, localIndex const blockColIndex=0)
Get the operator corresponding to a sub-block.
std::enable_if_t< traits::VectorBasedTraits< Vector >::template HasMemberFunction_applyTranspose< OP > > applyTranspose(BlockVectorView< VECTOR > const &src, BlockVectorView< VECTOR > &dst) const
Apply the transpose of block operator to a block vector.
BlockOperatorView(BlockOperatorView< VECTOR, OPERATOR > const &x)=default
Copy constructor.
BlockOperatorView & operator=(BlockOperatorView const &)=delete
Deleted copy assignment.
void setPointer(localIndex const blockRowIndex, localIndex const blockColIndex, OPERATOR *op)
Set/replace a pointer to a block.
virtual localIndex numLocalRows() const override
Get the number of local rows.
OPERATOR const & block(localIndex const blockRowIndex, localIndex const blockColIndex) const
Get the operator corresponding to a sub-block.
OPERATOR & block(localIndex const blockRowIndex, localIndex const blockColIndex)
Get the operator corresponding to a sub-block.
Abstract view of a block vector.
VECTOR const & block(localIndex const blockIndex) const
Get a reference to the vector corresponding to block blockRowIndex.
Abstract base class for linear operators.
BlockVectorView< VECTOR > Vector
Alias for template parameter.
#define GEOS_LAI_ASSERT_GE(lhs, rhs)
Definition: common.hpp:70
#define GEOS_LAI_ASSERT_GT(lhs, rhs)
Definition: common.hpp:63
#define GEOS_LAI_ASSERT(expr)
Definition: common.hpp:35
Array< T, 2, PERMUTATION > array2d
Alias for 2D array.
Definition: DataTypes.hpp:191
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84