PreconditionerBlockJacobi.hpp

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2018-2019 Total, S.A
 * Copyright (c) 2019-     GEOSX Contributors
 * All right reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */

#ifndef GEOSX_LINEARALGEBRA_SOLVERS_PRECONDITIONERBLOCKJACOBI_HPP_
#define GEOSX_LINEARALGEBRA_SOLVERS_PRECONDITIONERBLOCKJACOBI_HPP_

#include "linearAlgebra/interfaces/LinearOperator.hpp"
#include "linearAlgebra/solvers/PreconditionerBase.hpp"
#include "linearAlgebra/interfaces/BlasLapackLA.hpp"

namespace geosx
{

template< typename LAI >
class PreconditionerBlockJacobi : public PreconditionerBase< LAI >
{
public:

  using Base = PreconditionerBase< LAI >;

  using Vector = typename Base::Vector;

  using Matrix = typename Base::Matrix;

  PreconditionerBlockJacobi( localIndex const & blockSize = 0 )
    : m_blockDiag{}
  {
    m_blockSize = blockSize;
  }

  virtual void compute( Matrix const & mat ) override
  {
    GEOSX_LAI_ASSERT( mat.ready() );
    GEOSX_LAI_ASSERT_GT( m_blockSize, 0 );
    GEOSX_LAI_ASSERT_EQ( mat.numLocalRows() % m_blockSize, 0 );
    GEOSX_LAI_ASSERT_EQ( mat.numLocalCols() % m_blockSize, 0 );

    PreconditionerBase< LAI >::compute( mat );

    m_blockDiag.createWithLocalSize( mat.numLocalRows(), mat.numLocalCols(), m_blockSize, mat.getComm() );
    m_blockDiag.open();

    array1d< globalIndex > idxBlk( m_blockSize );
    array2d< real64 > values( m_blockSize, m_blockSize );
    array2d< real64 > valuesInv( m_blockSize, m_blockSize );
    array1d< globalIndex > cols;
    array1d< real64 > vals;
    for( globalIndex i = mat.ilower(); i < mat.iupper(); i+=m_blockSize )
    {
      values.setValues< serialPolicy >( 0.0 );
      for( localIndex j = 0; j < m_blockSize; ++j )
      {
        globalIndex const iRow = i + LvArray::integerConversion< globalIndex >( j );
        idxBlk[j] = iRow;
        localIndex const rowLength = mat.globalRowLength( iRow );
        cols.resize( rowLength );
        vals.resize( rowLength );
        mat.getRowCopy( iRow, cols, vals );
        for( localIndex k = 0; k < rowLength; ++k )
        {
          localIndex const jCol = LvArray::integerConversion< localIndex >( cols[k]-i );
          if( cols[k] >= i && cols[k] < i+LvArray::integerConversion< globalIndex >( m_blockSize ) )
          {
            values( j, jCol ) = vals[k];
          }
        }
      }
      BlasLapackLA::matrixInverse( values, valuesInv );
      m_blockDiag.insert( idxBlk, idxBlk, valuesInv );
    }
    m_blockDiag.close();
  }

  virtual void compute( Matrix const & mat,
                        DofManager const & dofManager ) override
  {
    GEOSX_UNUSED_VAR( dofManager );
    compute( mat );
  }

  virtual void clear() override
  {
    m_blockDiag.reset();
  }

  virtual void apply( Vector const & src,
                      Vector & dst ) const override
  {
    GEOSX_LAI_ASSERT( m_blockDiag.ready() );
    GEOSX_LAI_ASSERT_EQ( this->numGlobalRows(), dst.globalSize() );
    GEOSX_LAI_ASSERT_EQ( this->numGlobalCols(), src.globalSize() );

    m_blockDiag.apply( src, dst );
  }

  virtual bool hasPreconditionerMatrix() const override
  {
    GEOSX_LAI_ASSERT( m_blockDiag.ready() );
    return true;
  }

  virtual Matrix const & preconditionerMatrix() const override
  {
    GEOSX_LAI_ASSERT( m_blockDiag.ready() );
    return m_blockDiag;
  }

private:

  Matrix m_blockDiag;

  localIndex m_blockSize = 0;
};

}

#endif //GEOSX_LINEARALGEBRA_SOLVERS_PRECONDITIONERBLOCKJACOBI_HPP_