Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
test_thyra_tpetra_ifpack2_issue_535.cpp
Go to the documentation of this file.
2#include "Teuchos_DefaultComm.hpp"
3#include "Tpetra_CrsMatrix.hpp"
4#include "Thyra_Ifpack2PreconditionerFactory.hpp"
5#include "Thyra_TpetraThyraWrappers.hpp"
6#include <cstdlib>
7
8using Teuchos::RCP;
9using Teuchos::rcp;
10typedef Tpetra::Map<int, int> map_type;
11typedef map_type::local_ordinal_type LO;
12typedef map_type::global_ordinal_type GO;
13typedef double ST;
14typedef Tpetra::CrsMatrix<ST, LO, GO> crs_matrix_type;
15typedef Tpetra::Operator<ST, LO, GO> op_type;
16typedef Tpetra::Vector<ST, LO, GO> vec_type;
17
18Teuchos::RCP<const crs_matrix_type>
19create_matrix (const Teuchos::RCP<const Teuchos::Comm<int> >& comm)
20{
21 const Tpetra::global_size_t numGlobalElements = 100;
22
23 const GO indexBase = 0;
24 auto map = rcp (new map_type (numGlobalElements, indexBase, comm));
25
26 const LO numMyElements = map->getLocalNumElements ();
27 auto myGlobalElements = map->getLocalElementList ();
28 auto A = rcp (new crs_matrix_type (map, 1));
29
30 for (LO lclRow = 0; lclRow < numMyElements; ++lclRow) {
31 const GO gblInd = map->getGlobalElement (lclRow);
32 const ST val = static_cast<ST> (1.0) / static_cast<ST> (gblInd + 1);
33 A->insertGlobalValues (gblInd,
34 Teuchos::tuple (gblInd),
35 Teuchos::tuple (val));
36 }
37 A->fillComplete ();
38 return A;
39}
40
41int
42main (int argc, char *argv[])
43{
44 using Teuchos::outArg;
45 using Teuchos::REDUCE_MIN;
46 using Teuchos::reduceAll;
47 using std::cerr;
48 using std::endl;
49 int lclSuccess = 1;
50 int gblSuccess = 0;
51
52 Teuchos::GlobalMPISession session (&argc, &argv, NULL);
53 auto out = Teuchos::VerboseObjectBase::getDefaultOStream ();
54 const auto comm = Teuchos::DefaultComm<int>::getComm ();
55 const int myRank = comm->getRank ();
56
57 *out << "Creating matrix" << endl;
58 RCP<const crs_matrix_type> A;
59 try {
60 A = create_matrix (comm);
61 }
62 catch (std::exception& e) {
63 lclSuccess = 0;
64 std::ostringstream os;
65 os << "Proc " << myRank << ": create_matrix(comm) threw an "
66 "exception: " << e.what () << endl;
67 cerr << os.str ();
68 }
69 catch (...) {
70 lclSuccess = 0;
71 std::ostringstream os;
72 os << "Proc " << myRank << ": create_matrix(comm) threw an "
73 "exception not a subclass of std::exception." << endl;
74 cerr << os.str ();
75 }
76 gblSuccess = 0;
77 reduceAll<int, int> (*comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
78 if (gblSuccess != 1) {
79 if (myRank == 0) {
80 *out << "create_matrix(comm) threw an exception on some process."
81 << endl;
82 }
83 *out << "End Result: TEST FAILED" << endl;
84 return EXIT_FAILURE;
85 }
86
87 // Make sure that A is nonnull on all processes.
88 if (A.is_null ()) {
89 lclSuccess = 0;
90 }
91 gblSuccess = 0;
92 reduceAll<int, int> (*comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
93 if (gblSuccess != 1) {
94 if (myRank == 0) {
95 *out << "The result of create_matrix(comm) is null on at least one "
96 "process." << endl;
97 }
98 *out << "End Result: TEST FAILED" << endl;
99 return EXIT_FAILURE;
100 }
101
102 *out << "Creating vectors" << endl;
103 vec_type b (A->getRangeMap ());
104 b.putScalar (1.0);
105 vec_type x (A->getDomainMap ());
106 x.putScalar (0.0);
107
108 *out << "Creating Stratimikos linear solver builder" << endl;
110 auto p = Teuchos::rcp (new Teuchos::ParameterList ());
111 try {
112 builder.setParameterList (p);
113 }
114 catch (std::exception& e) {
115 lclSuccess = 0;
116 std::ostringstream os;
117 os << "Proc " << myRank << ": builder.setParameterList(p) threw an "
118 "exception: " << e.what () << endl;
119 cerr << os.str ();
120 }
121 catch (...) {
122 lclSuccess = 0;
123 std::ostringstream os;
124 os << "Proc " << myRank << ": builder.setParameterList(p) threw an "
125 "exception not a subclass of std::exception." << endl;
126 cerr << os.str ();
127 }
128 gblSuccess = 0;
129 reduceAll<int, int> (*comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
130 if (gblSuccess != 1) {
131 if (myRank == 0) {
132 *out << "builder.setParameterList(p) threw an exception on some process."
133 << endl;
134 }
135 *out << "End Result: TEST FAILED" << endl;
136 return EXIT_FAILURE;
137 }
138
139 *out << "Calling builder.createLinearSolveStrategy" << endl;
140 auto lowsFactory = builder.createLinearSolveStrategy ("");
141 lowsFactory->setVerbLevel (Teuchos::VERB_LOW);
142
143 *out << "Calling Thyra::createConstLinearOp" << endl;
144 const op_type& opA = *A;
145 auto thyraA = Thyra::createConstLinearOp (Teuchos::rcpFromRef (opA));
146
147 // using Teuchos::rcp_implicit_cast;
148 // auto thyraA = Thyra::createConstLinearOp (rcp_implicit_cast<const op_type> (A));
149
150 // Make sure that thyraA is nonnull on all processes.
151 if (thyraA.is_null ()) {
152 lclSuccess = 0;
153 }
154 gblSuccess = 0;
155 reduceAll<int, int> (*comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
156 if (gblSuccess != 1) {
157 if (myRank == 0) {
158 *out << "The result of Thyra::createConstLinearOp is null on at least "
159 "one process." << endl;
160 }
161 *out << "End Result: TEST FAILED" << endl;
162 return EXIT_FAILURE;
163 }
164
165 *out << "Creating Thyra Ifpack2 factory" << endl;
166 RCP<Thyra::PreconditionerFactoryBase<double> > factory =
167 rcp (new Thyra::Ifpack2PreconditionerFactory<crs_matrix_type> ());
168
169 *out << "Creating Ifpack2 preconditioner using factory" << endl;
170 typedef Thyra::PreconditionerBase<ST> prec_type;
171 RCP<prec_type> prec;
172 try {
173 prec = factory->createPrec ();
174 }
175 catch (std::exception& e) {
176 lclSuccess = 0;
177 std::ostringstream os;
178 os << "Proc " << myRank << ": factory->createPrec() threw an "
179 "exception: " << e.what () << endl;
180 cerr << os.str ();
181 }
182 catch (...) {
183 lclSuccess = 0;
184 std::ostringstream os;
185 os << "Proc " << myRank << ": factory->createPrec() threw an "
186 "exception not a subclass of std::exception." << endl;
187 cerr << os.str ();
188 }
189 gblSuccess = 0;
190 reduceAll<int, int> (*comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
191 if (gblSuccess != 1) {
192 if (myRank == 0) {
193 *out << "factory->createPrec() threw an exception on some process."
194 << endl;
195 }
196 *out << "End Result: TEST FAILED" << endl;
197 return EXIT_FAILURE;
198 }
199
200 // Make sure that prec is nonnull on all processes.
201 if (prec.is_null ()) {
202 lclSuccess = 0;
203 }
204 gblSuccess = 0;
205 reduceAll<int, int> (*comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
206 if (gblSuccess != 1) {
207 if (myRank == 0) {
208 *out << "The result of factory->createPrec() is null on at least one "
209 "process." << endl;
210 }
211 *out << "End Result: TEST FAILED" << endl;
212 return EXIT_FAILURE;
213 }
214
215 Thyra::initializePrec (*factory, thyraA, prec.ptr ()); // segfault!
216
217 // If this test starts failing, please check ifpack2/adapters/thyra/Thyra_Ifpack2PreconditionerFactory_def.hpp
218 // and corresponding spot in stratimikos/adapters/belos/tpetra/Thyra_BelosTpetraPreconditionerFactory_def.hpp
219 // See PR #11222 for most recent adjustments to these files. -JAL
220
221 *out << "End Result: TEST PASSED" << endl;
222 return EXIT_SUCCESS;
223}
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...
RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > createLinearSolveStrategy(const std::string &linearSolveStrategyName) const
void setParameterList(RCP< ParameterList > const &paramList)
Tpetra::Map< int, int > map_type
int main(int argc, char *argv[])
map_type::global_ordinal_type GO
Tpetra::Vector< ST, LO, GO > vec_type
Teuchos::RCP< const crs_matrix_type > create_matrix(const Teuchos::RCP< const Teuchos::Comm< int > > &comm)
Tpetra::CrsMatrix< ST, LO, GO > crs_matrix_type
map_type::local_ordinal_type LO
Tpetra::Operator< ST, LO, GO > op_type