Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
belos_tpetra_thyra_lowsf_prec_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 // and a preconditioner. These matrices are then converted into Thyra
45 // linear operators 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_PreconditionerFactoryHelpers.hpp"
56 #include "Thyra_TpetraLinearOp.hpp"
57 #include "Thyra_DefaultPreconditioner.hpp"
58 
59 // Tpetra includes
60 #include "Tpetra_ElementSpace.hpp"
61 #include "Tpetra_VectorSpace.hpp"
62 #include "Tpetra_CisMatrix.hpp"
63 
64 // Teuchos includes
65 #include "Teuchos_ParameterList.hpp"
66 #include "Teuchos_VerboseObject.hpp"
67 #include "Teuchos_GlobalMPISession.hpp"
68 
69 #ifdef HAVE_MPI
70 # include "Tpetra_MpiPlatform.hpp"
71 #else
72 # include "Tpetra_SerialPlatform.hpp"
73 #endif
74 
75 int main(int argc, char* argv[])
76 {
77  //
78  // Get a default output stream from the Teuchos::VerboseObjectBase
79  //
80  Teuchos::RCP<Teuchos::FancyOStream>
81  out = Teuchos::VerboseObjectBase::getDefaultOStream();
82 
83  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
84 
85  //
86  // Create the scalar-type typedefs
87  //
88 
89 #ifdef HAVE_COMPLEX
90  typedef std::complex<double> ST; // Scalar-type typedef
91 #elif HAVE_COMPLEX_H
92  typedef std::complex<double> ST; // Scalar-type typedef
93 #else
94  typedef double ST; // Scalar-type typedef
95 #endif
96  typedef Teuchos::ScalarTraits<ST>::magnitudeType MT; // Magnitude-type typdef
97  typedef int OT; // Ordinal-type typdef
98 
99  ST one = Teuchos::ScalarTraits<ST>::one();
100  ST zero = Teuchos::ScalarTraits<ST>::zero();
101 
102  //
103  // Create 2D Laplacian using Tpetra::CisMatrix
104  //
105 
106 #ifdef HAVE_MPI
107  MPI_Comm mpiComm = MPI_COMM_WORLD;
108  const Tpetra::MpiPlatform<OT,OT> ordinalPlatform(mpiComm);
109  const Tpetra::MpiPlatform<OT,ST> scalarPlatform(mpiComm);
110 #else
111  const Tpetra::SerialPlatform<OT,OT> ordinalPlatform;
112  const Tpetra::SerialPlatform<OT,ST> scalarPlatform;
113 #endif
114 
115  // Declare global dimension of the linear operator
116  OT globalDim = 500;
117 
118  // Create the element space and std::vector space
119  const Tpetra::ElementSpace<OT> elementSpace(globalDim,0,ordinalPlatform);
120  const Tpetra::VectorSpace<OT,ST> vectorSpace(elementSpace,scalarPlatform);
121 
122  // Allocate the Tpetra::CisMatrix object.
123  RCP<Tpetra::CisMatrix<OT,ST> >
124  tpetra_A = rcp(new Tpetra::CisMatrix<OT,ST>(vectorSpace));
125 
126  // Allocate the Tpetra::CisMatrix preconditioning object.
127  RCP<Tpetra::CisMatrix<OT,ST> >
128  tpetra_Prec = rcp(new Tpetra::CisMatrix<OT,ST>(vectorSpace));
129 
130  // Get the indexes of the rows on this processor
131  const int numMyElements = vectorSpace.getNumMyEntries();
132  const std::vector<int> &myGlobalElements = vectorSpace.elementSpace().getMyGlobalElements();
133 
134  // Fill the local matrix entries one row at a time.
135  const ST offDiag = -1.0, diag = 2.0;
136  int numEntries; ST values[3]; int indexes[3];
137  ST prec_value[1]; int prec_index[1];
138  prec_value[0] = one/diag; // Diagonal preconditioning
139  for( int k = 0; k < numMyElements; ++k ) {
140  const int rowIndex = myGlobalElements[k];
141  prec_index[0] = rowIndex;
142  if( rowIndex == 0 ) { // First row
143  numEntries = 2;
144  values[0] = diag; values[1] = offDiag;
145  indexes[0] = 0; indexes[1] = 1;
146  }
147  else if( rowIndex == globalDim - 1 ) { // Last row
148  numEntries = 2;
149  values[0] = offDiag; values[1] = diag;
150  indexes[0] = globalDim-2; indexes[1] = globalDim-1;
151  }
152  else { // Middle rows
153  numEntries = 3;
154  values[0] = offDiag; values[1] = diag; values[2] = offDiag;
155  indexes[0] = rowIndex-1; indexes[1] = rowIndex; indexes[2] = rowIndex+1;
156  }
157  tpetra_A->submitEntries(Tpetra::Insert,rowIndex,numEntries,values,indexes);
158  tpetra_Prec->submitEntries(Tpetra::Insert,rowIndex,1,prec_value,prec_index);
159  }
160 
161  // Finish the construction of the Tpetra::CisMatrix
162  tpetra_A->fillComplete();
163  tpetra_Prec->fillComplete();
164 
165  // Create a Thyra linear operator (A) using the Tpetra::CisMatrix (tpetra_A).
166  RCP<Thyra::LinearOpBase<ST> >
167  A = Teuchos::rcp( new Thyra::TpetraLinearOp<OT,ST>(
168  Teuchos::rcp_implicit_cast<Tpetra::Operator<OT,ST> >(tpetra_A) ) );
169 
170  // Create a Thyra linear operator (Prec) using the Tpetra::CisMatrix (tpetra_Prec)
171  RCP<Thyra::LinearOpBase<ST> >
172  Prec = Teuchos::rcp( new Thyra::TpetraLinearOp<OT,ST>(
173  Teuchos::rcp_implicit_cast<Tpetra::Operator<OT,ST> >(tpetra_Prec) ) );
174 
175  // Create a Thyra default preconditioner (DefPrec) using the Thyra linear operator (Prec)
176  RCP<const Thyra::DefaultPreconditioner<ST> >
177  DefPrec = Thyra::unspecifiedPrec<ST>(Prec);
178 
179  //
180  // Set the parameters for the Belos LOWS Factory and create a parameter list.
181  // NOTE: All the code below only uses Thyra and is independent of Tpetra.
182  //
183  int blockSize = 1; // can only be 1 right now.
184  int maxIterations = globalDim;
185  int maxRestarts = 0;
186  int gmresKrylovLength = globalDim;
187  int outputFrequency = 100;
188  bool outputMaxResOnly = true;
189  MT maxResid = 1e-3;
190 
191  Teuchos::RCP<Teuchos::ParameterList>
192  belosLOWSFPL = Teuchos::rcp( new Teuchos::ParameterList() );
193 
194  belosLOWSFPL->set("Solver Type","Block GMRES");
195 
196  Teuchos::ParameterList& belosLOWSFPL_solver =
197  belosLOWSFPL->sublist("Solver Types");
198 
199  Teuchos::ParameterList& belosLOWSFPL_gmres =
200  belosLOWSFPL_solver.sublist("Block GMRES");
201 
202  belosLOWSFPL_gmres.set("Maximum Iterations",int(maxIterations));
203  belosLOWSFPL_gmres.set("Convergence Tolerance",double(maxResid));
204  belosLOWSFPL_gmres.set("Maximum Restarts",int(maxRestarts));
205  belosLOWSFPL_gmres.set("Block Size",int(blockSize));
206  belosLOWSFPL_gmres.set("Num Blocks",int(gmresKrylovLength));
207  belosLOWSFPL_gmres.set("Output Frequency",int(outputFrequency));
208  belosLOWSFPL_gmres.set("Show Maximum Residual Norm Only",bool(outputMaxResOnly));
209 
210  // Whether the linear solver succeeded.
211  // (this will be set during the residual check at the end)
212  bool success = true;
213 
214  // Number of random right-hand sides we will be solving for.
215  int numRhs = 5;
216 
217  // Get the domain space for the Thyra linear operator
218  Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domain = A->domain();
219 
220  // Create the Belos LOWS factory.
221  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ST> >
222  belosLOWSFactory = Teuchos::rcp(new Thyra::BelosLinearOpWithSolveFactory<ST>());
223 
224  // Set the parameter list to specify the behavior of the factory.
225  belosLOWSFactory->setParameterList( belosLOWSFPL );
226 
227  // Set the output stream and the verbosity level (prints to std::cout by defualt)
228  // NOTE: Set to VERB_NONE for no output from the solver.
229  belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
230 
231  // Create a BelosLinearOpWithSolve object from the Belos LOWS factory.
232  Teuchos::RCP<Thyra::LinearOpWithSolveBase<ST> >
233  nsA = belosLOWSFactory->createOp();
234 
235  // Initialize the BelosLinearOpWithSolve object with the Thyra linear operator (A)
236  // and preconditioner (DefPrec).
237  Thyra::initializeOp<ST>( *belosLOWSFactory, A, &*nsA );
238 
239  // Create a right-hand side with numRhs vectors in it and randomize it.
240  Teuchos::RCP< Thyra::MultiVectorBase<ST> >
241  b = Thyra::createMembers(domain, numRhs);
242  Thyra::seed_randomize<ST>(0);
243  Thyra::randomize(-one, one, &*b);
244 
245  // Create an initial std::vector with numRhs vectors in it and initialize it to zero.
246  Teuchos::RCP< Thyra::MultiVectorBase<ST> >
247  x = Thyra::createMembers(domain, numRhs);
248  Thyra::assign(&*x, zero);
249 
250  // Perform solve using the linear operator to get the approximate solution of Ax=b,
251  // where b is the right-hand side and x is the left-hand side.
252  Thyra::SolveStatus<ST> solveStatus;
253  solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
254 
255  // Print out status of solve.
256  *out << "\nBelos LOWS Status: "<< solveStatus << std::endl;
257 
258  //
259  // Compute residual and ST check convergence.
260  //
261  std::vector<MT> norm_b(numRhs), norm_res(numRhs);
262  Teuchos::RCP< Thyra::MultiVectorBase<ST> >
263  y = Thyra::createMembers(domain, numRhs);
264 
265  // Compute the column norms of the right-hand side b.
266  Thyra::norms_2( *b, &norm_b[0] );
267 
268  // Compute y=A*x, where x is the solution from the linear solver.
269  A->apply( Thyra::NONCONJ_ELE, *x, &*y );
270 
271  // Compute A*x-b = y-b
272  Thyra::update( -one, *b, &*y );
273 
274  // Compute the column norms of A*x-b.
275  Thyra::norms_2( *y, &norm_res[0] );
276 
277  // Print out the final relative residual norms.
278  MT rel_res = 0.0;
279  *out << "Final relative residual norms" << std::endl;
280  for (int i=0; i<numRhs; ++i) {
281  rel_res = norm_res[i]/norm_b[i];
282  if (rel_res > maxResid)
283  success = false;
284  *out << "RHS " << i+1 << " : "
285  << std::setw(16) << std::right << rel_res << std::endl;
286  }
287 
288  return ( success ? 0 : 1 );
289 }
int main(int argc, char *argv[])
double ST
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.