Amesos Package Browser (Single Doxygen Collection) Development
simpleStratimikosSolve.cpp
Go to the documentation of this file.
00001 #include "simpleStratimikosSolve.hpp"
00002 
00003 #include "Epetra_CrsMatrix.h"
00004 #include "Epetra_MultiVector.h"
00005 #include "Teuchos_ParameterList.hpp"
00006 #include "Teuchos_RCP.hpp"
00007 #include "Thyra_EpetraLinearOp.hpp"
00008 #include "Thyra_SpmdVectorSpaceBase.hpp"
00009 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00010 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00011 #include "Thyra_LinearOpWithSolveFactoryBase.hpp"
00012 
00013 #include "Thyra_EpetraThyraWrappers.hpp"  // Contains create_MultiVector
00014 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"  // Contains LinearOpWithSolve ?
00015 
00016 /*
00017   simpleStratimikosSolve performs x = A \ b; 
00018   and returns 0 on success
00019  */
00020 
00021 int simpleStratimikosSolve( 
00022    Epetra_CrsMatrix const& epetra_A,  // non-persisting, non-changeable
00023    Epetra_MultiVector const& epetra_B,  // non-persisting, non-changeable
00024    Epetra_MultiVector *epetra_X,  // non-persisting, changeable
00025    Teuchos::ParameterList *paramList  // non-persisting, changeable
00026    )
00027   {
00028 
00029     using Teuchos::RCP;
00030     using Teuchos::rcp;
00031 
00032 
00033     //
00034     // C) The "Glue" code that takes Epetra objects and wraps them as Thyra
00035     // objects
00036     //
00037     // This next set of code wraps the Epetra objects that define the linear
00038     // system to be solved as Thyra objects so that they can be passed to the
00039     // linear solver.
00040     //
00041  
00042     // Create RCPs that will be used to hold the Thyra wrappers
00043 
00044     typedef RCP<const Thyra::LinearOpBase<double> > LinearOpPtr;
00045     typedef RCP<Thyra::MutliVectorBase<double> > MultiVectorPtr;
00046 
00047     LinearOpPtr A = Thyra::epetraLinearOp( epetra_A );
00048     VectorPtr   X = Thyra::create_Vector( rcp(epetra_X,false), A->domain() );
00049     VectorPtr   B = Thyra::create_Vector( rcp(&epetra_B,false), A->range() );
00050  
00051     // Note that above Thyra is only interacted with in the most trival of
00052     // ways.  For most users, Thyra will only be seen as a thin wrapper that
00053     // they need know little about in order to wrap their objects in order to
00054     // pass them to Thyra-enabled solvers.
00055 
00056      
00057     //
00058     // D) Thyra-specific code for solving the linear system
00059     //
00060     // Note that this code has no mention of any concrete implementation and
00061     // therefore can be used in any use case.
00062     //
00063  
00064     Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00065     // Set the parameter list
00066     linearSolverBuilder.setParameterList( rcp(paramList,false) ) ; 
00067     // Create a linear solver factory given information read from the
00068     // parameter list.
00069     RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
00070       lowsFactory = linearSolverBuilder.createLinearSolveStrategy("");
00071     // Setup the verbosity level
00072     lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
00073     // Create a linear solver based on the forward operator A
00074     RCP<Thyra::LinearOpWithSolveBase<double> >
00075       lows = Thyra::linearOpWithSolve(*lowsFactory,A);
00076       //      lows = Thyra::linearOpWithSolve(*lowsFactory,rcp(&A,false));
00077     // Solve the linear system (note: the initial guess in 'X' is critical)
00078     Thyra::SolveStatus<double>
00079       status = Thyra::solve(*lows,Thyra::NOTRANS,*B,&*X);
00080     // Write the linear solver parameters after they were read
00081     linearSolverBuilder.writeParamsFile(*lowsFactory);
00082 
00083 
00084 #if 0
00085     std::cout << __FILE__ << "::" << __LINE__ << 
00086       " paramlist  = " << *(linearSolverBuilder.getParameterList( ))   << std::endl ; 
00087 #endif
00088 
00089     return (status.solveStatus!=Thyra::SOLVE_STATUS_CONVERGED);
00090  
00091   }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines