Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
belos_tpetra_thyra_lowsf_lap.cpp
Go to the documentation of this file.
1 
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Stratimikos: Thyra-based strategies for linear solvers
6 // Copyright (2006) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 //
43 // This driver creates a 2D Laplacian operator as a Tpetra::CisMatrix.
44 // This matrix is then converted into a Thyra linear operator
45 // through the Thyra-Tpetra adapters.
46 //
47 // The right-hand sides are all random vectors.
48 // The initial guesses are all set to zero.
49 //
50 //
51 
52 // Thyra includes
53 #include "Thyra_BelosLinearOpWithSolveFactory.hpp"
54 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
55 #include "Thyra_TpetraLinearOp.hpp"
56 
57 // Tpetra includes
58 #include "Tpetra_ElementSpace.hpp"
59 #include "Tpetra_VectorSpace.hpp"
60 #include "Tpetra_CisMatrix.hpp"
61 
62 // Teuchos includes
63 #include "Teuchos_ParameterList.hpp"
64 #include "Teuchos_VerboseObject.hpp"
65 #include "Teuchos_GlobalMPISession.hpp"
66 
67 #ifdef HAVE_MPI
68 # include "Tpetra_MpiPlatform.hpp"
69 #else
70 # include "Tpetra_SerialPlatform.hpp"
71 #endif
72 
73 int main(int argc, char* argv[])
74 {
75  //
76  // Get a default output stream from the Teuchos::VerboseObjectBase
77  //
78  Teuchos::RCP<Teuchos::FancyOStream>
79  out = Teuchos::VerboseObjectBase::getDefaultOStream();
80 
81  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
82 
83  // --------------------------------------------------------------------------------
84  //
85  // Create the scalar-type typedefs
86  //
87  // --------------------------------------------------------------------------------
88 
89 #ifdef HAVE_COMPLEX
90  typedef std::complex<double> ST;
91 #elif HAVE_COMPLEX_H
92  typedef std::complex<double> ST;
93 #else
94  typedef double ST;
95 #endif
96  typedef Teuchos::ScalarTraits<ST>::magnitudeType MT;
97  typedef int OT;
98 
99  // --------------------------------------------------------------------------------
100  //
101  // Create 2D Laplacian using Tpetra::CisMatrix
102  //
103  // --------------------------------------------------------------------------------
104 
105 #ifdef HAVE_MPI
106  MPI_Comm mpiComm = MPI_COMM_WORLD;
107  const Tpetra::MpiPlatform<OT,OT> ordinalPlatform(mpiComm);
108  const Tpetra::MpiPlatform<OT,ST> scalarPlatform(mpiComm);
109 #else
110  const Tpetra::SerialPlatform<OT,OT> ordinalPlatform;
111  const Tpetra::SerialPlatform<OT,ST> scalarPlatform;
112 #endif
113 
114  // Declare global dimension of the linear operator
115  OT globalDim = 500;
116 
117  // Create the element space and std::vector space
118  const Tpetra::ElementSpace<OT> elementSpace(globalDim,0,ordinalPlatform);
119  const Tpetra::VectorSpace<OT,ST> vectorSpace(elementSpace,scalarPlatform);
120 
121  // Allocate the Tpetra::CisMatrix object.
122  RCP<Tpetra::CisMatrix<OT,ST> >
123  tpetra_A = rcp(new Tpetra::CisMatrix<OT,ST>(vectorSpace));
124 
125  // Get the indexes of the rows on this processor
126  const int numMyElements = vectorSpace.getNumMyEntries();
127  const std::vector<int> &myGlobalElements = vectorSpace.elementSpace().getMyGlobalElements();
128 
129  // Fill the local matrix entries one row at a time.
130  const ST offDiag = -1.0, diag = 2.0;
131  int numEntries; ST values[3]; int indexes[3];
132  for( int k = 0; k < numMyElements; ++k ) {
133  const int rowIndex = myGlobalElements[k];
134  if( rowIndex == 0 ) { // First row
135  numEntries = 2;
136  values[0] = diag; values[1] = offDiag;
137  indexes[0] = 0; indexes[1] = 1;
138  }
139  else if( rowIndex == globalDim - 1 ) { // Last row
140  numEntries = 2;
141  values[0] = offDiag; values[1] = diag;
142  indexes[0] = globalDim-2; indexes[1] = globalDim-1;
143  }
144  else { // Middle rows
145  numEntries = 3;
146  values[0] = offDiag; values[1] = diag; values[2] = offDiag;
147  indexes[0] = rowIndex-1; indexes[1] = rowIndex; indexes[2] = rowIndex+1;
148  }
149  tpetra_A->submitEntries(Tpetra::Insert,rowIndex,numEntries,values,indexes);
150  }
151 
152  // Finish the construction of the Tpetra::CisMatrix
153  tpetra_A->fillComplete();
154 
155  // Create a Thyra linear operator (A) using the Tpetra::CisMatrix (tpetra_A).
156  RCP<Thyra::LinearOpBase<ST> >
157  A = Teuchos::rcp( new Thyra::TpetraLinearOp<OT,ST>(
158  Teuchos::rcp_implicit_cast<Tpetra::Operator<OT,ST> >(tpetra_A) ) );
159 
160  // Set the parameters for the Belos LOWS Factory and create a parameter list.
161  // NOTE: All the code below only uses Thyra and is independent of Tpetra.
162  //
163  int blockSize = 1; // can only be 1 right now.
164  int maxIterations = globalDim;
165  int maxRestarts = 0;
166  int gmresKrylovLength = globalDim;
167  int outputFrequency = 100;
168  bool outputMaxResOnly = true;
169  MT maxResid = 1e-3;
170 
171  ST one = Teuchos::ScalarTraits<ST>::one();
172  ST zero = Teuchos::ScalarTraits<ST>::zero();
173 
174  Teuchos::RCP<Teuchos::ParameterList>
175  belosLOWSFPL = Teuchos::rcp( new Teuchos::ParameterList() );
176 
177  belosLOWSFPL->set("Solver Type","Block GMRES");
178 
179  Teuchos::ParameterList& belosLOWSFPL_solver =
180  belosLOWSFPL->sublist("Solver Types");
181 
182  Teuchos::ParameterList& belosLOWSFPL_gmres =
183  belosLOWSFPL_solver.sublist("Block GMRES");
184 
185  belosLOWSFPL_gmres.set("Maximum Iterations",int(maxIterations));
186  belosLOWSFPL_gmres.set("Convergence Tolerance",double(maxResid));
187  belosLOWSFPL_gmres.set("Maximum Restarts",int(maxRestarts));
188  belosLOWSFPL_gmres.set("Block Size",int(blockSize));
189  belosLOWSFPL_gmres.set("Num Blocks",int(gmresKrylovLength));
190  belosLOWSFPL_gmres.set("Output Frequency",int(outputFrequency));
191  belosLOWSFPL_gmres.set("Show Maximum Residual Norm Only",bool(outputMaxResOnly));
192 
193  // Whether the linear solver succeeded.
194  // (this will be set during the residual check at the end)
195  bool success = true;
196 
197  // --------------------------------------------------------------------------------
198  //
199  // Create the right-hand sides and solution vectors
200  //
201  // --------------------------------------------------------------------------------
202 
203  // Number of random right-hand sides we will be solving for.
204  int numRhs = 5;
205 
206  // Create smart pointer to right-hand side and solution std::vector to be filled in below.
207  Teuchos::RCP<Thyra::MultiVectorBase<ST> > x, b;
208 
209  if (numRhs==1) {
210  //
211  // In this case we can construct vectors using Tpetra and just "wrap" them in Thyra objects.
212  //
213 
214  // Create RHS std::vector
215  Teuchos::RCP<Tpetra::Vector<OT,ST> > tpetra_b =
216  Teuchos::rcp( new Tpetra::Vector<OT,ST>(vectorSpace) );
217 
218  // Randomize RHS std::vector
219  tpetra_b->setAllToRandom();
220 
221  // Wrap Tpetra std::vector as Thyra std::vector
222  b = Thyra::create_Vector(tpetra_b);
223 
224  // Create solution (LHS) std::vector
225  Teuchos::RCP<Tpetra::Vector<OT,ST> > tpetra_x =
226  Teuchos::rcp( new Tpetra::Vector<OT,ST>(vectorSpace) );
227 
228  // Initialize solution to zero
229  tpetra_x->setAllToScalar( zero );
230 
231  // Wrap Tpetra std::vector as Thyra std::vector
232  x = Thyra::create_Vector(tpetra_x);
233 
234  }
235  else {
236  //
237  // In this case we can construct the multivector using Thyra and extract columns of
238  // the multivector as Tpetra std::vector, which can then be filled. This is because
239  // Tpetra does not have a multivector object and Thyra will emulate the multivector.
240  //
241 
242  // Get the domain space for the Thyra linear operator
243  Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domain = A->domain();
244 
245  // Create a right-hand side and solution vectors with numRhs vectors in it.
246  x = Thyra::createMembers(domain, numRhs);
247  b = Thyra::createMembers(domain, numRhs);
248 
249  // Extract the Tpetra std::vector from the columns of the multivector and fill them with
250  // random numbers (b) or zeros (x).
251 
252  for ( int j=0; j<numRhs; ++j ) {
253  //
254  // Get the j-th column from b as a Tpetra std::vector and randomize it.
255  //
256  Teuchos::RCP<Tpetra::Vector<OT,ST> >
257  tpetra_b_j = Thyra::get_Tpetra_Vector(vectorSpace,b->col(j));
258  tpetra_b_j->setAllToRandom();
259  //
260  // Get the j-th column from x as a Tpetra std::vector and set it to zero.
261  //
262  Teuchos::RCP<Tpetra::Vector<OT,ST> >
263  tpetra_x_j = Thyra::get_Tpetra_Vector(vectorSpace,x->col(j));
264  tpetra_x_j->setAllToScalar( zero );
265  //
266  // NOTE: Tpetra vectors have element access via the [] operator.
267  // So and additional approach for filling in a Tpetra std::vector is:
268  //
269  for ( int i=0; i<numMyElements; ++i ) {
270  //
271  // Get the global index.
272  //
273  const int rowIndex = myGlobalElements[ i ];
274 
275  // Set the entry to zero.
276  (*tpetra_x_j)[ rowIndex ] = zero;
277  }
278  }
279  }
280 
281  // --------------------------------------------------------------------------------
282  //
283  // Create the linear operator with solve factory/object and solve the linear
284  // system using the iterative solvers in Belos.
285  //
286  // NOTE: This is the only part of the code that solely uses Thyra.
287  //
288  // --------------------------------------------------------------------------------
289 
290  // Create the Belos LOWS factory.
291  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ST> >
292  belosLOWSFactory = Teuchos::rcp(new Thyra::BelosLinearOpWithSolveFactory<ST>());
293 
294  // Set the output stream and the verbosity level (prints to std::cout by defualt)
295  // NOTE: Set to VERB_NONE for no output from the solver.
296  belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
297 
298  // Set the parameter list to specify the behavior of the factory.
299  belosLOWSFactory->setParameterList( belosLOWSFPL );
300 
301  // Create a BelosLinearOpWithSolve object from the Belos LOWS factory.
302  Teuchos::RCP<Thyra::LinearOpWithSolveBase<ST> >
303  nsA = belosLOWSFactory->createOp();
304 
305  // Initialize the BelosLinearOpWithSolve object with the Thyra linear operator.
306  Thyra::initializeOp<ST>( *belosLOWSFactory, A, &*nsA );
307 
308  // Perform solve using the linear operator to get the approximate solution of Ax=b,
309  // where b is the right-hand side and x is the left-hand side.
310  Thyra::SolveStatus<ST> solveStatus;
311  solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
312 
313  // Print out status of solve.
314  *out << "\nBelos LOWS Status: "<< solveStatus << std::endl;
315 
316  // --------------------------------------------------------------------------------
317  // Compute residual and check convergence.
318  // NOTE: We will do this by extracting the Tpetra vectors from the multivector.
319  //
320  // NOTE 2: We are using scalar traits here instead of the magnitude type
321  // because Tpetra defines its norm methods to return the scalar type.
322  // --------------------------------------------------------------------------------
323  std::vector<ST> norm_b(numRhs), norm_res(numRhs);
324 
325  for ( int j=0; j<numRhs; ++j ) {
326  //
327  // Get the j-th columns from x and b as Tpetra vectors.
328  //
329  Teuchos::RCP<Tpetra::Vector<OT,ST> >
330  tpetra_x_j = Thyra::get_Tpetra_Vector(vectorSpace,x->col(j));
331 
332  Teuchos::RCP<Tpetra::Vector<OT,ST> >
333  tpetra_b_j = Thyra::get_Tpetra_Vector(vectorSpace,b->col(j));
334 
335  // Compute the column norms of the right-hand side b.
336  norm_b[j] = tpetra_b_j->norm2();
337 
338  // Compute y=A*x, where x is the solution from the linear solver.
339  Tpetra::Vector<OT,ST> y(vectorSpace);
340  tpetra_A->apply( *tpetra_x_j, y );
341 
342  // Compute b-A*x = b-y
343  y.update( one, *tpetra_b_j, -one );
344 
345  // Compute the column norms of A*x-b.
346  norm_res[j] = y.norm2();
347  }
348 
349  // Print out the final relative residual norms.
350  MT rel_res = 0.0;
351  *out << "Final relative residual norms" << std::endl;
352  for (int i=0; i<numRhs; ++i) {
353  rel_res = Teuchos::ScalarTraits<ST>::real(norm_res[i]/norm_b[i]);
354  if (rel_res > maxResid)
355  success = false;
356  *out << "RHS " << i+1 << " : "
357  << std::setw(16) << std::right << rel_res << std::endl;
358  }
359 
360  return ( success ? 0 : 1 );
361 }
double ST
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.
int main(int argc, char *argv[])