InverseNormalOperator.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_INVERSENORMALOPERATOR_HPP_
#define GEOS_LINEARALGEBRA_UTILITIES_INVERSENORMALOPERATOR_HPP_
 
#include "common/LinearOperator.hpp"
#include "common/PreconditionerBase.hpp"
#include "common/traits.hpp"
 
namespace geos
{
 
namespace internal
{
 
template< typename LAI, template< typename > class SOLVER >
bool constexpr HasApplyTranspose = traits::VectorBasedTraits< typename LAI::ParallelVector >::template HasMemberFunction_applyTranspose< SOLVER< LAI > >;
 
template< typename LAI, template< typename > class SOLVER >
class ExplicitTransposeInverse
{
public:
 
  using Matrix = typename LAI::ParallelMatrix;
  using Vector = typename LAI::ParallelVector;
  using Solver = SOLVER< LAI >;
 
  ExplicitTransposeInverse( Matrix const & matrix, Solver const & solver ):
    m_transposeSolver( solver.parameters() )
  {
    matrix.transpose( m_transposeMatrix );
    m_transposeSolver.setup( m_transposeMatrix );
  }
 
  void apply( Solver const &, Vector const & src, Vector & dst ) const
  {
    m_transposeSolver.apply( src, dst );
  }
 
public:
 
  Matrix m_transposeMatrix;
  Solver m_transposeSolver;
};
 
template< typename LAI, template< typename > class SOLVER >
class ImplicitTransposeInverse
{
protected:
 
  using Matrix = typename LAI::ParallelMatrix;
  using Vector = typename LAI::ParallelVector;
  using Solver = SOLVER< LAI >;
 
  ImplicitTransposeInverse( Matrix const &, Solver const & )
  {}
 
  void apply( Solver const & solver, Vector const & src, Vector & dst ) const
  {
    solver.applyTranspose( src, dst );
  }
};
 
template< typename LAI, template< typename > class SOLVER >
using TransposeInverse = std::conditional_t< HasApplyTranspose< LAI, SOLVER >,
                                             ImplicitTransposeInverse< LAI, SOLVER >,
                                             ExplicitTransposeInverse< LAI, SOLVER > >;
 
}
 
template< typename LAI, template< typename > class SOLVER >
class InverseNormalOperator : public LinearOperator< typename LAI::ParallelVector >
{
public:
 
  using Base = LinearOperator< typename LAI::ParallelVector >;
 
  using Matrix = typename LAI::ParallelMatrix;
 
  using Vector = typename Base::Vector;
 
  using Solver = SOLVER< LAI >;
 
  explicit InverseNormalOperator( Matrix const & mat,
                                  Solver const & solver ):
    m_matrix( mat ),
    m_solver( solver ),
    m_transposeInverse( mat, solver )
  {}
 
  virtual ~InverseNormalOperator() override = default;
 
  void apply( Vector const & src, Vector & dst ) const override
  {
    m_solver.apply( src, dst );
    m_transposeInverse.apply( m_solver, dst, dst );
  }
 
  globalIndex numGlobalRows() const override
  {
    return m_matrix.numGlobalCols();
  }
 
  globalIndex numGlobalCols() const override
  {
    return m_matrix.numGlobalCols();
  }
 
  localIndex numLocalRows() const override
  {
    return m_matrix.numLocalCols();
  }
 
  localIndex numLocalCols() const override
  {
    return m_matrix.numLocalCols();
  }
 
  MPI_Comm comm() const override
  {
    return m_matrix.comm();
  }
 
private:
 
  using TransposeInverse = internal::TransposeInverse< LAI, SOLVER >;
 
  Matrix const & m_matrix;
 
  Solver const & m_solver;
 
  TransposeInverse m_transposeInverse;
};
 
} // namespace geos
 
#endif //GEOS_LINEARALGEBRA_UTILITIES_INVERSENORMALOPERATOR_HPP_