Ifpack2_DiagonalFilter_def.hpp

// @HEADER
// *****************************************************************************
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//
// Copyright 2009 NTESS and the Ifpack2 contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef IFPACK2_DIAGONALFILTER_DEF_HPP
#define IFPACK2_DIAGONALFILTER_DEF_HPP
#include "Ifpack2_DiagonalFilter_decl.hpp"
#include <vector>

#include "Teuchos_Comm.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_Vector.hpp"


namespace Ifpack2 {

template<class MatrixType>
DiagonalFilter<MatrixType>::
DiagonalFilter (const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Matrix,
                typename Teuchos::ScalarTraits<Scalar>::magnitudeType AbsoluteThreshold,
                typename Teuchos::ScalarTraits<Scalar>::magnitudeType RelativeThreshold):
  A_(Matrix),
  AbsoluteThreshold_(AbsoluteThreshold),
  RelativeThreshold_(RelativeThreshold)
{
  pos_.resize(getLocalNumRows());
  val_=Teuchos::rcp(new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A_->getRowMap()));

  nonconst_local_inds_host_view_type Indices("Indices",getLocalMaxNumRowEntries());
  nonconst_values_host_view_type Values("Values",getLocalMaxNumRowEntries());
  size_t NumEntries;
  magnitudeType mysign;


  for (size_t MyRow = 0 ; MyRow < getLocalNumRows() ; ++MyRow) {
    pos_[MyRow] = -1;
    A_->getLocalRowCopy(MyRow,Indices,Values,NumEntries);

    for (size_t i = 0 ; i < NumEntries ; ++i) {
      if ((size_t)Indices[i] == MyRow) {
        pos_[MyRow] = i;
        if (Teuchos::ScalarTraits<Scalar>::real(Values[i]) < 0) {
          mysign=-Teuchos::ScalarTraits<magnitudeType>::one();
        }
        else {
          mysign=Teuchos::ScalarTraits<magnitudeType>::one();
        }
        val_->replaceLocalValue (MyRow, Values[i] * (RelativeThreshold_ - 1) +
                                 AbsoluteThreshold_ * mysign);
        break;
      }
    }
  }
}

template<class MatrixType>
DiagonalFilter<MatrixType>::~DiagonalFilter() { }

template<class MatrixType>
Teuchos::RCP<const Teuchos::Comm<int> > DiagonalFilter<MatrixType>::getComm() const
{
  return A_->getComm();
}


template<class MatrixType>
Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
                               typename MatrixType::global_ordinal_type,
                               typename MatrixType::node_type> >
DiagonalFilter<MatrixType>::getRowMap() const
{
  return A_->getRowMap();
}

template<class MatrixType>
Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
                               typename MatrixType::global_ordinal_type,
                               typename MatrixType::node_type> >
DiagonalFilter<MatrixType>::getColMap() const
{
  return A_->getColMap();
}

template<class MatrixType>
Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
                               typename MatrixType::global_ordinal_type,
                               typename MatrixType::node_type> >
DiagonalFilter<MatrixType>::getDomainMap() const
{
  return A_->getDomainMap();
}

template<class MatrixType>
Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
                               typename MatrixType::global_ordinal_type,
                               typename MatrixType::node_type> >
DiagonalFilter<MatrixType>::getRangeMap() const
{
  return A_->getRangeMap();
}

template<class MatrixType>
Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type,
                                     typename MatrixType::global_ordinal_type,
                                     typename MatrixType::node_type> >
DiagonalFilter<MatrixType>::getGraph() const
{
  return A_->getGraph();
}

template<class MatrixType>
global_size_t DiagonalFilter<MatrixType>::getGlobalNumRows() const
{
  return A_->getGlobalNumRows();
}

template<class MatrixType>
global_size_t DiagonalFilter<MatrixType>::getGlobalNumCols() const
{
  return A_->getGlobalNumCols();
}

template<class MatrixType>
size_t DiagonalFilter<MatrixType>::getLocalNumRows() const
{
  return A_->getLocalNumRows();
}

template<class MatrixType>
size_t DiagonalFilter<MatrixType>::getLocalNumCols() const
{
  return A_->getLocalNumCols();
}

template<class MatrixType>
typename MatrixType::global_ordinal_type DiagonalFilter<MatrixType>::getIndexBase() const
{
  return A_->getIndexBase();
}

template<class MatrixType>
global_size_t DiagonalFilter<MatrixType>::getGlobalNumEntries() const
{
  return A_->getGlobalNumEntries();
}

template<class MatrixType>
size_t DiagonalFilter<MatrixType>::getLocalNumEntries() const
{
  return A_->getLocalNumEntries();
}

template<class MatrixType>
size_t DiagonalFilter<MatrixType>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
{
  return A_->getNumEntriesInGlobalRow(globalRow);
}

template<class MatrixType>
size_t DiagonalFilter<MatrixType>::getNumEntriesInLocalRow(LocalOrdinal localRow) const
{
  return A_->getNumEntriesInLocalRow(localRow);
}

template<class MatrixType>
size_t DiagonalFilter<MatrixType>::getGlobalMaxNumRowEntries() const
{
  return A_->getGlobalMaxNumRowEntries();
}

template<class MatrixType>
size_t DiagonalFilter<MatrixType>::getLocalMaxNumRowEntries() const
{
  return A_->getLocalMaxNumRowEntries();
}

template<class MatrixType>
typename MatrixType::local_ordinal_type DiagonalFilter<MatrixType>::getBlockSize() const
{
  return A_->getBlockSize();
}

template<class MatrixType>
bool DiagonalFilter<MatrixType>::hasColMap() const
{
  return A_->hasColMap();
}

template<class MatrixType>
bool DiagonalFilter<MatrixType>::isLocallyIndexed() const
{
  return A_->isLocallyIndexed();
}

template<class MatrixType>
bool DiagonalFilter<MatrixType>::isGloballyIndexed() const
{
  return A_->isGloballyIndexed();
}

template<class MatrixType>
bool DiagonalFilter<MatrixType>::isFillComplete() const
{
  return A_->isFillComplete();
}

template<class MatrixType>
void DiagonalFilter<MatrixType>::
  getGlobalRowCopy (GlobalOrdinal GlobalRow,
                   nonconst_global_inds_host_view_type &Indices,
                   nonconst_values_host_view_type &Values,
                   size_t& NumEntries) const

{
  Teuchos::ArrayRCP< const Scalar > myvals=val_->get1dView();
  LocalOrdinal LocalRow=getRowMap()->getLocalElement(GlobalRow);

  A_->getGlobalRowCopy(GlobalRow, Indices,Values,NumEntries);

  if (pos_[LocalRow] != -1)
    Values[pos_[LocalRow]] += myvals[LocalRow];
}


template<class MatrixType>
void DiagonalFilter<MatrixType>::
 getLocalRowCopy (LocalOrdinal LocalRow,
                   nonconst_local_inds_host_view_type &Indices,
                   nonconst_values_host_view_type &Values,
                   size_t& NumEntries) const
{
  Teuchos::ArrayRCP< const Scalar > myvals=val_->get1dView();

  A_->getLocalRowCopy(LocalRow, Indices,Values,NumEntries);

  if (pos_[LocalRow] != -1)
    Values[pos_[LocalRow]] += myvals[LocalRow];
}


template<class MatrixType>
void DiagonalFilter<MatrixType>::getGlobalRowView(GlobalOrdinal /* GlobalRow */,
                                                  global_inds_host_view_type &/*indices*/,
                                                  values_host_view_type &/*values*/) const
{
  throw std::runtime_error("Ifpack2::DiagonalFilter: does not support getGlobalRowView.");
}


template<class MatrixType>
void DiagonalFilter<MatrixType>::getLocalRowView(LocalOrdinal /* LocalRow */,
    local_inds_host_view_type & /*indices*/,
    values_host_view_type & /*values*/) const
{
  throw std::runtime_error("Ifpack2::DiagonalFilter: does not support getLocalRowView.");
}


template<class MatrixType>
void DiagonalFilter<MatrixType>::getLocalDiagCopy(Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const
{
  // Returns the matrix's actual diagonal, rather than the "filtered" diagonal.
  // This is dubious, but it duplicates the functionality of Old Ifpack.
  return A_->getLocalDiagCopy(diag);
}

template<class MatrixType>
void DiagonalFilter<MatrixType>::leftScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* x */)
{
  throw std::runtime_error("Ifpack2::DiagonalFilter does not support leftScale.");
}

template<class MatrixType>
void DiagonalFilter<MatrixType>::rightScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* x */)
{
  throw std::runtime_error("Ifpack2::DiagonalFilter does not support rightScale.");
}

template<class MatrixType>
void DiagonalFilter<MatrixType>::
apply (const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
       Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
       Teuchos::ETransp mode,
       Scalar alpha,
       Scalar beta) const
{
  // FIXME (mfh 27 Jul 2015) If this is a "diagonal filter," why do we
  // need to apply the whole matrix???

  Scalar one = Teuchos::ScalarTraits<Scalar>::one();
  A_->apply(X,Y,mode,alpha,beta);
  Y.elementWiseMultiply(one,*val_,X,one);
}

template<class MatrixType>
bool DiagonalFilter<MatrixType>::hasTransposeApply() const
{
  return A_->hasTransposeApply();
}

template<class MatrixType>
bool DiagonalFilter<MatrixType>::supportsRowViews() const
{
  return false;
}

template<class MatrixType>
typename DiagonalFilter<MatrixType>::mag_type DiagonalFilter<MatrixType>::getFrobeniusNorm() const
{
  throw std::runtime_error("Ifpack2::DiagonalFilter does not implement getFrobeniusNorm.");
}

} // namespace Ifpack2

#define IFPACK2_DIAGONALFILTER_INSTANT(S,LO,GO,N) \
  template class Ifpack2::DiagonalFilter< Tpetra::RowMatrix<S, LO, GO, N> >;

#endif