Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
belos_epetra_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// Epetra_CrsMatrix. This matrix is then converted into a Thyra linear operator
45// through the Thyra-Epetra adapters.
46// The right-hand-side from the HB file is used instead of random vectors.
47// The initial guesses are all set to zero.
48//
49//
50
51// Thyra includes
52#include "Thyra_BelosLinearOpWithSolveFactory.hpp"
53#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
54#include "Thyra_EpetraLinearOp.hpp"
55
56// Epetra includes
57#include "Epetra_SerialComm.h"
58#include "EpetraExt_readEpetraLinearSystem.h"
59
60// Ifpack includes
61#ifdef HAVE_BELOS_IFPACK
63#endif
64
65// Teuchos includes
66#include "Teuchos_ParameterList.hpp"
67#include "Teuchos_VerboseObject.hpp"
68
69int main(int argc, char* argv[])
70{
71 //
72 // Get a default output stream from the Teuchos::VerboseObjectBase
73 //
74 Teuchos::RCP<Teuchos::FancyOStream>
75 out = Teuchos::VerboseObjectBase::getDefaultOStream();
76 //
77 // Set the parameters for the Belos LOWS Factory and create a parameter list.
78 //
79 int blockSize = 2;
80 int maxIterations = 400;
81 int maxRestarts = 25;
82 int gmresKrylovLength = 25;
83 int outputFrequency = 1;
84 bool outputMaxResOnly = true;
85 double maxResid = 1e-6;
86
87 Teuchos::RCP<Teuchos::ParameterList>
88 belosLOWSFPL = Teuchos::rcp( new Teuchos::ParameterList() );
89
90 belosLOWSFPL->set("Solver Type","Block GMRES");
91
92 Teuchos::ParameterList& belosLOWSFPL_solver =
93 belosLOWSFPL->sublist("Solver Types");
94
95 Teuchos::ParameterList& belosLOWSFPL_gmres =
96 belosLOWSFPL_solver.sublist("Block GMRES");
97
98 belosLOWSFPL_gmres.set("Maximum Iterations",int(maxIterations));
99 belosLOWSFPL_gmres.set("Convergence Tolerance",double(maxResid));
100 belosLOWSFPL_gmres.set("Maximum Restarts",int(maxRestarts));
101 belosLOWSFPL_gmres.set("Block Size",int(blockSize));
102 belosLOWSFPL_gmres.set("Num Blocks",int(gmresKrylovLength));
103 belosLOWSFPL_gmres.set("Output Frequency",int(outputFrequency));
104 belosLOWSFPL_gmres.set("Show Maximum Residual Norm Only",bool(outputMaxResOnly));
105
106#ifdef HAVE_BELOS_IFPACK
107 //
108 // Set the parameters for the Ifpack Preconditioner Factory and create parameter list
109 //
110 Teuchos::ParameterList &ifpackPFSL = belosLOWSFPL->sublist("IfpackPreconditionerFactory");
111
112 ifpackPFSL.set("Overlap",int(2));
113 ifpackPFSL.set("Prec Type","ILUT");
114#endif
115
116 // Whether the linear solver succeeded.
117 // (this will be set during the residual check at the end)
118 bool success = true;
119
120 // Number of random right-hand sides we will be solving for.
121 int numRhs = 5;
122
123 // Name of input matrix file
124 std::string matrixFile = "orsirr1.hb";
125
126 // Read in the matrix file (can be *.mtx, *.hb, etc.)
127 Epetra_SerialComm comm;
128 Teuchos::RCP<Epetra_CrsMatrix> epetra_A;
129 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
130
131 // Create a Thyra linear operator (A) using the Epetra_CrsMatrix (epetra_A).
132 Teuchos::RCP<const Thyra::LinearOpBase<double> >
133 A = Thyra::epetraLinearOp(epetra_A);
134
135 // Get the domain space for the Thyra linear operator
136 Teuchos::RCP<const Thyra::VectorSpaceBase<double> > domain = A->domain();
137
138 // Create the Belos LOWS factory.
139 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
140 belosLOWSFactory = Teuchos::rcp(new Thyra::BelosLinearOpWithSolveFactory<double>());
141
142#ifdef HAVE_BELOS_IFPACK
143
144 // Set the preconditioner factory for the LOWS factory.
145 belosLOWSFactory->setPreconditionerFactory(
146 Teuchos::rcp(new Thyra::IfpackPreconditionerFactory())
147 ,"IfpackPreconditionerFactory"
148 );
149#endif
150
151 // Set the parameter list to specify the behavior of the factory.
152 belosLOWSFactory->setParameterList( belosLOWSFPL );
153
154 // Set the output stream and the verbosity level (prints to std::cout by defualt)
155 belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
156
157 // Create a BelosLinearOpWithSolve object from the Belos LOWS factory.
158 Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> >
159 nsA = belosLOWSFactory->createOp();
160
161 // Initialize the BelosLinearOpWithSolve object with the Thyra linear operator.
162 Thyra::initializeOp<double>( *belosLOWSFactory, A, &*nsA );
163
164 // Create a right-hand side with numRhs vectors in it and randomize it.
165 Teuchos::RCP< Thyra::MultiVectorBase<double> >
166 b = Thyra::createMembers(domain, numRhs);
167 Thyra::seed_randomize<double>(0);
168 Thyra::randomize(-1.0, 1.0, &*b);
169
170 // Create an initial std::vector with numRhs vectors in it and initialize it to zero.
171 Teuchos::RCP< Thyra::MultiVectorBase<double> >
172 x = Thyra::createMembers(domain, numRhs);
173 Thyra::assign(&*x, 0.0);
174
175 // Perform solve using the linear operator to get the approximate solution of Ax=b,
176 // where b is the right-hand side and x is the left-hand side.
177 Thyra::SolveStatus<double> solveStatus;
178 solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
179
180 // Print out status of solve.
181 *out << "\nBelos LOWS Status: "<< solveStatus << std::endl;
182
183 //
184 // Compute residual and double check convergence.
185 //
186 std::vector<double> norm_b(numRhs), norm_res(numRhs);
187 Teuchos::RCP< Thyra::MultiVectorBase<double> >
188 y = Thyra::createMembers(domain, numRhs);
189
190 // Compute the column norms of the right-hand side b.
191 Thyra::norms_2( *b, &norm_b[0] );
192
193 // Compute y=A*x, where x is the solution from the linear solver.
194 A->apply( Thyra::NONCONJ_ELE, *x, &*y );
195
196 // Compute A*x-b = y-b
197 Thyra::update( -1.0, *b, &*y );
198
199 // Compute the column norms of A*x-b.
200 Thyra::norms_2( *y, &norm_res[0] );
201
202 // Print out the final relative residual norms.
203 double rel_res = 0.0;
204 *out << "Final relative residual norms" << std::endl;
205 for (int i=0; i<numRhs; ++i) {
206 rel_res = norm_res[i]/norm_b[i];
207 if (rel_res > maxResid)
208 success = false;
209 *out << "RHS " << i+1 << " : "
210 << std::setw(16) << std::right << rel_res << std::endl;
211 }
212
213 return ( success ? 0 : 1 );
214}
int main(int argc, char *argv[])
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.
Concrete preconditioner factory subclass based on Ifpack.