Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Amesos_TestSolver.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Amesos: Direct Sparse Solver Package
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#include "Teuchos_ParameterList.hpp"
31#include <string>
32// #include "Trilinos_Util_ReadTriples2Epetra.h"
33#include "Trilinos_Util.h"
34// #include "Trilinos_Util_ReadMatrixMarket2Epetra.h"
35#include "Epetra_LocalMap.h"
36#include "Epetra_Map.h"
37#include "Epetra_Vector.h"
38#include "Epetra_Export.h"
39#include "Epetra_CrsMatrix.h"
40#include "Epetra_LinearProblem.h"
41#include "Amesos_Time.h"
42
43#ifdef TEST_SPOOLES
44#include "SpoolesOO.h"
45#endif
46#ifdef HAVE_AMESOS_SUPERLU
47#include "Amesos_Superlu.h"
48#endif
49#ifdef HAVE_AMESOS_SUPERLUDIST
50#include "Amesos_Superludist.h"
51#endif
52#ifdef HAVE_AMESOS_SLUD
53#include "SuperludistOO.h"
54#endif
55#ifdef HAVE_AMESOS_SLUS
56#include "Epetra_SLU.h"
57#endif
58#ifdef HAVE_AMESOS_SLUD2
59#include "Superludist2_OO.h"
60#endif
61#ifdef TEST_AZTEC
62#include "AztecOO.h"
63#endif
64#ifdef HAVE_AMESOS_DSCPACK
65#include "Amesos_Dscpack.h"
66#endif
67#ifdef HAVE_AMESOS_LAPACK
68#include "Amesos_Lapack.h"
69#endif
70#ifdef HAVE_AMESOS_UMFPACK
71#include "Amesos_Umfpack.h"
72#endif
73#ifdef HAVE_AMESOS_KLU
74#include "Amesos_Klu.h"
75#endif
76#ifdef HAVE_AMESOS_SCALAPACK
77#include "Amesos_Scalapack.h"
78#endif
79#ifdef HAVE_AMESOS_TAUCS
80#include "Amesos_Taucs.h"
81#endif
82#ifdef HAVE_AMESOS_PARDISO
83#include "Amesos_Pardiso.h"
84#endif
85#ifdef HAVE_AMESOS_PARAKLETE
86#include "Amesos_Paraklete.h"
87#endif
88#if defined(HAVE_AMESOS_MUMPS) && defined(HAVE_MPI)
89#include "Amesos_Mumps.h"
90#endif
91
92#include "Amesos_TestSolver.h"
93#include "CrsMatrixTranspose.h"
95
96//
97// Amesos_TestSolver.cpp reads in a matrix in Harwell-Boeing format,
98// calls one of the sparse direct solvers and computes the error
99// and residual.
100//
101// It reads the matrix in on a single processor and can pass that
102// matrix to the solver or it can convert that matrix to a
103// distributed matrix and pass the distributed matrix to the solver.
104//
105// Amesos_TestSolver can test either A x = b or A^T x = b.
106// This can be a bit confusing because sparse direct solvers
107// use compressed column storage - the transpose of Trilinos'
108// sparse row storage.
109//
110// Matrices:
111// readA - Serial. As read from the file.
112// transposeA - Serial. The transpose of readA.
113// serialA - if (transpose) then transposeA else readA
114// distributedA - readA distributed to all processes
115// passA - if ( distributed ) then distributedA else serialA
116//
117//
118
119int Amesos_TestSolver( Epetra_Comm &Comm, char *matrix_file,
120 SparseSolverType SparseSolver,
121 bool transpose,
122 int special, AMESOS_MatrixType matrix_type ) {
123
124
125 Epetra_Map * readMap;
126 Epetra_CrsMatrix * readA;
127 Epetra_Vector * readx;
128 Epetra_Vector * readb;
129 Epetra_Vector * readxexact;
130
131 std::string FileName = matrix_file ;
132 int FN_Size = FileName.size() ;
133 std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
134 std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
135 bool NonContiguousMap = false;
136
137 if ( LastFiveBytes == ".triU" ) {
138 // Call routine to read in unsymmetric Triplet matrix
139 NonContiguousMap = true;
140 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, false, Comm, readMap, readA, readx,
141 readb, readxexact, NonContiguousMap ) );
142 } else {
143 if ( LastFiveBytes == ".triS" ) {
144 NonContiguousMap = true;
145 // Call routine to read in symmetric Triplet matrix
146 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, true, Comm, readMap, readA, readx,
147 readb, readxexact, NonContiguousMap ) );
148 } else {
149 if ( LastFourBytes == ".mtx" ) {
150 EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( matrix_file, Comm, readMap,
151 readA, readx, readb, readxexact) );
152 } else {
153 // Call routine to read in HB problem
154 Trilinos_Util_ReadHb2Epetra( matrix_file, Comm, readMap, readA, readx,
155 readb, readxexact) ;
156 }
157 }
158 }
159
160 Epetra_CrsMatrix transposeA(Copy, *readMap, 0);
161 Epetra_CrsMatrix *serialA ;
162
163 if ( transpose ) {
164 assert( CrsMatrixTranspose( readA, &transposeA ) == 0 );
165 serialA = &transposeA ;
166 } else {
167 serialA = readA ;
168 }
169
170 Epetra_RowMatrix * passA = 0;
171 Epetra_Vector * passx = 0;
172 Epetra_Vector * passb = 0;
173 Epetra_Vector * passxexact = 0;
174 Epetra_Vector * passresid = 0;
175 Epetra_Vector * passtmp = 0;
176
177 // Create uniform distributed map
178 Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
179 Epetra_Map* map_;
180
181 if( NonContiguousMap ) {
182 //
183 // map gives us NumMyElements and MyFirstElement;
184 //
185 int NumGlobalElements = readMap->NumGlobalElements();
186 int NumMyElements = map.NumMyElements();
187 int MyFirstElement = map.MinMyGID();
188 std::vector<int> MapMap_( NumGlobalElements );
189 readMap->MyGlobalElements( &MapMap_[0] ) ;
190 Comm.Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ;
191 map_ = new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);
192 } else {
193 map_ = new Epetra_Map( map ) ;
194 }
195
196
197 Epetra_CrsMatrix A(Copy, *map_, 0);
198
199
200 const Epetra_Map &OriginalMap = serialA->RowMatrixRowMap() ;
201 assert( OriginalMap.SameAs(*readMap) );
202 Epetra_Export exporter(OriginalMap, *map_);
203 Epetra_Export exporter2(OriginalMap, *map_);
204 Epetra_Export MatrixExporter(OriginalMap, *map_);
205 Epetra_CrsMatrix AwithDiag(Copy, *map_, 0);
206
207 Epetra_Vector x(*map_);
208 Epetra_Vector b(*map_);
209 Epetra_Vector xexact(*map_);
210 Epetra_Vector resid(*map_);
211 Epetra_Vector readresid(*readMap);
212 Epetra_Vector tmp(*map_);
213 Epetra_Vector readtmp(*readMap);
214
215 // Epetra_Vector xcomp(*map_); // X as computed by the solver
216 bool distribute_matrix = ( matrix_type == AMESOS_Distributed ) ;
217 if ( distribute_matrix ) {
218 // Create Exporter to distribute read-in matrix and vectors
219 //
220 // Initialize x, b and xexact to the values read in from the file
221 //
222 x.Export(*readx, exporter, Add);
223 b.Export(*readb, exporter, Add);
224 xexact.Export(*readxexact, exporter, Add);
225 Comm.Barrier();
226
227 A.Export(*serialA, exporter, Add);
228 assert(A.FillComplete()==0);
229
230 Comm.Barrier();
231
232 passA = &A;
233
234 passx = &x;
235 passb = &b;
236 passxexact = &xexact;
237 passresid = &resid;
238 passtmp = &tmp;
239
240 } else {
241
242 passA = serialA;
243 passx = readx;
244 passb = readb;
245 passxexact = readxexact;
246 passresid = &readresid;
247 passtmp = &readtmp;
248 }
249
250 Epetra_MultiVector CopyB( *passb ) ;
251
252
253 double Anorm = passA->NormInf() ;
254 SparseDirectTimingVars::SS_Result.Set_Anorm(Anorm) ;
255
256 Epetra_LinearProblem Problem( (Epetra_RowMatrix *) passA,
257 (Epetra_MultiVector *) passx,
258 (Epetra_MultiVector *) passb );
259
260
261 for ( int i = 0; i < 1+special ; i++ ) {
262 Epetra_Time TotalTime( Comm ) ;
263
264 if ( false ) {
265 // TEST_UMFPACK is never set by configure
266#ifdef HAVE_AMESOS_SUPERLUDIST
267 } else if ( SparseSolver == SUPERLUDIST ) {
268 Teuchos::ParameterList ParamList ;
269 ParamList.set( "MaxProcs", -3 );
270 Amesos_Superludist A_Superludist( Problem ) ;
271
272 //ParamList.set( "Redistribute", true );
273 //ParamList.set( "AddZeroToDiag", true );
274 Teuchos::ParameterList& SuperludistParams = ParamList.sublist("Superludist") ;
275 ParamList.set( "MaxProcs", -3 );
276
277 EPETRA_CHK_ERR( A_Superludist.SetParameters( ParamList ) );
278 EPETRA_CHK_ERR( A_Superludist.SetUseTranspose( transpose ) );
279 EPETRA_CHK_ERR( A_Superludist.SymbolicFactorization( ) );
280 EPETRA_CHK_ERR( A_Superludist.NumericFactorization( ) );
281 EPETRA_CHK_ERR( A_Superludist.Solve( ) );
282#endif
283#ifdef HAVE_AMESOS_DSCPACK
284 } else if ( SparseSolver == DSCPACK ) {
285
286 Teuchos::ParameterList ParamList ;
287 ParamList.set( "MaxProcs", -3 );
288
289 Amesos_Dscpack A_dscpack( Problem ) ;
290 EPETRA_CHK_ERR( A_dscpack.SetParameters( ParamList ) );
291 EPETRA_CHK_ERR( A_dscpack.SymbolicFactorization( ) );
292 EPETRA_CHK_ERR( A_dscpack.NumericFactorization( ) );
293 EPETRA_CHK_ERR( A_dscpack.Solve( ) );
294#endif
295#ifdef HAVE_AMESOS_SCALAPACK
296 } else if ( SparseSolver == SCALAPACK ) {
297
298 Teuchos::ParameterList ParamList ;
299 Amesos_Scalapack A_scalapack( Problem ) ;
300 ParamList.set( "MaxProcs", -3 );
301 EPETRA_CHK_ERR( A_scalapack.SetParameters( ParamList ) );
302 EPETRA_CHK_ERR( A_scalapack.SetUseTranspose( transpose ) );
303 EPETRA_CHK_ERR( A_scalapack.SymbolicFactorization( ) );
304 EPETRA_CHK_ERR( A_scalapack.NumericFactorization( ) );
305 EPETRA_CHK_ERR( A_scalapack.Solve( ) );
306
307#endif
308#ifdef HAVE_AMESOS_TAUCS
309 } else if ( SparseSolver == TAUCS ) {
310
311 Teuchos::ParameterList ParamList ;
312 Amesos_Taucs A_taucs( Problem ) ;
313 ParamList.set( "MaxProcs", -3 );
314 EPETRA_CHK_ERR( A_taucs.SetParameters( ParamList ) );
315 EPETRA_CHK_ERR( A_taucs.SetUseTranspose( transpose ) );
316 EPETRA_CHK_ERR( A_taucs.SymbolicFactorization( ) );
317 EPETRA_CHK_ERR( A_taucs.NumericFactorization( ) );
318 EPETRA_CHK_ERR( A_taucs.Solve( ) );
319
320#endif
321#ifdef HAVE_AMESOS_PARDISO
322 } else if ( SparseSolver == PARDISO ) {
323
324 Teuchos::ParameterList ParamList ;
325 Amesos_Pardiso A_pardiso( Problem ) ;
326 ParamList.set( "MaxProcs", -3 );
327 EPETRA_CHK_ERR( A_pardiso.SetParameters( ParamList ) );
328 EPETRA_CHK_ERR( A_pardiso.SetUseTranspose( transpose ) );
329 EPETRA_CHK_ERR( A_pardiso.SymbolicFactorization( ) );
330 EPETRA_CHK_ERR( A_pardiso.NumericFactorization( ) );
331 EPETRA_CHK_ERR( A_pardiso.Solve( ) );
332
333#endif
334#ifdef HAVE_AMESOS_PARAKLETE
335 } else if ( SparseSolver == PARAKLETE ) {
336
337 Teuchos::ParameterList ParamList ;
338 Amesos_Paraklete A_paraklete( Problem ) ;
339 ParamList.set( "MaxProcs", -3 );
340 EPETRA_CHK_ERR( A_paraklete.SetParameters( ParamList ) );
341 EPETRA_CHK_ERR( A_paraklete.SetUseTranspose( transpose ) );
342 EPETRA_CHK_ERR( A_paraklete.SymbolicFactorization( ) );
343 EPETRA_CHK_ERR( A_paraklete.NumericFactorization( ) );
344 EPETRA_CHK_ERR( A_paraklete.Solve( ) );
345
346#endif
347#if defined(HAVE_AMESOS_MUMPS) && defined(HAVE_MPI)
348 } else if ( SparseSolver == MUMPS ) {
349
350 Teuchos::ParameterList ParamList ;
351 Amesos_Mumps A_mumps( Problem ) ;
352 ParamList.set( "MaxProcs", -3 );
353 EPETRA_CHK_ERR( A_mumps.SetParameters( ParamList ) );
354 EPETRA_CHK_ERR( A_mumps.SetUseTranspose( transpose ) );
355 EPETRA_CHK_ERR( A_mumps.SymbolicFactorization( ) );
356 EPETRA_CHK_ERR( A_mumps.NumericFactorization( ) );
357 EPETRA_CHK_ERR( A_mumps.Solve( ) );
358
359#endif
360#ifdef HAVE_AMESOS_SUPERLU
361 } else if ( SparseSolver == SUPERLU ) {
362
363 Teuchos::ParameterList ParamList ;
364 Amesos_Superlu A_superlu( Problem ) ;
365 ParamList.set( "MaxProcs", -3 );
366 EPETRA_CHK_ERR( A_superlu.SetParameters( ParamList ) );
367 EPETRA_CHK_ERR( A_superlu.SetUseTranspose( transpose ) );
368 EPETRA_CHK_ERR( A_superlu.SymbolicFactorization( ) );
369 EPETRA_CHK_ERR( A_superlu.NumericFactorization( ) );
370 EPETRA_CHK_ERR( A_superlu.Solve( ) );
371
372#endif
373#ifdef HAVE_AMESOS_LAPACK
374 } else if ( SparseSolver == LAPACK ) {
375
376 Teuchos::ParameterList ParamList ;
377 Amesos_Lapack A_lapack( Problem ) ;
378 ParamList.set( "MaxProcs", -3 );
379 EPETRA_CHK_ERR( A_lapack.SetParameters( ParamList ) );
380 EPETRA_CHK_ERR( A_lapack.SetUseTranspose( transpose ) );
381 EPETRA_CHK_ERR( A_lapack.SymbolicFactorization( ) );
382 EPETRA_CHK_ERR( A_lapack.NumericFactorization( ) );
383 EPETRA_CHK_ERR( A_lapack.Solve( ) );
384#endif
385#ifdef HAVE_AMESOS_UMFPACK
386 } else if ( SparseSolver == UMFPACK ) {
387
388 Teuchos::ParameterList ParamList ;
389 Amesos_Umfpack A_umfpack( Problem ) ;
390 ParamList.set( "MaxProcs", -3 );
391 EPETRA_CHK_ERR( A_umfpack.SetParameters( ParamList ) );
392 EPETRA_CHK_ERR( A_umfpack.SetUseTranspose( transpose ) );
393 EPETRA_CHK_ERR( A_umfpack.SymbolicFactorization( ) );
394 EPETRA_CHK_ERR( A_umfpack.NumericFactorization( ) );
395 EPETRA_CHK_ERR( A_umfpack.Solve( ) );
396#endif
397#ifdef HAVE_AMESOS_KLU
398 } else if ( SparseSolver == KLU ) {
399
400
401 using namespace Teuchos;
402
403 Amesos_Time AT;
404 int setupTimePtr = -1, symTimePtr = -1, numTimePtr = -1, refacTimePtr = -1, solveTimePtr = -1;
405 AT.CreateTimer(Comm, 2);
406 AT.ResetTimer(0);
407
408 Teuchos::ParameterList ParamList ;
409 // ParamList.set("OutputLevel",2);
410 Amesos_Klu A_klu( Problem );
411 ParamList.set( "MaxProcs", -3 );
412 ParamList.set( "TrustMe", false );
413 // ParamList.set( "Refactorize", true );
414 EPETRA_CHK_ERR( A_klu.SetParameters( ParamList ) ) ;
415 EPETRA_CHK_ERR( A_klu.SetUseTranspose( transpose ) );
416 setupTimePtr = AT.AddTime("Setup", setupTimePtr, 0);
417 EPETRA_CHK_ERR( A_klu.SymbolicFactorization( ) );
418 symTimePtr = AT.AddTime("Symbolic", symTimePtr, 0);
419 EPETRA_CHK_ERR( A_klu.NumericFactorization( ) );
420 numTimePtr = AT.AddTime("Numeric", numTimePtr, 0);
421 EPETRA_CHK_ERR( A_klu.NumericFactorization( ) );
422 refacTimePtr = AT.AddTime("Refactor", refacTimePtr, 0);
423 // for ( int i=0; i<100000 ; i++ )
424 EPETRA_CHK_ERR( A_klu.Solve( ) );
425 solveTimePtr = AT.AddTime("Solve", solveTimePtr, 0);
426
427 double SetupTime = AT.GetTime(setupTimePtr);
428 double SymbolicTime = AT.GetTime(symTimePtr);
429 double NumericTime = AT.GetTime(numTimePtr);
430 double RefactorTime = AT.GetTime(refacTimePtr);
431 double SolveTime = AT.GetTime(solveTimePtr);
432
433 std::cout << __FILE__ << "::" << __LINE__ << " SetupTime = " << SetupTime << std::endl ;
434 std::cout << __FILE__ << "::" << __LINE__ << " SymbolicTime = " << SymbolicTime - SetupTime << std::endl ;
435 std::cout << __FILE__ << "::" << __LINE__ << " NumericTime = " << NumericTime - SymbolicTime<< std::endl ;
436 std::cout << __FILE__ << "::" << __LINE__ << " RefactorTime = " << RefactorTime - NumericTime << std::endl ;
437 std::cout << __FILE__ << "::" << __LINE__ << " SolveTime = " << SolveTime - RefactorTime << std::endl ;
438
439#endif
440 } else {
441 SparseDirectTimingVars::log_file << "Solver not implemented yet" << std::endl ;
442 std::cerr << "\n\n#################### Requested solver not available on this platform ##################### ATS\n" << std::endl ;
443 std::cout << " SparseSolver = " << SparseSolver << std::endl ;
444 std::cerr << " SparseSolver = " << SparseSolver << std::endl ;
445 }
446
447 SparseDirectTimingVars::SS_Result.Set_Total_Time( TotalTime.ElapsedTime() );
448 } // end for (int i=0; i<special; i++ )
449
450 //
451 // Compute the error = norm(xcomp - xexact )
452 //
453 double error;
454 passresid->Update(1.0, *passx, -1.0, *passxexact, 0.0);
455
456 passresid->Norm2(&error);
457 SparseDirectTimingVars::SS_Result.Set_Error(error) ;
458
459 // passxexact->Norm2(&error ) ;
460 // passx->Norm2(&error ) ;
461
462 //
463 // Compute the residual = norm(Ax - b)
464 //
465 double residual ;
466
467 passA->Multiply( transpose, *passx, *passtmp);
468 passresid->Update(1.0, *passtmp, -1.0, *passb, 0.0);
469 // passresid->Update(1.0, *passtmp, -1.0, CopyB, 0.0);
470 passresid->Norm2(&residual);
471
472 SparseDirectTimingVars::SS_Result.Set_Residual(residual) ;
473
474 double bnorm;
475 passb->Norm2( &bnorm ) ;
476 SparseDirectTimingVars::SS_Result.Set_Bnorm(bnorm) ;
477
478 double xnorm;
479 passx->Norm2( &xnorm ) ;
480 SparseDirectTimingVars::SS_Result.Set_Xnorm(xnorm) ;
481
482 delete readA;
483 delete readx;
484 delete readb;
485 delete readxexact;
486 delete readMap;
487 delete map_;
488
489 Comm.Barrier();
490
491 return 0;
492}
int Amesos_TestSolver(Epetra_Comm &Comm, char *matrix_file, SparseSolverType SparseSolver, bool transpose, int special, AMESOS_MatrixType matrix_type)
AMESOS_MatrixType
@ AMESOS_Distributed
SparseSolverType
@ TAUCS
@ PARAKLETE
@ SCALAPACK
@ LAPACK
@ UMFPACK
@ SUPERLUDIST
@ MUMPS
@ DSCPACK
@ SUPERLU
@ PARDISO
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
#define EPETRA_CHK_ERR(xxx)
Amesos_Dscpack: An object-oriented wrapper for Dscpack.
Amesos_Klu: A serial, unblocked code ideal for getting started and for very sparse matrices,...
Definition Amesos_Klu.h:116
Amesos_Lapack: an interface to LAPACK.
Amesos_Mumps: An object-oriented wrapper for the double precision version of MUMPS.
Amesos_Paraklete: A serial, unblocked code ideal for getting started and for very sparse matrices,...
Amesos_Pardiso: Interface to the PARDISO package.
Amesos_Scalapack: A serial and parallel dense solver. For now, we implement only the unsymmetric ScaL...
Amesos_Superlu: Amesos interface to Xioye Li's SuperLU 3.0 serial code.
Amesos_Superludist: An object-oriented wrapper for Superludist.
Amesos_Taucs: An interface to the TAUCS package.
Amesos_Time: Container for timing information.
Definition Amesos_Time.h:51
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
Definition Amesos_Time.h:64
Class Amesos_Umfpack: An object-oriented wrapper for UMFPACK.
static std::ofstream log_file
static SparseSolverResult SS_Result