Commit 5c88774c authored by grandper's avatar grandper
Browse files

Merge branch '201202_LinearSolverPrec'

parents 41b17c9b 88ba1540
......@@ -13,6 +13,7 @@ SET(algorithm_HEADERS
algorithm/EigenSolver.hpp
algorithm/LinearSolver.hpp
algorithm/PreconditionerML.hpp
algorithm/PreconditionerLinearSolver.hpp
CACHE INTERNAL "")
SET(algorithm_SOURCES
......@@ -25,6 +26,7 @@ SET(algorithm_SOURCES
algorithm/SolverAztecOO.cpp
algorithm/EigenSolver.cpp
algorithm/LinearSolver.cpp
algorithm/PreconditionerLinearSolver.cpp
CACHE INTERNAL "")
......
......@@ -119,7 +119,7 @@ LinearSolver::solve( vectorPtr_Type solutionPtr )
return -1;
}
// Setup the Solver Operator?? Really here??
// Setup the Solver Operator
setupSolverOperator();
// Reset status informations
......@@ -184,19 +184,15 @@ LinearSolver::solve( vectorPtr_Type solutionPtr )
if( M_quitOnFailure && failure )
exit( -1 );
// Reset the solver to free the internal pointers
M_solverOperator->resetSolver();
// If the number of iterations reaches the threshold of maxIterForReuse
// we reset the preconditioners to force to solver to recompute it next
// time
if( numIters > M_maxItersForReuse )
resetPreconditioner();
// <!-- TO BE RECODED IF POSSIBLE
// AztecOO and Belos contain pointers
// to some operators.
// ML is crashing for this reason.
M_solverOperator.reset();
// -->
return numIters;
}
......@@ -287,7 +283,9 @@ void
LinearSolver::resetPreconditioner()
{
if( M_preconditioner )
{
M_preconditioner->resetPreconditioner();
}
}
bool
......@@ -319,6 +317,48 @@ LinearSolver::showMe( std::ostream& output ) const
}
}
void
LinearSolver::setupSolverOperator()
{
// Creation of a solver if there exists any
if( !M_solverOperator )
{
switch( M_solverType )
{
case Belos:
M_solverOperator.reset( Operators::SolverOperatorFactory::instance().createObject( "Belos" ) );
break;
case AztecOO:
M_solverOperator.reset( Operators::SolverOperatorFactory::instance().createObject( "AztecOO" ) );
break;
default:
M_displayer->leaderPrint( "SLV- ERROR: The type of solver is not recognized!\n" );
exit( 1 );
break;
}
}
// Set the preconditioner operator in the SolverOperator object
if( M_preconditioner )
{
M_solverOperator->setPreconditioner( M_preconditioner->preconditionerPtr() );
}
else
{
M_solverOperator->setPreconditioner( M_preconditionerOperator );
}
// Set the operator in the SolverOperator object
M_solverOperator->setOperator( M_operator );
// Set the tolerance if it has been set
if( M_tolerance > 0 )
M_solverOperator->setTolerance( M_tolerance );
// Set the parameter inside the solver
M_solverOperator->setParameters( M_parameterList.sublist( "Solver: Operator List" ) );
}
// ===================================================
// Set Methods
// ===================================================
......@@ -365,14 +405,18 @@ LinearSolver::setPreconditioner( preconditionerPtr_Type preconditionerPtr )
void
LinearSolver::setPreconditioner( operatorPtr_Type preconditionerPtr )
{
// Does the solverOperator exists?
// If a LifeV::Preconditioner exists it must be deleted
M_preconditioner.reset();
M_preconditionerOperator = preconditionerPtr;
}
void
LinearSolver::setBaseMatrixForPreconditioner( matrixPtr_Type baseMatrixPtr )
{
M_baseMatrixForPreconditioner = baseMatrixPtr;
}
void
LinearSolver::setParameters( const Teuchos::ParameterList& list )
{
......@@ -520,44 +564,4 @@ LinearSolver::hasConverged() const
// Private Methods
// ===================================================
void
LinearSolver::setupSolverOperator()
{
// If a SolverOperator already exists we simply clean it!
if( M_solverOperator )
{
M_solverOperator.reset();
}
switch( M_solverType )
{
case Belos:
M_solverOperator.reset( Operators::SolverOperatorFactory::instance().createObject( "Belos" ) );
break;
case AztecOO:
M_solverOperator.reset( Operators::SolverOperatorFactory::instance().createObject( "AztecOO" ) );
break;
default:
M_displayer->leaderPrint( "SLV- ERROR: The type of solver is not recognized!\n" );
exit( 1 );
break;
}
// Set the operator in the SolverOperator object
M_solverOperator->setOperator( M_operator );
// Set the preconditioner operator in the SolverOperator object
if( M_preconditioner )
M_solverOperator->setPreconditioner( M_preconditioner->preconditionerPtr() );
else
M_solverOperator->setPreconditioner( M_preconditionerOperator );
// Set the tolerance if it has been set
if( M_tolerance > 0 )
M_solverOperator->setTolerance( M_tolerance );
// Set the parameter inside the solver
M_solverOperator->setParameters( M_parameterList.sublist( "Solver: Operator List", true, "" ) );
}
} // namespace LifeV
......@@ -195,6 +195,9 @@ public:
//! Print informations about the solver
void showMe( std::ostream& output = std::cout ) const;
//! Setup the solver operator to be used
void setupSolverOperator();
//@}
//! @name Set Method
......@@ -244,6 +247,12 @@ public:
*/
void setPreconditioner( operatorPtr_Type preconditionerPtr );
//! Method to set a matrix on which the preconditioner should be created
/*!
@param baseMatrixPtr matrix on which the preconditioner should be created
*/
void setBaseMatrixForPreconditioner( matrixPtr_Type baseMatrixPtr );
//! Method to setup the solver using Teuchos::ParameterList
/*!
@param list Teuchos::ParameterList object
......@@ -338,9 +347,6 @@ private:
//! @name Private Methods
//@{
//! Setup the solver operator to be used
void setupSolverOperator();
//@}
operatorPtr_Type M_operator;
......
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/*!
@file
@brief LinearSolver preconditioner
@author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
@maintainer Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
@date 17-09-2011
*/
#include <lifev/core/algorithm/PreconditionerLinearSolver.hpp>
#include <lifev/core/util/LifeDebug.hpp>
#include <lifev/core/util/LifeChrono.hpp>
#include <lifev/core/filter/GetPot.hpp>
namespace LifeV
{
// ===================================================
// Constructors & Destructor
// ===================================================
PreconditionerLinearSolver::PreconditionerLinearSolver( boost::shared_ptr<Epetra_Comm> comm ) :
Preconditioner ( comm ),
M_printSubiterationCount( false ),
M_precName ( "" ),
M_precDataSection ( "" )
{
}
PreconditionerLinearSolver::~PreconditionerLinearSolver()
{
M_solver.reset();
M_preconditioner.reset();
}
// ===================================================
// Methods
// ===================================================
void
PreconditionerLinearSolver::createParametersList( list_Type& list,
const GetPot& dataFile,
const std::string& section,
const std::string& subSection )
{
createLinearSolverList( list, dataFile, section, subSection, M_displayer.comm()->MyPID() == 0 );
}
void
PreconditionerLinearSolver::createLinearSolverList( list_Type& list,
const GetPot& dataFile,
const std::string& section,
const std::string& subsection,
const bool& verbose )
{
bool displayList = dataFile( ( section + "/displayList" ).data(), false);
// If this option is true, the solver will print the iteration count
const std::string solverParamFile = dataFile( ( section + "/" + subsection + "/parameters_file" ).data(), "none" );
list = *( Teuchos::getParametersFromXmlFile( solverParamFile ) );
if ( displayList && verbose )
{
std::cout << "PreconditionerLinearSolver parameters list:" << std::endl;
std::cout << "-----------------------------" << std::endl;
list.print( std::cout );
std::cout << "-----------------------------" << std::endl;
}
}
Int
PreconditionerLinearSolver::buildPreconditioner( operator_type& matrix )
{
// Setup the solver
M_solver.reset( new solver_Type( this->M_displayer.comm() ) );
M_solver->setParameters( M_list.sublist( "LinearSolver" ) );
M_solver->setOperator( matrix );
// Setup the preconditioner for the solver
M_preconditioner.reset( PRECFactory::instance().createObject( M_precName ) );
ASSERT( M_preconditioner.get() != 0, " Preconditioner not set" );
M_preconditioner->setDataFromGetPot( M_dataFile, M_precDataSection );
M_solver->setPreconditioner( M_preconditioner );
M_solver->buildPreconditioner();
M_solver->setupSolverOperator();
M_solver->solver()->setUsedForPreconditioning( M_printSubiterationCount );
this->M_preconditionerCreated = true;
return 0;
}
void
PreconditionerLinearSolver::resetPreconditioner()
{
M_solver.reset();
this->M_preconditionerCreated = false;
}
Real
PreconditionerLinearSolver::condest()
{
return 0.0;
}
void
PreconditionerLinearSolver::showMe( std::ostream& output ) const
{
M_solver->showMe(output);
}
// ===================================================
// Epetra Operator Interface Methods
// ===================================================
Int
PreconditionerLinearSolver::SetUseTranspose( const bool useTranspose )
{
return M_solver->solver()->SetUseTranspose(useTranspose);
}
bool
PreconditionerLinearSolver::UseTranspose()
{
return M_solver->solver()->UseTranspose();
}
Int
PreconditionerLinearSolver::Apply( const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const
{
return M_solver->solver()->Apply( X, Y );
}
Int
PreconditionerLinearSolver::ApplyInverse( const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const
{
if( M_solver )
{
M_solver->solver()->ApplyInverse( X, Y );
if( M_printSubiterationCount )
{
M_displayer.leaderPrint( "> ", M_solver->numIterations(), " subiterations\n" );
}
}
return 0;
}
const Epetra_Map&
PreconditionerLinearSolver::OperatorRangeMap() const
{
return M_solver->solver()->OperatorRangeMap();
}
const Epetra_Map&
PreconditionerLinearSolver::OperatorDomainMap() const
{
return M_solver->solver()->OperatorDomainMap();
}
// ===================================================
// Set Methods
// ===================================================
void
PreconditionerLinearSolver::setDataFromGetPot ( const GetPot& dataFile, const std::string& section )
{
createLinearSolverList( M_list, dataFile, section, "LinearSolver", M_displayer.comm()->MyPID() == 0 );
M_printSubiterationCount = this->M_list.get( "Print Subiteration Count", false );
M_precName = this->M_list.get( "Preconditioner", "ML" );
M_precDataSection = this->M_list.get( "Preconditioner Data Section", "" );
M_dataFile = dataFile;
}
void
PreconditionerLinearSolver::setSolver( SolverAztecOO& /*solver*/ )
{
}
// ===================================================
// Get Methods
// ===================================================
//! Return true if the preconditioner is set
bool
PreconditionerLinearSolver::isPreconditionerSet() const
{
return M_solver;
}
PreconditionerLinearSolver::prec_raw_type*
PreconditionerLinearSolver::preconditioner()
{
return M_solver->solver().get();
}
PreconditionerLinearSolver::prec_type
PreconditionerLinearSolver::preconditionerPtr()
{
return M_solver->solver();
}
PreconditionerLinearSolver::solverPtr_Type
PreconditionerLinearSolver::solverPtr()
{
return M_solver;
}
std::string
PreconditionerLinearSolver::preconditionerType()
{
return "LinearSolver preconditioner";
}
// ===================================================
// Private Methods
// ===================================================
} // namespace LifeV
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/*!
@file
@brief Preconditioner Solver Belos
@author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
@maintainer Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
@date 17-09-2011
*/
#ifndef PRECONDITIONERSOLVERBELOS_HPP
#define PRECONDITIONERSOLVERBELOS_HPP 1
// Tell the compiler to ignore specific kind of warnings:
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"
#include <Teuchos_ParameterList.hpp>
#include "Teuchos_XMLParameterListHelpers.hpp"
#include <Epetra_Operator.h>
// Tell the compiler to ignore specific kind of warnings:
#pragma GCC diagnostic warning "-Wunused-variable"
#pragma GCC diagnostic warning "-Wunused-parameter"
#include <lifev/core/array/MapEpetra.hpp>
#include <lifev/core/array/VectorEpetra.hpp>
#include <lifev/core/array/MatrixEpetra.hpp>
#include <lifev/core/util/Displayer.hpp>
#include <lifev/core/algorithm/Preconditioner.hpp>
#include <lifev/core/algorithm/LinearSolver.hpp>
#include <lifev/core/filter/GetPot.hpp>
class GetPot;
namespace LifeV
{
//! PreconditionerLinearSolver - Class to wrap linear solver
/*!
@author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
*/
class PreconditionerLinearSolver
: public Preconditioner
{
public:
//! @name Public Types
//@{
typedef Epetra_Operator prec_raw_type;
typedef boost::shared_ptr<prec_raw_type> prec_type;
typedef LinearSolver solver_Type;
typedef boost::shared_ptr<solver_Type> solverPtr_Type;
typedef MatrixEpetra<Real> operator_raw_type;
typedef boost::shared_ptr<operator_raw_type> operator_type;
typedef Preconditioner preconditioner_Type;
typedef boost::shared_ptr<preconditioner_Type> preconditionerPtr_Type;
typedef Displayer::comm_Type comm_Type;
typedef Displayer::commPtr_Type commPtr_Type;
typedef Teuchos::ParameterList list_Type;
//@}
//! @name Constructors & Destructor
//@{
//! Default constructor
/*!
* @param comm The communicator.
*/
#ifdef HAVE_MPI
PreconditionerLinearSolver( boost::shared_ptr<Epetra_Comm> comm = boost::shared_ptr<Epetra_Comm>( new Epetra_MpiComm( MPI_COMM_WORLD ) ) );
#else
PreconditionerLinearSolver( boost::shared_ptr<Epetra_Comm> comm = boost::shared_ptr<Epetra_Comm>( new Epetra_SerialComm ) );
#endif
//! Destructor
virtual ~PreconditionerLinearSolver();
//@}
//! @name Methods
//@{
//! Create the list of parameters of the preconditioner
/*!
@param list A Parameter list to be filled
@param dataFile A GetPot object containing the data about the preconditioner
@param section The section in "dataFile" where to find data about the preconditioner
@param subSection The subsection in "dataFile" where to find data about the preconditioner
*/
virtual void createParametersList( list_Type& list,
const GetPot& dataFile,
const std::string& section,
const std::string& subSection );
//! Create the list of parameters of the preconditioner
/*!
@param list A Parameter list to be filled
@param dataFile A GetPot object containing the data about the preconditioner
@param section The section in "dataFile" where to find data about the preconditioner
@param subSection The subsection in "dataFile" where to find data about the preconditioner
*/
static void createLinearSolverList( list_Type& list,
const GetPot& dataFile,
const std::string& section,
const std::string& subSection = "SolverLinear",
const bool& verbose = true );
//! Build a preconditioner based on the given matrix
/*!
@param matrix Matrix upon which construct the preconditioner
*/
virtual Int buildPreconditioner( operator_type& matrix );
//! Reset the preconditioner
virtual void resetPreconditioner();
//! Return An estimation of the condition number of the preconditioner
virtual Real condest();
//! Show informations about the preconditioner
virtual void showMe( std::ostream& output = std::cout ) const;
//@}
//! @name Epetra Operator Interface Methods
//@{