Teko_InverseLibrary.cpp

/*
// @HEADER
// 
// ***********************************************************************
// 
//      Teko: A package for block and physics based preconditioning
//                  Copyright 2010 Sandia Corporation 
//  
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//  
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//  
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//  
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//  
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission. 
//  
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//  
// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
// 
// ***********************************************************************
// 
// @HEADER
 
*/
 
#include "Teko_InverseLibrary.hpp"
 
#include "Teko_SolveInverseFactory.hpp"
#include "Teko_PreconditionerInverseFactory.hpp"
#include "Teko_BlockPreconditionerFactory.hpp"
 
#include "Teko_NeumannSeriesPreconditionerFactory.hpp"
#include "Teuchos_AbstractFactoryStd.hpp"
#include "Teko_Utilities.hpp"
 
#include <algorithm>
 
using Teuchos::RCP;
using Teuchos::rcp;
 
namespace Teko {

void addToStratimikosBuilder(const RCP<Stratimikos::DefaultLinearSolverBuilder> & builder)
{
   typedef Thyra::PreconditionerFactoryBase<double> PrecFactory;
 
   RCP<const Teuchos::ParameterList> parameters = builder->getValidParameters();
 
   if(!parameters->sublist("Preconditioner Types").isSublist("Neumann Series")) {
     RCP<const Teuchos::AbstractFactory<Thyra::PreconditionerFactoryBase<double> > > factory;
     
     factory = Teuchos::abstractFactoryStd<PrecFactory,Teko::NeumannSeriesPreconditionerFactory<double> >();
     builder->setPreconditioningStrategyFactory(factory,"Neumann Series");
   }
 
}
 
InverseLibrary::InverseLibrary()
{
   Teko_DEBUG_SCOPE("InverseLibrary::InverseLibrary", 10);
 
   defaultBuilder_ = Teuchos::rcp(new Stratimikos::DefaultLinearSolverBuilder());
   addToStratimikosBuilder(defaultBuilder_);
 
   // setup some valid Stratimikos parameters
 
   // set valid solve factory names
   stratValidSolver_.push_back("Belos"); 
   stratValidSolver_.push_back("Amesos"); 
   stratValidSolver_.push_back("Amesos2"); 
   stratValidSolver_.push_back("AztecOO"); 
 
   // set valid preconditioner factory name
   stratValidPrecond_.push_back("ML"); 
   stratValidPrecond_.push_back("Ifpack"); 
   stratValidPrecond_.push_back("Neumann Series"); 
   stratValidPrecond_.push_back("MueLu"); 
   stratValidPrecond_.push_back("Ifpack2"); 
 
   // set valid Teko preconditioner factory names
   PreconditionerFactory::getPreconditionerFactoryNames(blockValidPrecond_);
 
   Teko_DEBUG_MSG_BEGIN(10)
      DEBUG_STREAM << "Loaded \"block\" preconditioners = ";
      for(std::size_t i=0;i<blockValidPrecond_.size();i++)
         DEBUG_STREAM << blockValidPrecond_[i] << ", ";
      DEBUG_STREAM << std::endl;
   Teko_DEBUG_MSG_END()
}
 
InverseLibrary::InverseLibrary(const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> & strat)
  : defaultBuilder_(strat)
{
   Teko_DEBUG_SCOPE("InverseLibrary::InverseLibrary", 10);
 
   RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList(*defaultBuilder_->getValidParameters()));
   Teuchos::ParameterList lst(pl->sublist("Linear Solver Types"));
   Teuchos::ParameterList pft(pl->sublist("Preconditioner Types"));
 
   Teuchos::ParameterList::ConstIterator itr;
 
   // set valid solve factory names
   for(itr=lst.begin();itr!=lst.end();++itr)
      stratValidSolver_.push_back(itr->first); 
 
   Teko_DEBUG_MSG_BEGIN(10)
      DEBUG_STREAM << "Loaded \"Stratimikos\" solvers = ";
      for(std::size_t i=0;i<stratValidSolver_.size();i++)
         DEBUG_STREAM << stratValidSolver_[i] << ", ";
      DEBUG_STREAM << std::endl;
   Teko_DEBUG_MSG_END()
 
   // set valid prec factory names
   for(itr=pft.begin();itr!=pft.end();++itr)
      stratValidPrecond_.push_back(itr->first); 
 
   Teko_DEBUG_MSG_BEGIN(10)
      DEBUG_STREAM << "Loaded \"Stratimikos\" preconditioners = ";
      for(std::size_t i=0;i<stratValidPrecond_.size();i++)
         DEBUG_STREAM << stratValidPrecond_[i] << ", ";
      DEBUG_STREAM << std::endl;
   Teko_DEBUG_MSG_END()
 
   // set valid Teko preconditioner factory names
   PreconditionerFactory::getPreconditionerFactoryNames(blockValidPrecond_);
 
   Teko_DEBUG_MSG_BEGIN(10)
      DEBUG_STREAM << "Loaded \"block\" preconditioners = ";
      for(std::size_t i=0;i<blockValidPrecond_.size();i++)
         DEBUG_STREAM << blockValidPrecond_[i] << ", ";
      DEBUG_STREAM << std::endl;
   Teko_DEBUG_MSG_END()
}

void InverseLibrary::addInverse(const std::string & label,const Teuchos::ParameterList & pl)
{
   // strip out the label
   const std::string type = pl.get<std::string>("Type");
 
   // copy the parameter list so we can modify it
   Teuchos::ParameterList settingsList;
   settingsList.set(type,pl);
   settingsList.sublist(type).remove("Type");
 
   // is this a Stratimikos preconditioner or solver
   if(std::find(stratValidPrecond_.begin(),stratValidPrecond_.end(),type)!=stratValidPrecond_.end()) {
      // this is a Stratimikos preconditioner factory
      addStratPrecond(label,type,settingsList);
   }
   else if(std::find(stratValidSolver_.begin(),stratValidSolver_.end(),type)!=stratValidSolver_.end()) {
      // this is a Stratimikos preconditioner factory
      addStratSolver(label,type,settingsList);
   }
   else if(std::find(blockValidPrecond_.begin(),blockValidPrecond_.end(),type)!=blockValidPrecond_.end()) {
      // this is a Teko preconditioner factory
      addBlockPrecond(label,type,settingsList);
   }
   else {
      Teuchos::FancyOStream & os = *Teko::getOutputStream();
      os << "ERROR: Could not find inverse type \"" << type 
         << "\" required by inverse name \"" << label << "\"" << std::endl;
      TEUCHOS_ASSERT(false);
   }
}

void InverseLibrary::addStratSolver(const std::string & label,const std::string & type,const Teuchos::ParameterList & pl)
{
   // add some additional parameters onto the list
   RCP<Teuchos::ParameterList> stratList = rcp(new Teuchos::ParameterList());
   stratList->set("Linear Solver Type",type);
   stratList->set("Linear Solver Types",pl);
   stratList->set("Preconditioner Type","None");
 
   stratSolver_[label] = stratList;
}

void InverseLibrary::addStratPrecond(const std::string & label,const std::string & type,const Teuchos::ParameterList & pl)
{
   // add some additional parameters onto the list
   RCP<Teuchos::ParameterList> stratList = rcp(new Teuchos::ParameterList());
   stratList->set("Preconditioner Type",type);
   stratList->set("Preconditioner Types",pl);
 
   stratPrecond_[label] = stratList;
}

void InverseLibrary::addBlockPrecond(const std::string & label,const std::string & type,const Teuchos::ParameterList & pl)
{
   // add some additional parameters onto the list
   RCP<Teuchos::ParameterList> blockList = rcp(new Teuchos::ParameterList());
   blockList->set("Preconditioner Type",type);
   blockList->set("Preconditioner Settings",pl.sublist(type));
 
   // add the Teko preconditioner parameter list into the library
   blockPrecond_[label] = blockList;
}

Teuchos::RCP<const Teuchos::ParameterList> InverseLibrary::getParameterList(const std::string & label) const
{
   std::map<std::string,RCP<const Teuchos::ParameterList> >::const_iterator itr;
   
   // check preconditioners
   itr = stratPrecond_.find(label);
   if(itr!=stratPrecond_.end()) return itr->second;
    
   // check solvers
   itr = stratSolver_.find(label);
   if(itr!=stratSolver_.end()) return itr->second;
   
   // check solvers
   itr = blockPrecond_.find(label);
   if(itr!=blockPrecond_.end()) return itr->second;
 
   return Teuchos::null;
}

Teuchos::RCP<InverseFactory> InverseLibrary::getInverseFactory(const std::string & label) const
{
   Teko_DEBUG_SCOPE("InverseLibrary::getInverseFactory",10);
 
   std::map<std::string,RCP<const Teuchos::ParameterList> >::const_iterator itr;
 
   bool isStratSolver=false,isStratPrecond=false,isBlockPrecond=false;
 
   // is this a Stratimikos solver?
   itr = stratPrecond_.find(label);
   isStratPrecond = itr!=stratPrecond_.end();
 
   // is this a Stratimikos preconditioner?
   if(not isStratPrecond) {
      itr = stratSolver_.find(label);
      isStratSolver = itr!=stratSolver_.end();
   }
 
   // must be a "block" preconditioner
   if(not (isStratSolver || isStratPrecond)) {
      itr = blockPrecond_.find(label);
      isBlockPrecond = itr!=blockPrecond_.end();
   }
 
   Teko_DEBUG_MSG("Inverse \"" << label << "\" is of type " 
             << "strat prec = " << isStratPrecond << ", "
             << "strat solv = " << isStratSolver << ", " 
             << "block prec = " << isBlockPrecond,3);
 
   // Must be one of Strat solver, strat preconditioner, block preconditioner
   if(not (isStratSolver || isStratPrecond || isBlockPrecond)) {
      RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
 
      *out << "InverseLibrary::getInverseFactory could not find \"" << label << "\" ... aborting\n";
      *out << "Choose one of: " << std::endl; 
 
      *out << "   Stratimikos preconditioners = ";
      for(itr=stratPrecond_.begin();itr!=stratPrecond_.end();++itr)
         *out << "      \"" << itr->first << "\"\n";
      *out << std::endl;
 
      *out << "   Stratimikos solvers = ";
      for(itr=stratSolver_.begin();itr!=stratSolver_.end();++itr)
         *out << "      \"" << itr->first << "\"\n";
      *out << std::endl;
 
      *out << "   Block preconditioners = ";
      for(itr=blockPrecond_.begin();itr!=blockPrecond_.end();++itr)
         *out << "      \"" << itr->first << "\"\n";
      *out << std::endl;
 
      TEUCHOS_ASSERT(isStratSolver || isStratPrecond || isBlockPrecond);
   }
   
   RCP<const Teuchos::ParameterList> pl = itr->second;
 
   // build inverse factory
   if(isStratPrecond) {
      // remove required parameters
      RCP<Teuchos::ParameterList> plCopy = rcp(new Teuchos::ParameterList(*pl));
      std::string type = plCopy->get<std::string>("Preconditioner Type");
      RCP<Teuchos::ParameterList> xtraParams;
      if(plCopy->sublist("Preconditioner Types").sublist(type).isParameter("Required Parameters")) {
         xtraParams = rcp(new Teuchos::ParameterList(
               plCopy->sublist("Preconditioner Types").sublist(type).sublist("Required Parameters"))); 
         plCopy->sublist("Preconditioner Types").sublist(type).remove("Required Parameters"); 
      }
 
      // print some debuggin info
      Teko_DEBUG_MSG_BEGIN(10);
         DEBUG_STREAM << "Printing parameter list: " << std::endl; 
         Teko_DEBUG_PUSHTAB(); plCopy->print(DEBUG_STREAM); Teko_DEBUG_POPTAB();
 
         if(xtraParams!=Teuchos::null) {
            DEBUG_STREAM << "Printing extra parameters: " << std::endl; 
            Teko_DEBUG_PUSHTAB(); xtraParams->print(DEBUG_STREAM); Teko_DEBUG_POPTAB();
         }
      Teko_DEBUG_MSG_END();
 
      // Stratimikos::DefaultLinearSolverBuilder strat;
      // addToStratimikosBuilder(strat);
      defaultBuilder_->setParameterList(plCopy);
 
      // try to build a preconditioner factory
      RCP<Thyra::PreconditionerFactoryBase<double> > precFact = defaultBuilder_->createPreconditioningStrategy(type);
 
      // string must map to a preconditioner
      RCP<Teko::PreconditionerInverseFactory> precInvFact 
            = rcp(new PreconditionerInverseFactory(precFact,xtraParams,getRequestHandler()));
      precInvFact->setupParameterListFromRequestHandler();
      return precInvFact;
   }
   else if(isStratSolver) {
      RCP<Teuchos::ParameterList> solveList = rcp(new Teuchos::ParameterList(*pl));
      std::string type = solveList->get<std::string>("Linear Solver Type");
 
      // get preconditioner name, remove "Use Preconditioner" parameter
      Teuchos::ParameterList & solveSettings = solveList->sublist("Linear Solver Types").sublist(type);
      std::string precKeyWord = "Use Preconditioner";
      std::string precName = "None";
      if(solveSettings.isParameter(precKeyWord)) {
         precName = solveSettings.get<std::string>(precKeyWord);
         solveSettings.remove(precKeyWord);
      }
 
      // build Thyra preconditioner factory
      RCP<Thyra::PreconditionerFactoryBase<double> > precFactory;
      if(precName!="None") {
         // we will manually set the preconditioner, so set this to null
         solveList->set<std::string>("Preconditioner Type","None");
         
         // build inverse that preconditioner corresponds to
         RCP<PreconditionerInverseFactory> precInvFactory 
               = Teuchos::rcp_dynamic_cast<PreconditionerInverseFactory>(getInverseFactory(precName));
 
         // extract preconditioner factory from preconditioner _inverse_ factory
         precFactory = precInvFactory->getPrecFactory();
      }
 
      // Stratimikos::DefaultLinearSolverBuilder strat;
      // addToStratimikosBuilder(strat);
      defaultBuilder_->setParameterList(solveList);
 
      // try to build a solver factory
      RCP<Thyra::LinearOpWithSolveFactoryBase<double> > solveFact = defaultBuilder_->createLinearSolveStrategy(type);
      if(precFactory!=Teuchos::null)
         solveFact->setPreconditionerFactory(precFactory,precName);
 
      // if its around, build a InverseFactory
      return rcp(new SolveInverseFactory(solveFact));
   }
   else if(isBlockPrecond) {
      try {
         std::string type = pl->get<std::string>("Preconditioner Type");
         const Teuchos::ParameterList & settings = pl->sublist("Preconditioner Settings");
   
         // build preconditioner factory from the string
         RCP<PreconditionerFactory> precFact 
               = PreconditionerFactory::buildPreconditionerFactory(type,settings,Teuchos::rcpFromRef(*this));
    
         TEUCHOS_ASSERT(precFact!=Teuchos::null);
   
         // return the inverse factory object
         return rcp(new PreconditionerInverseFactory(precFact,getRequestHandler()));   
      }
      catch(std::exception & e) {
         RCP<Teuchos::FancyOStream> out = Teko::getOutputStream();
         
         *out << "Teko: \"getInverseFactory\" failed, Parameter List =\n";
         pl->print(*out);
 
         *out << "*** THROWN EXCEPTION ***\n";
         *out << e.what() << std::endl;
         *out << "************************\n";
         
         throw e;
      }
   }
 
   TEUCHOS_ASSERT(false);
}

void InverseLibrary::PrintAvailableInverses(std::ostream & os) const
{
   std::map<std::string,Teuchos::RCP<const Teuchos::ParameterList> >::const_iterator itr; 
 
   os << "Stratimikos Solvers: " << std::endl;
   os << "********************************" << std::endl;
   for(itr=stratSolver_.begin();itr!=stratSolver_.end();++itr) {
      os << "name = \"" << itr->first << "\"" << std::endl;
      itr->second->print(os);
      os << std::endl;
   }
 
   os << "Stratimikos Preconditioners: " << std::endl;
   os << "********************************" << std::endl;
   for(itr=stratPrecond_.begin();itr!=stratPrecond_.end();++itr) {
      os << "name = \"" << itr->first << "\"" << std::endl;
      itr->second->print(os);
      os << std::endl;
   }
 
   os << "Teko Preconditioners: " << std::endl;
   os << "********************************" << std::endl;
   for(itr=blockPrecond_.begin();itr!=blockPrecond_.end();++itr) {
      os << "name = \"" << itr->first << "\"" << std::endl;
      itr->second->print(os);
      os << std::endl;
   }
}

RCP<InverseLibrary> InverseLibrary::buildFromParameterList(const Teuchos::ParameterList & pl,bool useStratDefaults)
{
   // build from Stratimikos or allocate a new inverse library
   RCP<InverseLibrary> invLib;
   if(useStratDefaults)
      invLib = InverseLibrary::buildFromStratimikos();
   else
      invLib = rcp(new InverseLibrary());
 
   // to convert the void* like entry
   Teuchos::ParameterList * temp = 0;
 
   // loop over all entries in parameter list
   Teuchos::ParameterList::ConstIterator itr;
   for(itr=pl.begin();itr!=pl.end();++itr) {
      // get current entry
      std::string label             = itr->first;
      Teuchos::ParameterList & list = itr->second.getValue(temp);
      
      // add to library
      invLib->addInverse(label,list);
   }
   
   return invLib;
}

RCP<InverseLibrary> InverseLibrary::buildFromParameterList(const Teuchos::ParameterList & pl,
                                                           const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> & strat)
{
   // if strat is set to null, use the defaults
   if(strat==Teuchos::null)
     return buildFromParameterList(pl,true);
 
   // build from Stratimikos or allocate a new inverse library
   RCP<InverseLibrary> invLib = InverseLibrary::buildFromStratimikos(strat);
 
   // to convert the void* like entry
   Teuchos::ParameterList * temp = 0;
 
   // loop over all entries in parameter list
   Teuchos::ParameterList::ConstIterator itr;
   for(itr=pl.begin();itr!=pl.end();++itr) {
      // get current entry
      std::string label             = itr->first;
      Teuchos::ParameterList & list = itr->second.getValue(temp);
      
      // add to library
      invLib->addInverse(label,list);
   }
   
   return invLib;
}

Teuchos::RCP<InverseLibrary> InverseLibrary::buildFromStratimikos()
{
   RCP<InverseLibrary> invLib = rcp(new InverseLibrary());
 
   // get default inveres in Stratimikos
   RCP<Stratimikos::DefaultLinearSolverBuilder> strat = invLib->defaultBuilder_;
   RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList(*strat->getValidParameters()));
   Teuchos::ParameterList lst(pl->sublist("Linear Solver Types"));
   Teuchos::ParameterList pft(pl->sublist("Preconditioner Types"));
 
   Teuchos::ParameterList::ConstIterator itr;
   Teuchos::ParameterList * temp = 0;
 
   // loop over all entries in solver list
   for(itr=lst.begin();itr!=lst.end();++itr) {
      // get current entry
      std::string label             = itr->first;
      Teuchos::ParameterList & list = itr->second.getValue(temp);
      list.set("Type",label);
      
      // add to library
      invLib->addInverse(label,list);
   }
 
   // loop over all entries in preconditioner list
   for(itr=pft.begin();itr!=pft.end();++itr) {
      // get current entry
      std::string label             = itr->first;
      Teuchos::ParameterList & list = itr->second.getValue(temp);
      list.set("Type",label);
      
      // add to library
      invLib->addInverse(label,list);
   }
 
   return invLib;
}

Teuchos::RCP<InverseLibrary> InverseLibrary::buildFromStratimikos(const Stratimikos::DefaultLinearSolverBuilder & strat)
{
   RCP<InverseLibrary> invLib = rcp(new InverseLibrary());
 
   // get default inveres in Stratimikos
   RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList(*strat.getValidParameters()));
   Teuchos::ParameterList lst(pl->sublist("Linear Solver Types"));
   Teuchos::ParameterList pft(pl->sublist("Preconditioner Types"));
 
   Teuchos::ParameterList::ConstIterator itr;
   Teuchos::ParameterList * temp = 0;
 
   // loop over all entries in solver list
   for(itr=lst.begin();itr!=lst.end();++itr) {
      // get current entry
      std::string label             = itr->first;
      Teuchos::ParameterList & list = itr->second.getValue(temp);
      list.set("Type",label);
      
      // add to library
      invLib->addInverse(label,list);
   }
 
   // loop over all entries in preconditioner list
   for(itr=pft.begin();itr!=pft.end();++itr) {
      // get current entry
      std::string label             = itr->first;
      Teuchos::ParameterList & list = itr->second.getValue(temp);
      list.set("Type",label);
      
      // add to library
      invLib->addInverse(label,list);
   }
 
   return invLib;
}

Teuchos::RCP<InverseLibrary> InverseLibrary::buildFromStratimikos(const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> & strat)
{
   RCP<InverseLibrary> invLib = rcp(new InverseLibrary(strat));
 
   // get default inveres in Stratimikos
   RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList(*strat->getValidParameters()));
   Teuchos::ParameterList lst(pl->sublist("Linear Solver Types"));
   Teuchos::ParameterList pft(pl->sublist("Preconditioner Types"));
 
   Teuchos::ParameterList::ConstIterator itr;
   Teuchos::ParameterList * temp = 0;
 
   // loop over all entries in solver list
   for(itr=lst.begin();itr!=lst.end();++itr) {
      // get current entry
      std::string label             = itr->first;
      Teuchos::ParameterList & list = itr->second.getValue(temp);
      list.set("Type",label);
      
      // add to library
      invLib->addInverse(label,list);
   }
 
   // loop over all entries in preconditioner list
   for(itr=pft.begin();itr!=pft.end();++itr) {
      // get current entry
      std::string label             = itr->first;
      Teuchos::ParameterList & list = itr->second.getValue(temp);
      list.set("Type",label);
      
      // add to library
      invLib->addInverse(label,list);
   }
 
   return invLib;
}
 
} // end namespace Teko