Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_Condest.hpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#ifndef IFPACK2_CONDEST_HPP
44#define IFPACK2_CONDEST_HPP
45
46#include "Ifpack2_ConfigDefs.hpp"
47#include "Ifpack2_CondestType.hpp"
49#include <Teuchos_Ptr.hpp>
50
51namespace Ifpack2 {
52
92template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93typename Teuchos::ScalarTraits<Scalar>::magnitudeType
95 const Ifpack2::CondestType CT,
96 const int MaxIters = 1550,
97 const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& Tol = Teuchos::as<Scalar> (1e-9),
98 const Teuchos::Ptr<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& matrix_in = Teuchos::null)
99{
100 using Teuchos::Ptr;
101 typedef Teuchos::ScalarTraits<Scalar> STS;
102 typedef typename STS::magnitudeType MT;
103 typedef Teuchos::ScalarTraits<MT> STM;
104 typedef Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> row_matrix_type;
105 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> vec_type;
106
107 MT condNumEst = -STS::magnitude( STS::one() );
108
109 // Users may either provide a matrix for which to estimate the
110 // condition number, or use the Preconditioner's built-in matrix.
111 Ptr<const row_matrix_type> matrix = matrix_in;
112 if (matrix_in == Teuchos::null) {
113 matrix = TIFP.getMatrix ().ptr ();
114 TEUCHOS_TEST_FOR_EXCEPTION(
115 matrix == Teuchos::null,
116 std::logic_error,
117 "Ifpack2::Condest: Both the input matrix (matrix_in) and the Ifpack2 "
118 "preconditioner's matrix are null, so we have no matrix with which to "
119 "compute a condition number estimate. This probably indicates a bug "
120 "in Ifpack2, since no Ifpack2::Preconditioner subclass should accept a"
121 "null matrix.");
122 }
123
124 if (CT == Ifpack2::Cheap) {
125 vec_type ones (TIFP.getDomainMap ()); // Vector of ones
126 ones.putScalar (STS::one ());
127 vec_type onesResult (TIFP.getRangeMap ()); // A*ones
128 onesResult.putScalar (STS::zero ());
129 TIFP.apply (ones, onesResult);
130 condNumEst = onesResult.normInf (); // max (abs (A*ones))
131 TEUCHOS_TEST_FOR_EXCEPTION(
132 STM::isnaninf (condNumEst),
133 std::runtime_error,
134 "Ifpack2::Condest: $\\|A*[1, ..., 1]^T\\|_{\\infty}$ = " << condNumEst << " is NaN or Inf.");
135 } else if (CT == Ifpack2::CG) {
136 TEUCHOS_TEST_FOR_EXCEPTION(
137 true, std::logic_error,
138 "Ifpack2::Condest: Condition number estimation using CG is not currently supported.");
139 } else if (CT == Ifpack2::GMRES) {
140 TEUCHOS_TEST_FOR_EXCEPTION(
141 true, std::logic_error,
142 "Ifpack2::Condest: Condition number estimation using GMRES is not currently supported.");
143 }
144 return condNumEst;
145}
146
147}//namespace Ifpack2
148
149#endif // IFPACK2_CONDEST_HPP
150
Interface for all Ifpack2 preconditioners.
Definition Ifpack2_Preconditioner.hpp:108
virtual Teuchos::RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getMatrix() const =0
The input matrix given to the constructor.
virtual void apply(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const =0
Apply the preconditioner to X, putting the result in Y.
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
The range Map of this operator.
virtual Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
The domain Map of this operator.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74
CondestType
Ifpack2::CondestType: enum to define the type of condition number estimate.
Definition Ifpack2_CondestType.hpp:50
@ CG
Uses AztecOO's CG.
Definition Ifpack2_CondestType.hpp:52
@ Cheap
cheap estimate
Definition Ifpack2_CondestType.hpp:51
@ GMRES
Uses AztecOO's GMRES.
Definition Ifpack2_CondestType.hpp:53