Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
simpleStratimikosSolve.cpp
Go to the documentation of this file.
2
3#include "Epetra_CrsMatrix.h"
4#include "Epetra_MultiVector.h"
5#include "Teuchos_ParameterList.hpp"
6#include "Teuchos_RCP.hpp"
7#include "Thyra_EpetraLinearOp.hpp"
8#include "Thyra_SpmdVectorSpaceBase.hpp"
9#include "Thyra_DefaultSpmdVectorSpace.hpp"
10#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
11#include "Thyra_LinearOpWithSolveFactoryBase.hpp"
12
13#include "Thyra_EpetraThyraWrappers.hpp" // Contains create_MultiVector
14#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp" // Contains LinearOpWithSolve ?
15
16/*
17 simpleStratimikosSolve performs x = A \ b;
18 and returns 0 on success
19 */
20
22 Epetra_CrsMatrix const& epetra_A, // non-persisting, non-changeable
23 Epetra_MultiVector const& epetra_B, // non-persisting, non-changeable
24 Epetra_MultiVector *epetra_X, // non-persisting, changeable
25 Teuchos::ParameterList *paramList // non-persisting, changeable
26 )
27 {
28
29 using Teuchos::RCP;
30 using Teuchos::rcp;
31
32
33 //
34 // C) The "Glue" code that takes Epetra objects and wraps them as Thyra
35 // objects
36 //
37 // This next set of code wraps the Epetra objects that define the linear
38 // system to be solved as Thyra objects so that they can be passed to the
39 // linear solver.
40 //
41
42 // Create RCPs that will be used to hold the Thyra wrappers
43
44 typedef RCP<const Thyra::LinearOpBase<double> > LinearOpPtr;
45 typedef RCP<Thyra::MutliVectorBase<double> > MultiVectorPtr;
46
47 LinearOpPtr A = Thyra::epetraLinearOp( epetra_A );
48 VectorPtr X = Thyra::create_Vector( rcp(epetra_X,false), A->domain() );
49 VectorPtr B = Thyra::create_Vector( rcp(&epetra_B,false), A->range() );
50
51 // Note that above Thyra is only interacted with in the most trival of
52 // ways. For most users, Thyra will only be seen as a thin wrapper that
53 // they need know little about in order to wrap their objects in order to
54 // pass them to Thyra-enabled solvers.
55
56
57 //
58 // D) Thyra-specific code for solving the linear system
59 //
60 // Note that this code has no mention of any concrete implementation and
61 // therefore can be used in any use case.
62 //
63
64 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
65 // Set the parameter list
66 linearSolverBuilder.setParameterList( rcp(paramList,false) ) ;
67 // Create a linear solver factory given information read from the
68 // parameter list.
69 RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
70 lowsFactory = linearSolverBuilder.createLinearSolveStrategy("");
71 // Setup the verbosity level
72 lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
73 // Create a linear solver based on the forward operator A
74 RCP<Thyra::LinearOpWithSolveBase<double> >
75 lows = Thyra::linearOpWithSolve(*lowsFactory,A);
76 // lows = Thyra::linearOpWithSolve(*lowsFactory,rcp(&A,false));
77 // Solve the linear system (note: the initial guess in 'X' is critical)
78 Thyra::SolveStatus<double>
79 status = Thyra::solve(*lows,Thyra::NOTRANS,*B,&*X);
80 // Write the linear solver parameters after they were read
81 linearSolverBuilder.writeParamsFile(*lowsFactory);
82
83
84#if 0
85 std::cout << __FILE__ << "::" << __LINE__ <<
86 " paramlist = " << *(linearSolverBuilder.getParameterList( )) << std::endl ;
87#endif
88
89 return (status.solveStatus!=Thyra::SOLVE_STATUS_CONVERGED);
90
91 }
int simpleStratimikosSolve(Epetra_CrsMatrix const &epetra_A, Epetra_MultiVector const &epetra_B, Epetra_MultiVector *epetra_X, Teuchos::ParameterList *paramList)