Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
compare_solvers.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Amesos: An Interface to Direct Solvers
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25//
26// ***********************************************************************
27// @HEADER
28
29#include "Amesos_ConfigDefs.h"
30// This example needs Galeri to generate the linear system.
31#ifdef HAVE_MPI
32#include "mpi.h"
33#include "Epetra_MpiComm.h"
34#else
35#include "Epetra_SerialComm.h"
36#endif
37#include "Epetra_Vector.h"
38#include "Epetra_Time.h"
39#include "Epetra_RowMatrix.h"
40#include "Epetra_CrsMatrix.h"
41#include "Amesos.h"
42#include "Amesos_BaseSolver.h"
43#include "Teuchos_ParameterList.hpp"
44#include "Galeri_Maps.h"
45#include "Galeri_CrsMatrices.h"
46
47using namespace Teuchos;
48using namespace Galeri;
49
50#include "Trilinos_Util.h"
51#include "Trilinos_Util_ReadMatrixMarket2Epetra.h"
52#include "CrsMatrixTranspose.h"
53#include "Teuchos_RCP.hpp"
54#include "Epetra_Export.h"
55
56int MyCreateCrsMatrix( const char *in_filename, const Epetra_Comm &Comm,
57 Epetra_Map *& readMap,
58 const bool transpose, const bool distribute,
59 bool& symmetric, Epetra_CrsMatrix *& Matrix ) {
60
61 Epetra_CrsMatrix * readA = 0;
62 Epetra_Vector * readx = 0;
63 Epetra_Vector * readb = 0;
64 Epetra_Vector * readxexact = 0;
65
66 //
67 // This hack allows TestOptions to be run from either the test/TestOptions/ directory or from
68 // the test/ directory (as it is in nightly testing and in make "run-tests")
69 //
70 FILE *in_file = fopen( in_filename, "r");
71
72 const char *filename;
73 if (in_file == NULL )
74 filename = &in_filename[1] ; // Strip off ithe "." from
75 // "../" and try again
76 else {
77 filename = in_filename ;
78 fclose( in_file );
79 }
80
81 symmetric = false ;
82 std::string FileName (filename);
83
84 int FN_Size = FileName.size() ;
85 std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
86 std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
87
88 if ( LastFiveBytes == ".TimD" ) {
89 // Call routine to read in a file in a Zero Based File in Tim Davis format
90 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx,
91 readb, readxexact, false, true, true ) );
92 symmetric = false;
93 } else {
94 if ( LastFiveBytes == ".triU" ) {
95 // Call routine to read in unsymmetric Triplet matrix
96 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx,
97 readb, readxexact) );
98 symmetric = false;
99 } else {
100 if ( LastFiveBytes == ".triS" ) {
101 // Call routine to read in symmetric Triplet matrix
102 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, true, Comm, readMap, readA, readx,
103 readb, readxexact) );
104 symmetric = true;
105 } else {
106 if ( LastFourBytes == ".mtx" ) {
107 EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( filename, Comm, readMap,
108 readA, readx, readb, readxexact) );
109 FILE* in_file = fopen( filename, "r");
110 assert (in_file != NULL) ; // Checked in Trilinos_Util_CountMatrixMarket()
111 const int BUFSIZE = 800 ;
112 char buffer[BUFSIZE] ;
113 fgets( buffer, BUFSIZE, in_file ) ; // Pick symmetry info off of this string
114 std::string headerline1 = buffer;
115#ifdef TFLOP
116 if ( headerline1.find("symmetric") < BUFSIZE ) symmetric = true;
117#else
118 if ( headerline1.find("symmetric") != std::string::npos) symmetric = true;
119
120#endif
121 fclose(in_file);
122
123 } else {
124 // Call routine to read in HB problem
125 Trilinos_Util_ReadHb2Epetra( filename, Comm, readMap, readA, readx,
126 readb, readxexact) ;
127 if ( LastFourBytes == ".rsa" ) symmetric = true ;
128 }
129 }
130 }
131 }
132
133
134 if ( readb ) delete readb;
135 if ( readx ) delete readx;
136 if ( readxexact ) delete readxexact;
137
138 Epetra_CrsMatrix *serialA ;
139 Epetra_CrsMatrix *transposeA;
140
141 if ( transpose ) {
142 transposeA = new Epetra_CrsMatrix( Copy, *readMap, 0 );
143 assert( CrsMatrixTranspose( readA, transposeA ) == 0 );
144 serialA = transposeA ;
145 delete readA;
146 readA = 0 ;
147 } else {
148 serialA = readA ;
149 }
150
151 assert( (void *) &serialA->Graph() ) ;
152 assert( (void *) &serialA->RowMap() ) ;
153 assert( serialA->RowMap().SameAs(*readMap) ) ;
154
155 if ( distribute ) {
156 // Create uniform distributed map
157 Epetra_Map DistMap(readMap->NumGlobalElements(), 0, Comm);
158
159 // Create Exporter to distribute read-in matrix and vectors
160 Epetra_Export exporter( *readMap, DistMap );
161
162 Epetra_CrsMatrix *Amat = new Epetra_CrsMatrix( Copy, DistMap, 0 );
163 Amat->Export(*serialA, exporter, Add);
164 assert(Amat->FillComplete()==0);
165
166 Matrix = Amat;
167 //
168 // Make sure that deleting Amat->RowMap() will delete map
169 //
170 // Bug: We can't manage to delete map his way anyway,
171 // and this fails on tranposes, so for now I just accept
172 // the memory loss.
173 // assert( &(Amat->RowMap()) == map ) ;
174 delete readMap;
175 readMap = 0 ;
176 delete serialA;
177 } else {
178
179 Matrix = serialA;
180 }
181
182
183 return 0;
184}
185
186
187
188
189
190// ===================== //
191// M A I N D R I V E R //
192// ===================== //
193//
194// This example compares all the available Amesos solvers
195// for the solution of the same linear system.
196//
197// The example can be run in serial and in parallel.
198//
199// Author: Marzio Sala, SNL 9214
200// Last modified: Oct-05.
201
202int main(int argc, char *argv[])
203{
204#ifdef HAVE_MPI
205 MPI_Init(&argc, &argv);
206 Epetra_MpiComm Comm(MPI_COMM_WORLD);
207#else
208 Epetra_SerialComm Comm;
209#endif
210
211 bool verbose = (Comm.MyPID() == 0);
212 double TotalResidual = 0.0;
213
214 // Create the Map, defined as a grid, of size nx x ny x nz,
215 // subdivided into mx x my x mz cubes, each assigned to a
216 // different processor.
217
218#ifndef FILENAME_SPECIFIED_ON_COMMAND_LINE
219 ParameterList GaleriList;
220 GaleriList.set("nx", 4);
221 GaleriList.set("ny", 4);
222 GaleriList.set("nz", 4 * Comm.NumProc());
223 GaleriList.set("mx", 1);
224 GaleriList.set("my", 1);
225 GaleriList.set("mz", Comm.NumProc());
226 Epetra_Map* Map = CreateMap("Cartesian3D", Comm, GaleriList);
227
228 // Create a matrix, in this case corresponding to a 3D Laplacian
229 // discretized using a classical 7-point stencil. Please refer to
230 // the Galeri documentation for an overview of available matrices.
231 //
232 // NOTE: matrix must be symmetric if DSCPACK is used.
233
234 Epetra_CrsMatrix* Matrix = CreateCrsMatrix("Laplace3D", Map, GaleriList);
235#else
236 bool transpose = false ;
237 bool distribute = false ;
238 bool symmetric ;
239 Epetra_CrsMatrix *Matrix = 0 ;
240 Epetra_Map *Map = 0 ;
241 MyCreateCrsMatrix( argv[1], Comm, Map, transpose, distribute, symmetric, Matrix ) ;
242
243#endif
244
245 // build vectors, in this case with 1 vector
246 Epetra_MultiVector LHS(*Map, 1);
247 Epetra_MultiVector RHS(*Map, 1);
248
249 // create a linear problem object
250 Epetra_LinearProblem Problem(Matrix, &LHS, &RHS);
251
252 // use this list to set up parameters, now it is required
253 // to use all the available processes (if supported by the
254 // underlying solver). Uncomment the following two lines
255 // to let Amesos print out some timing and status information.
256 ParameterList List;
257 List.set("PrintTiming",true);
258 List.set("PrintStatus",true);
259 List.set("MaxProcs",Comm.NumProc());
260
261 std::vector<std::string> SolverType;
262 SolverType.push_back("Amesos_Paraklete");
263 SolverType.push_back("Amesos_Klu");
264 Comm.Barrier() ;
265#if 1
266 SolverType.push_back("Amesos_Lapack");
267 SolverType.push_back("Amesos_Umfpack");
268 SolverType.push_back("Amesos_Pardiso");
269 SolverType.push_back("Amesos_Taucs");
270 SolverType.push_back("Amesos_Superlu");
271 SolverType.push_back("Amesos_Superludist");
272 SolverType.push_back("Amesos_Mumps");
273 SolverType.push_back("Amesos_Dscpack");
274 SolverType.push_back("Amesos_Scalapack");
275#endif
276 Epetra_Time Time(Comm);
277
278 // this is the Amesos factory object that will create
279 // a specific Amesos solver.
280 Amesos Factory;
281
282 // Cycle over all solvers.
283 // Only installed solvers will be tested.
284 for (unsigned int i = 0 ; i < SolverType.size() ; ++i)
285 {
286 // Check whether the solver is available or not
287 if (Factory.Query(SolverType[i]))
288 {
289 // 1.- set exact solution (constant vector)
290 LHS.PutScalar(1.0);
291
292 // 2.- create corresponding rhs
293 Matrix->Multiply(false, LHS, RHS);
294
295 // 3.- randomize solution vector
296 LHS.Random();
297
298 // 4.- create the amesos solver object
299 Amesos_BaseSolver* Solver = Factory.Create(SolverType[i], Problem);
300 assert (Solver != 0);
301
302 Solver->SetParameters(List);
303 Solver->SetUseTranspose( true) ;
304
305 // 5.- factorize and solve
306
307
308 Comm.Barrier() ;
309 if (verbose)
310 std::cout << std::endl
311 << "Solver " << SolverType[i]
312 << ", verbose = " << verbose << std::endl ;
313 Comm.Barrier() ;
314
315
316 Time.ResetStartTime();
318 if (verbose)
319 std::cout << std::endl
320 << "Solver " << SolverType[i]
321 << ", symbolic factorization time = "
322 << Time.ElapsedTime() << std::endl;
323 Comm.Barrier() ;
324
326 if (verbose)
327 std::cout << "Solver " << SolverType[i]
328 << ", numeric factorization time = "
329 << Time.ElapsedTime() << std::endl;
330 Comm.Barrier() ;
331
332 AMESOS_CHK_ERR(Solver->Solve());
333 if (verbose)
334 std::cout << "Solver " << SolverType[i]
335 << ", solve time = "
336 << Time.ElapsedTime() << std::endl;
337 Comm.Barrier() ;
338
339 // 6.- compute difference between exact solution and Amesos one
340 // (there are other ways of doing this in Epetra, but let's
341 // keep it simple)
342 double d = 0.0, d_tot = 0.0;
343 for (int j = 0 ; j< LHS.Map().NumMyElements() ; ++j)
344 d += (LHS[0][j] - 1.0) * (LHS[0][j] - 1.0);
345
346 Comm.SumAll(&d,&d_tot,1);
347 if (verbose)
348 std::cout << "Solver " << SolverType[i] << ", ||x - x_exact||_2 = "
349 << sqrt(d_tot) << std::endl;
350
351 // 7.- delete the object
352 delete Solver;
353
354 TotalResidual += d_tot;
355 }
356 }
357
358 delete Matrix;
359 delete Map;
360
361 if (TotalResidual > 1e-9)
362 exit(EXIT_FAILURE);
363
364#ifdef HAVE_MPI
365 MPI_Finalize();
366#endif
367
368 return(EXIT_SUCCESS);
369} // end of main()
static bool verbose
Definition: Amesos.cpp:67
#define AMESOS_CHK_ERR(a)
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
#define EPETRA_CHK_ERR(xxx)
int CreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)
Definition: TestOptions.cpp:81
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
virtual int Solve()=0
Solves A X = B (or AT x = B)
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition: Amesos.h:44
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Definition: Amesos.cpp:69
bool Query(const char *ClassType)
Queries whether a given interface is available or not.
Definition: Amesos.cpp:193
int main(int argc, char *argv[])
int MyCreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)