Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Condest.hpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #ifndef IFPACK2_CONDEST_HPP
00031 #define IFPACK2_CONDEST_HPP
00032 
00033 #include "Ifpack2_ConfigDefs.hpp"
00034 #include "Ifpack2_CondestType.hpp"
00035 #include "Ifpack2_Preconditioner.hpp"
00036 #include <Teuchos_Ptr.hpp>
00037 
00038 namespace Ifpack2 {
00039 
00071 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00072 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00073 Condest (const Ifpack2::Preconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>& TIFP,
00074    const Ifpack2::CondestType CT,
00075    const int MaxIters = 1550,
00076    const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& Tol = Teuchos::as<Scalar> (1e-9),
00077    const Teuchos::Ptr<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& matrix_in = Teuchos::null)
00078 {
00079   using Teuchos::Ptr;
00080   typedef Teuchos::ScalarTraits<Scalar> STS;
00081   typedef typename STS::magnitudeType MT;
00082   typedef Teuchos::ScalarTraits<MT> STM;
00083   typedef Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> row_matrix_type;
00084   typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> vec_type;
00085 
00086   MT condNumEst = -STS::magnitude( STS::one() );
00087 
00088   // Users may either provide a matrix for which to estimate the
00089   // condition number, or use the Preconditioner's built-in matrix.
00090   Ptr<const row_matrix_type> matrix = matrix_in;
00091   if (matrix_in == Teuchos::null) {
00092     matrix = TIFP.getMatrix ().ptr ();
00093     TEUCHOS_TEST_FOR_EXCEPTION(
00094       matrix == Teuchos::null,
00095       std::logic_error,
00096       "Ifpack2::Condest: Both the input matrix (matrix_in) and the Ifpack2 "
00097       "preconditioner's matrix are null, so we have no matrix with which to "
00098       "compute a condition number estimate.  This probably indicates a bug "
00099       "in Ifpack2, since no Ifpack2::Preconditioner subclass should accept a"
00100       "null matrix.");
00101   }
00102 
00103   if (CT == Ifpack2::Cheap) {
00104     vec_type ones (TIFP.getDomainMap ()); // Vector of ones
00105     ones.putScalar (STS::one ());
00106     vec_type onesResult (TIFP.getRangeMap ()); // A*ones
00107     onesResult.putScalar (STS::zero ());
00108     TIFP.apply (ones, onesResult);
00109     condNumEst = onesResult.normInf (); // max (abs (A*ones))
00110     TEUCHOS_TEST_FOR_EXCEPTION(
00111       STM::isnaninf (condNumEst),
00112       std::runtime_error,
00113       "Ifpack2::Condest: $\\|A*[1, ..., 1]^T\\|_{\\infty}$ = " << condNumEst << " is NaN or Inf.");
00114   } else if (CT == Ifpack2::CG) {
00115     TEUCHOS_TEST_FOR_EXCEPTION(
00116       true, std::logic_error, 
00117       "Ifpack2::Condest: Condition number estimation using CG is not currently supported.");
00118   } else if (CT == Ifpack2::GMRES) {
00119     TEUCHOS_TEST_FOR_EXCEPTION(
00120       true, std::logic_error, 
00121       "Ifpack2::Condest: Condition number estimation using GMRES is not currently supported.");
00122   }
00123   return condNumEst;
00124 }
00125 
00126 }//namespace Ifpack2
00127 
00128 #endif // IFPACK2_CONDEST_HPP
00129 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends