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 
00040 template<class Scalar,class LocalOrdinal,class GlobalOrdinal, class Node>
00041 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00042 Condest(const Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>& TIFP,
00043                const Ifpack2::CondestType CT,
00044                const int MaxIters = 1550,
00045                const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& Tol = 1e-9,
00046                const Teuchos::Ptr<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &matrix_in = Teuchos::null)
00047 {
00048   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
00049   magnitudeType ConditionNumberEstimate = -1.0;
00050   Teuchos::Ptr<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > matrix = matrix_in;
00051   if (matrix_in == Teuchos::null) {
00052     matrix = TIFP.getMatrix().ptr();
00053   }
00054 
00055   if (CT == Ifpack2::Cheap) {
00056 
00057     // Create a vector with all values equal to one
00058     Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Ones(TIFP.getDomainMap());
00059     Ones.putScalar(1.0);
00060     // Create the vector of results
00061     Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> OnesResult(TIFP.getRangeMap());
00062     // Compute the effect of the solve on the vector of ones
00063     TIFP.apply(Ones, OnesResult);
00064     ConditionNumberEstimate = OnesResult.normInf();
00065   }
00066   else if (CT == Ifpack2::CG) {
00067 
00068 #ifdef HAVE_IFPACK2_AZTECOO
00069 this code not yet converted!!!
00070     Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> LHS(IFP.getDomainMap());
00071     LHS.PutScalar(0.0);
00072     Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> RHS(IFP.getRangeMap());
00073     RHS.Random();
00074     Tpetra_LinearProblem Problem;
00075     Problem.SetOperator(matrix);
00076     Problem.SetLHS(&LHS);
00077     Problem.SetRHS(&RHS);
00078 
00079     AztecOO Solver(Problem);
00080     Solver.SetAztecOption(AZ_output,AZ_none);
00081     Solver.SetAztecOption(AZ_solver,AZ_cg_condnum);
00082     Solver.Iterate(MaxIters,Tol);
00083 
00084     const double* status = Solver.GetAztecStatus();
00085     ConditionNumberEstimate = status[AZ_condnum];
00086 #endif
00087 
00088   } else if (CT == Ifpack2::GMRES) {
00089 
00090 #ifdef HAVE_IFPACK2_AZTECOO
00091 this code not yet converted!!!
00092     Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> LHS(IFP.getDomainMap());
00093     LHS.PutScalar(0.0);
00094     Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> RHS(IFP.getRangeMap());
00095     RHS.Random();
00096     Tpetra_LinearProblem Problem;
00097     Problem.SetOperator(matrix);
00098     Problem.SetLHS(&LHS);
00099     Problem.SetRHS(&RHS);
00100 
00101     AztecOO Solver(Problem);
00102     Solver.SetAztecOption(AZ_solver,AZ_gmres_condnum);
00103     Solver.SetAztecOption(AZ_output,AZ_none);
00104     // the following can be problematic for large problems,
00105     // but any restart would destroy useful information about
00106     // the condition number.
00107     Solver.SetAztecOption(AZ_kspace,MaxIters);
00108     Solver.Iterate(MaxIters,Tol);
00109 
00110     const double* status = Solver.GetAztecStatus();
00111     ConditionNumberEstimate = status[AZ_condnum];
00112 #endif
00113   }
00114 
00115   return(ConditionNumberEstimate);
00116 }
00117 
00118 }//namespace Ifpack2
00119 
00120 #endif // IFPACK2_CONDEST_HPP
00121 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends