Amesos Package Browser (Single Doxygen Collection) Development
stratimikos_example.cpp
Go to the documentation of this file.
00001 //
00002 //
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //         Stratimikos: Thyra-based strategies for linear solvers
00007 //                Copyright (2006) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00027 // 
00028 // ***********************************************************************
00029 // @HEADER
00030 
00031 #include "Stratimikos_Config.h"
00032 #include "Teuchos_ConfigDefs.hpp"
00033 #include "Teuchos_FancyOStream.hpp"
00034 #include "simpleStratimikosSolve.hpp"
00035 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00036 #include "Thyra_EpetraLinearOp.hpp"
00037 //  #include "Thyra_LinearOpWithSolveFactoryExamples.hpp"
00038 #include "EpetraExt_readEpetraLinearSystem.h"
00039 #include "Teuchos_ParameterList.hpp"
00040 
00041 namespace Teuchos { class ParameterList; }
00042 
00043 #ifdef HAVE_MPI
00044 #  include "Epetra_MpiComm.h"
00045 #else
00046 #  include "Epetra_SerialComm.h"
00047 #endif
00048 
00049 bool TestSingleStratimikosSolver(
00050   Teuchos::ParameterList                  *paramList_inout
00051   ,const bool                             dumpAll
00052   ,Teuchos::FancyOStream                  *out
00053   )
00054 {
00055 
00056   using Teuchos::rcp;
00057   using Teuchos::RCP;
00058   using Teuchos::OSTab;
00059   using Teuchos::ParameterList;
00060   using Teuchos::getParameter;
00061   bool result, success = true;
00062 
00063   try {
00064 
00065     TEUCHOS_TEST_FOR_EXCEPT(!paramList_inout);
00066 
00067     RCP<ParameterList>
00068       paramList = rcp(paramList_inout,false);
00069 
00070 #if 0
00071     if(out) {
00072       *out << "\nEchoing input parameters ...\n";
00073       paramList->print(*out,1,true,false);
00074     }
00075 #endif
00076 
00077     // Create list of valid parameter sublists
00078     Teuchos::ParameterList validParamList("TestSingleStratimikosSolver");
00079     validParamList.set("Matrix File","fileName");
00080     validParamList.sublist("Linear Solver Builder");
00081     
00082     Teuchos::ParameterList &StratimikosList = paramList_inout->sublist("Linear Solver Builder");
00083 
00084 #if 0 
00085 
00086     this fails because validParamList has only been set one level deep - see 5-10 lines above
00087 
00088     if(out) *out << "\nValidating top-level input parameters ...\n";
00089     StratimikosList.validateParameters(validParamList.sublist("Linear Solver Builder"),0);
00090 #endif
00091 
00092 #if 0
00093     paramList->validateParameters(validParamList,0);
00094 
00095     if(out) *out << "\nValidating top-level input parameters ...\n";
00096     paramList->validateParameters(StratimikosList,0);
00097 #endif
00098 
00099     const std::string
00100       &matrixFile = getParameter<std::string>(*paramList,"Matrix File");
00101 
00102     if(out) *out << "\nReading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
00103   
00104 #ifdef HAVE_MPI
00105     Epetra_MpiComm comm(MPI_COMM_WORLD);
00106 #else
00107     Epetra_SerialComm comm;
00108 #endif
00109     Teuchos::RCP<Epetra_CrsMatrix> epetra_A;
00110     EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
00111 
00112     const int num_vecs = 1 ; 
00113     const Epetra_Map& A_domainmap = epetra_A->DomainMap() ; 
00114     Epetra_MultiVector X( A_domainmap, num_vecs ) ; 
00115     const Epetra_Map& A_rangemap = epetra_A->RangeMap() ; 
00116     Epetra_MultiVector B( A_rangemap, num_vecs ) ; 
00117     B.Random(); 
00118 
00119     int status = simpleStratimikosSolve( *epetra_A, B, &X, &StratimikosList );
00120     assert( status==0 ) ; 
00121 
00122 
00123     int nrows =epetra_A->NumGlobalRows();
00124     Epetra_MultiVector Residual( B ) ;
00125     epetra_A->Apply( X, Residual );
00126     Residual.Update( 1.0, B, -1.0 );
00127     double ResidualNorm;
00128     double BNorm;
00129     Residual.Norm2( &ResidualNorm );
00130     B.Norm2( &BNorm );
00131 
00132     const Teuchos::ParameterList& LST_PL =  StratimikosList.sublist("Linear Solver Types"); 
00133     const Teuchos::ParameterList& Amesos_PL =  LST_PL.sublist("Amesos"); 
00134     std::string SolverType = Amesos_PL.get<std::string>("Solver Type");
00135     assert( SolverType == "Lapack" ) ; // for now - remove later 
00136     const Teuchos::ParameterList& AmesosSettings_PL =  Amesos_PL.sublist("Amesos Settings"); 
00137     const Teuchos::ParameterList& SolverType_PL =  AmesosSettings_PL.sublist(SolverType); 
00138     double AddToDiag = -13.0;
00139     if (  SolverType == "Lapack" ) {
00140       AddToDiag = SolverType_PL.get<double>("AddToDiag");
00141     }
00142     assert( AddToDiag >= 0.0 ) ; 
00143     double MaxError = 1e-15*BNorm + 10.0 * AddToDiag ; 
00144     double MinError = 0.02 * ( 1e-15*BNorm + AddToDiag ); 
00145     success = ResidualNorm < nrows * MaxError && 
00146       ResidualNorm > MinError ; 
00147 
00148 #if 0
00149     std::cout << __FILE__ << "::" << __LINE__ 
00150    << " ResidualNorm = " << ResidualNorm  
00151    << " MinError = " << MinError 
00152    << " MaxError = " << MaxError << std::endl ; 
00153     std::cout << " B = " ; 
00154     B.Print( std::cout ) ; 
00155     std::cout << " X = " ; 
00156     X.Print( std::cout ) ; 
00157     std::cout << " epetra_A = " ; 
00158     epetra_A->Print( std::cout ) ; 
00159     std::cout << " success = " << success << " ResidualNorm = " <<  ResidualNorm << std::endl ; 
00160 #endif
00161 
00162     if(false && out) {
00163       *out << "\nPrinting the parameter list (showing what was used) ...\n";
00164       paramList->print(*out,1,true,true);
00165     }
00166     
00167   }
00168   catch( const std::exception &excpt ) {
00169     std::cerr << "*** Caught standard exception : " << excpt.what() << std::endl;
00170     success = false;
00171   }
00172   
00173   return success;
00174   
00175 }
00176 
00177 #include "Teuchos_GlobalMPISession.hpp"
00178 #include "Teuchos_VerboseObject.hpp"
00179 #include "Teuchos_CommandLineProcessor.hpp"
00180 #include "Teuchos_XMLParameterListHelpers.hpp"
00181 #include "Teuchos_StandardCatchMacros.hpp"
00182 
00183 int main(int argc, char* argv[])
00184 {
00185 
00186   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00187   
00188   using Teuchos::CommandLineProcessor;
00189 
00190   bool success = true;
00191   bool verbose = true;
00192 
00193   Teuchos::RCP<Teuchos::FancyOStream>
00194     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00195 
00196   try {
00197 
00198     //
00199     // Read options from command-line
00200     //
00201     
00202     std::string     inputFile              = "stratimikos_amesos_lapack.xml";
00203     std::string     extraParams            = "";
00204     bool            dumpAll                = false;
00205 
00206     CommandLineProcessor  clp(false); // Don't throw exceptions
00207     std::string default_inputFile ;
00208     sprintf( &default_inputFile[0], "Input file [%s]", &inputFile[0] );
00209     clp.setOption( "input-file", &inputFile, &default_inputFile[0], false );
00210     clp.setOption( "extra-params", &extraParams, "Extra parameters overriding the parameters read in from --input-file");
00211     clp.setDocString(
00212       "Testing program for Trilinos (and non-Trilinos) linear solvers access through Thyra."
00213       );
00214 
00215     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00216     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00217 
00218     Teuchos::ParameterList paramList;
00219     if(verbose) *out << "\nReading parameters from XML file \""<<inputFile<<"\" ...\n";
00220     Teuchos::updateParametersFromXmlFile(inputFile,&paramList);
00221     if(extraParams.length()) {
00222       if(verbose) *out << "\nAppending extra parameters from the XML string \""<<extraParams<<"\" ...\n";
00223       Teuchos::updateParametersFromXmlString(extraParams,&paramList);
00224     }
00225     
00226     success
00227       = TestSingleStratimikosSolver(
00228         &paramList,dumpAll,verbose?&*out:0
00229         );
00230     
00231   }
00232   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success)
00233   
00234   if (verbose) {
00235     if(success)  *out << "\nCongratulations! All of the tests checked out!\n";
00236     else         *out << "\nOh no! At least one of the tests failed!\n";
00237   }
00238 
00239   return ( success ? 0 : 1 );
00240 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines