symmetric_lumped_schur_matrix.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <kernel/base_header.hpp>
#include <kernel/util/assertion.hpp>
 
namespace FEAT
{
  namespace Global
  {
    template
    <
      typename LumpedMatrixA_,
      typename MatrixB_,
      typename MatrixD_,
      typename FilterA_
    >
    class SymmetricLumpedSchurMatrix
    {
    public:
      typedef LumpedMatrixA_ LumpedMatrixTypeA;
      typedef MatrixB_ MatrixTypeB;
      typedef MatrixD_ MatrixTypeD;
      typedef FilterA_ FilterTypeA;
 
      typedef typename LumpedMatrixA_::DataType DataType;
      typedef typename MatrixD_::VectorTypeL VectorTypeL;
      typedef typename MatrixB_::VectorTypeR VectorTypeR;
 
      typedef LumpedMatrixA_ VectorTypeML;
      typedef LumpedMatrixA_ VectorTypeMR;
 
      typedef typename MatrixD_::GateRowType GateRowType;
      typedef typename MatrixB_::GateColType GateColType;
 
      static constexpr bool is_global = true;
      static constexpr bool is_local = false;
 
      LumpedMatrixA_ inv_lumped_matrix_a;
      const MatrixB_& matrix_b;
      const MatrixD_& matrix_d;
      const FilterA_& filter_a;
 
    private:
      // These two need to be modified even when apply() (which is const) is called
      mutable VectorTypeML _vec_ml;
      mutable VectorTypeMR _vec_mr;
 
    public:
      SymmetricLumpedSchurMatrix(
        const LumpedMatrixA_& lumped_matrix_a_,
        const MatrixB_& matrix_b_,
        const MatrixD_& matrix_d_,
        const FilterA_& filter_a_) :
        inv_lumped_matrix_a(lumped_matrix_a_.clone(LAFEM::CloneMode::Layout)),
        matrix_b(matrix_b_),
        matrix_d(matrix_d_),
        filter_a(filter_a_),
        _vec_ml(lumped_matrix_a_.clone(LAFEM::CloneMode::Layout)),
        _vec_mr(lumped_matrix_a_.clone(LAFEM::CloneMode::Layout))
      {
        ASSERT(matrix_d.columns() == matrix_b.rows());
        ASSERT(matrix_d.used_elements() == matrix_b.used_elements());
 
        inv_lumped_matrix_a.component_invert(lumped_matrix_a_);
      }
 
      virtual ~SymmetricLumpedSchurMatrix()
      {
      }
 
      const Dist::Comm* get_comm() const
      {
        return inv_lumped_matrix_a.get_comm();
      }
 
      void update_lumped_a(const LumpedMatrixA_& lumped_matrix_a_)
      {
        // If these were initialized empty (as it frequently happens with Global containers), adjust the sizes
        if(_vec_ml.local().size() == Index(0))
        {
          _vec_ml.local().clone(lumped_matrix_a_.local(), LAFEM::CloneMode::Layout);
        }
        if(_vec_mr.local().size() == Index(0))
        {
          _vec_mr.local().clone(lumped_matrix_a_.local(), LAFEM::CloneMode::Layout);
        }
        if(inv_lumped_matrix_a.local().size() == Index(0))
        {
          inv_lumped_matrix_a.local().clone(lumped_matrix_a_.local(), LAFEM::CloneMode::Layout);
        }
 
        inv_lumped_matrix_a.component_invert(lumped_matrix_a_);
 
      }
 
      SymmetricLumpedSchurMatrix clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return SymmetricSchurMatrix(inv_lumped_matrix_a.clone(mode), matrix_b.clone(mode), matrix_d.clone(mode),
        filter_a.clone(mode));
      }
 
      VectorTypeL create_vector_l() const
      {
        return matrix_d.create_vector_l();
      }
 
      VectorTypeR create_vector_r() const
      {
        return matrix_b.create_vector_r();
      }
 
      Index columns() const
      {
        return matrix_b.columns();
      }
 
      Index rows() const
      {
        return matrix_d.columns();
      }
 
      Index used_elements() const
      {
        return inv_lumped_matrix_a.used_elements() + matrix_b.used_elements() + matrix_d.used_elements();
      }
 
      std::size_t bytes() const
      {
        return _vec_ml.bytes() + _vec_mr.bytes();
      }
 
      void extract_diag(VectorTypeL& diag, bool sync=true) const
      {
        // If we have a gate and it contains neighbors, we have to do it the complicated way
        if(diag.get_gate() != nullptr && !diag.get_gate()->_ranks.empty() )
        {
          diag.format();
 
          typename MatrixB_::LocalMatrixType matrix_b1(matrix_b.convert_to_1());
 
          matrix_b1.add_trace_double_mat_mult(diag.local(), matrix_d.local(), inv_lumped_matrix_a.local(), DataType(1));
        }
        else
        {
          matrix_d.local().row_norm2sqr(diag.local(), inv_lumped_matrix_a.local());
        }
 
        if(sync)
        {
          diag.sync_0();
        }
 
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x) const
      {
        matrix_b.apply(_vec_mr, x);
 
        filter_a.filter_def(_vec_mr);
        _vec_ml.component_product(_vec_mr, inv_lumped_matrix_a);
        filter_a.filter_cor(_vec_ml);
 
        matrix_d.apply(r, _vec_ml);
      }
 
      void preproc_rhs(VectorTypeML& r, const VectorTypeMR& x) const
      {
        _vec_mr.copy(x);
        filter_a.filter_def(_vec_mr);
        _vec_ml.component_product(_vec_mr, inv_lumped_matrix_a);
        filter_a.filter_cor(_vec_ml);
 
        matrix_d.apply(r, _vec_ml);
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, const DataType alpha = DataType(1)) const
      {
        // r <- y + alpha*(D A^(-1) B)*x
        matrix_b.apply(_vec_mr, x);
 
        filter_a.filter_def(_vec_mr);
        _vec_ml.component_product(_vec_mr, inv_lumped_matrix_a);
        filter_a.filter_cor(_vec_ml);
 
        matrix_d.apply(r, _vec_ml, y, alpha);
      }
 
      //LocalMatrix convert_to_1() const
      //{
      //  LocalMatrix locmat = _matrix.clone(LAFEM::CloneMode::Weak);
      //  if((_row_gate != nullptr) && (_col_gate != nullptr))
      //    synch_matrix(locmat, *_row_gate->_comm, _row_gate->_ranks, _row_gate->_mirrors, _col_gate->_mirrors);
      //  return locmat;
      //}
 
    }; // class SymmetricLumpedSchurMatrix<...>
  } // namespace Global
} // namespace FEAT