Ifpack2_Condest.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_CONDEST_HPP
#define IFPACK2_CONDEST_HPP

#include "Ifpack2_ConfigDefs.hpp"
#include "Ifpack2_CondestType.hpp"
#include "Ifpack2_Preconditioner.hpp"
#include <Teuchos_Ptr.hpp>

namespace Ifpack2 {

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Condest (const Ifpack2::Preconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>& TIFP,
         const Ifpack2::CondestType CT,
         const int MaxIters = 1550,
         const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& Tol = Teuchos::as<Scalar> (1e-9),
         const Teuchos::Ptr<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& matrix_in = Teuchos::null)
{
  using Teuchos::Ptr;
  typedef Teuchos::ScalarTraits<Scalar> STS;
  typedef typename STS::magnitudeType MT;
  typedef Teuchos::ScalarTraits<MT> STM;
  typedef Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> row_matrix_type;
  typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> vec_type;

  MT condNumEst = -STS::magnitude( STS::one() );

  // Users may either provide a matrix for which to estimate the
  // condition number, or use the Preconditioner's built-in matrix.
  Ptr<const row_matrix_type> matrix = matrix_in;
  if (matrix_in == Teuchos::null) {
    matrix = TIFP.getMatrix ().ptr ();
    TEUCHOS_TEST_FOR_EXCEPTION(
      matrix == Teuchos::null,
      std::logic_error,
      "Ifpack2::Condest: Both the input matrix (matrix_in) and the Ifpack2 "
      "preconditioner's matrix are null, so we have no matrix with which to "
      "compute a condition number estimate.  This probably indicates a bug "
      "in Ifpack2, since no Ifpack2::Preconditioner subclass should accept a"
      "null matrix.");
  }

  if (CT == Ifpack2::Cheap) {
    vec_type ones (TIFP.getDomainMap ()); // Vector of ones
    ones.putScalar (STS::one ());
    vec_type onesResult (TIFP.getRangeMap ()); // A*ones
    onesResult.putScalar (STS::zero ());
    TIFP.apply (ones, onesResult);
    condNumEst = onesResult.normInf (); // max (abs (A*ones))
    TEUCHOS_TEST_FOR_EXCEPTION(
      STM::isnaninf (condNumEst),
      std::runtime_error,
      "Ifpack2::Condest: $\\|A*[1, ..., 1]^T\\|_{\\infty}$ = " << condNumEst << " is NaN or Inf.");
  } else if (CT == Ifpack2::CG) {
    TEUCHOS_TEST_FOR_EXCEPTION(
      true, std::logic_error,
      "Ifpack2::Condest: Condition number estimation using CG is not currently supported.");
  } else if (CT == Ifpack2::GMRES) {
    TEUCHOS_TEST_FOR_EXCEPTION(
      true, std::logic_error,
      "Ifpack2::Condest: Condition number estimation using GMRES is not currently supported.");
  }
  return condNumEst;
}

}//namespace Ifpack2

#endif // IFPACK2_CONDEST_HPP