Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
belos_tpetra_thyra_lowsf_hb.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 reads a problem from a Harwell-Boeing (HB) file into an
44 // Tpetra::CisMatrix. This matrix is then converted into a Thyra linear operator
45 // through the Thyra-Tpetra adapters.
46 //
47 // The right-hand-side from the HB file is used instead of 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 // I/O for Harwell-Boeing files
63 #ifdef HAVE_BELOS_TRIUTILS
64 #include "iohb.h"
65 #endif
66 
67 // Teuchos includes
68 #include "Teuchos_ParameterList.hpp"
69 #include "Teuchos_VerboseObject.hpp"
70 #include "Teuchos_GlobalMPISession.hpp"
71 
72 #ifdef HAVE_MPI
73 # include "Tpetra_MpiPlatform.hpp"
74 #else
75 # include "Tpetra_SerialPlatform.hpp"
76 #endif
77 
78 // My Tpetra::Operator include
79 #include "MyOperator.hpp"
80 
81 int main(int argc, char* argv[])
82 {
83  //
84  // Get a default output stream from the Teuchos::VerboseObjectBase
85  //
86  Teuchos::RCP<Teuchos::FancyOStream>
87  out = Teuchos::VerboseObjectBase::getDefaultOStream();
88 
89  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
90 
91 #ifdef HAVE_COMPLEX
92  typedef std::complex<double> ST; // Scalar-type typedef
93 #elif HAVE_COMPLEX_H
94  typedef std::complex<double> ST; // Scalar-type typedef
95 #else
96  typedef double ST; // Scalar-type typedef
97 #endif
98 
99  typedef Teuchos::ScalarTraits<ST>::magnitudeType MT; // Magnitude-type typedef
100  typedef int OT; // Ordinal-type typedef
101  ST one = Teuchos::ScalarTraits<ST>::one();
102  ST zero = Teuchos::ScalarTraits<ST>::zero();
103 
104 #ifdef HAVE_MPI
105  MPI_Comm mpiComm = MPI_COMM_WORLD;
106  const Tpetra::MpiPlatform<OT,OT> ordinalPlatform(mpiComm);
107  const Tpetra::MpiPlatform<OT,ST> scalarPlatform(mpiComm);
108 #else
109  const Tpetra::SerialPlatform<OT,OT> ordinalPlatform;
110  const Tpetra::SerialPlatform<OT,ST> scalarPlatform;
111 #endif
112 
113  //
114  // Get the data from the HB file
115  //
116 
117  // Name of input matrix file
118  std::string matrixFile = "mhd1280b.cua";
119 
120  int info=0;
121  int dim,dim2,nnz;
122  MT *dvals;
123  int *colptr,*rowind;
124  ST *cvals;
125  nnz = -1;
126  info = readHB_newmat_double(matrixFile.c_str(),&dim,&dim2,&nnz,
127  &colptr,&rowind,&dvals);
128 
129  if (info == 0 || nnz < 0) {
130  *out << "Error reading '" << matrixFile << "'" << std::endl;
131  }
132 #ifdef HAVE_MPI
133  MPI_Finalize();
134 #endif
135 
136  // Convert interleaved doubles to std::complex values
137  cvals = new ST[nnz];
138  for (int ii=0; ii<nnz; ii++) {
139  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
140  }
141 
142  // Declare global dimension of the linear operator
143  OT globalDim = dim;
144 
145  // Create the element space and std::vector space
146  const Tpetra::ElementSpace<OT> elementSpace(globalDim,0,ordinalPlatform);
147  const Tpetra::VectorSpace<OT,ST> vectorSpace(elementSpace,scalarPlatform);
148 
149  // Create my implementation of a Tpetra::Operator
150  RCP<Tpetra::Operator<OT,ST> >
151  tpetra_A = rcp( new MyOperator<OT,ST>(vectorSpace,dim,colptr,nnz,rowind,cvals) );
152 
153  // Create a Thyra linear operator (A) using the Tpetra::CisMatrix (tpetra_A).
154  RCP<Thyra::LinearOpBase<ST> >
155  A = Teuchos::rcp( new Thyra::TpetraLinearOp<OT,ST>(tpetra_A) );
156 
157  //
158  // Set the parameters for the Belos LOWS Factory and create a parameter list.
159  //
160  int blockSize = 1;
161  int maxIterations = globalDim;
162  int maxRestarts = 15;
163  int gmresKrylovLength = 50;
164  int outputFrequency = 100;
165  bool outputMaxResOnly = true;
166  MT maxResid = 1e-5;
167 
168  Teuchos::RCP<Teuchos::ParameterList>
169  belosLOWSFPL = Teuchos::rcp( new Teuchos::ParameterList() );
170 
171  belosLOWSFPL->set("Solver Type","Block GMRES");
172 
173  Teuchos::ParameterList& belosLOWSFPL_solver =
174  belosLOWSFPL->sublist("Solver Types");
175 
176  Teuchos::ParameterList& belosLOWSFPL_gmres =
177  belosLOWSFPL_solver.sublist("Block GMRES");
178 
179  belosLOWSFPL_gmres.set("Maximum Iterations",int(maxIterations));
180  belosLOWSFPL_gmres.set("Convergence Tolerance",MT(maxResid));
181  belosLOWSFPL_gmres.set("Maximum Restarts",int(maxRestarts));
182  belosLOWSFPL_gmres.set("Block Size",int(blockSize));
183  belosLOWSFPL_gmres.set("Num Blocks",int(gmresKrylovLength));
184  belosLOWSFPL_gmres.set("Output Frequency",int(outputFrequency));
185  belosLOWSFPL_gmres.set("Show Maximum Residual Norm Only",bool(outputMaxResOnly));
186 
187  // Whether the linear solver succeeded.
188  // (this will be set during the residual check at the end)
189  bool success = true;
190 
191  // Number of random right-hand sides we will be solving for.
192  int numRhs = 1;
193 
194  // Get the domain space for the Thyra linear operator
195  Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domain = A->domain();
196 
197  // Create the Belos LOWS factory.
198  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ST> >
199  belosLOWSFactory = Teuchos::rcp(new Thyra::BelosLinearOpWithSolveFactory<ST>());
200 
201  // Set the parameter list to specify the behavior of the factory.
202  belosLOWSFactory->setParameterList( belosLOWSFPL );
203 
204  // Set the output stream and the verbosity level (prints to std::cout by defualt)
205  // NOTE: Set to VERB_NONE for no output from the solver.
206  belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
207 
208  // Create a BelosLinearOpWithSolve object from the Belos LOWS factory.
209  Teuchos::RCP<Thyra::LinearOpWithSolveBase<ST> >
210  nsA = belosLOWSFactory->createOp();
211 
212  // Initialize the BelosLinearOpWithSolve object with the Thyra linear operator.
213  Thyra::initializeOp<ST>( *belosLOWSFactory, A, &*nsA );
214 
215  // Create a right-hand side with numRhs vectors in it.
216  Teuchos::RCP< Thyra::MultiVectorBase<ST> >
217  b = Thyra::createMembers(domain, numRhs);
218 
219  // Create an initial std::vector with numRhs vectors in it and initialize it to one.
220  Teuchos::RCP< Thyra::MultiVectorBase<ST> >
221  x = Thyra::createMembers(domain, numRhs);
222  Thyra::assign(&*x, one);
223 
224  // Initialize the right-hand side so that the solution is a std::vector of ones.
225  A->apply( Thyra::NONCONJ_ELE, *x, &*b );
226  Thyra::assign(&*x, zero);
227 
228  // Perform solve using the linear operator to get the approximate solution of Ax=b,
229  // where b is the right-hand side and x is the left-hand side.
230  Thyra::SolveStatus<ST> solveStatus;
231  solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
232 
233  // Print out status of solve.
234  *out << "\nBelos LOWS Status: "<< solveStatus << std::endl;
235 
236  //
237  // Compute residual and ST check convergence.
238  //
239  std::vector<MT> norm_b(numRhs), norm_res(numRhs);
240  Teuchos::RCP< Thyra::MultiVectorBase<ST> >
241  y = Thyra::createMembers(domain, numRhs);
242 
243  // Compute the column norms of the right-hand side b.
244  Thyra::norms_2( *b, &norm_b[0] );
245 
246  // Compute y=A*x, where x is the solution from the linear solver.
247  A->apply( Thyra::NONCONJ_ELE, *x, &*y );
248 
249  // Compute A*x-b = y-b
250  Thyra::update( -one, *b, &*y );
251 
252  // Compute the column norms of A*x-b.
253  Thyra::norms_2( *y, &norm_res[0] );
254 
255  // Print out the final relative residual norms.
256  MT rel_res = 0.0;
257  *out << "Final relative residual norms" << std::endl;
258  for (int i=0; i<numRhs; ++i) {
259  rel_res = norm_res[i]/norm_b[i];
260  if (rel_res > maxResid)
261  success = false;
262  *out << "RHS " << i+1 << " : "
263  << std::setw(16) << std::right << rel_res << std::endl;
264  }
265 
266  return ( success ? 0 : 1 );
267 }
int main(int argc, char *argv[])
double ST
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.
Simple example of a user&#39;s defined Tpetra::Operator class.
Definition: MyOperator.hpp:63