Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Amesos_TestMrhsSolver.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 "Trilinos_Util_ReadTriples2Epetra.h"
32//#include "Trilinos_Util_ReadMatrixMarket2Epetra.h"
33#include "Trilinos_Util.h"
34#include "Epetra_LocalMap.h"
35#include "Epetra_Map.h"
36#include "Epetra_Vector.h"
37#include "Epetra_MultiVector.h"
38#include "Epetra_Export.h"
39#include "Epetra_CrsMatrix.h"
40#include "Epetra_LinearProblem.h"
41#include "Epetra_Time.h"
42#ifdef HAVE_AMESOS_DSCPACK
43#include "Amesos_Dscpack.h"
44#endif
45#ifdef HAVE_AMESOS_LAPACK
46#include "Amesos_Lapack.h"
47#endif
48#ifdef HAVE_AMESOS_SLUS
49#include "Epetra_SLU.h"
50#endif
51#ifdef HAVE_AMESOS_UMFPACK
52#include "Amesos_Umfpack.h"
53#endif
54#ifdef HAVE_AMESOS_KLU
55#include "Amesos_Klu.h"
56#endif
57#ifdef HAVE_AMESOS_TAUCS
58#include "Amesos_Taucs.h"
59#endif
60#ifdef HAVE_AMESOS_PARDISO
61#include "Amesos_Pardiso.h"
62#endif
63#ifdef HAVE_AMESOS_PARAKLETE
64#include "Amesos_Paraklete.h"
65#endif
66#if defined(HAVE_AMESOS_MUMPS) && defined(HAVE_MPI)
67#include "Amesos_Mumps.h"
68#endif
69#ifdef HAVE_AMESOS_SCALAPACK
70#include "Amesos_Scalapack.h"
71#endif
72#ifdef HAVE_AMESOS_SLUD
73#include "SuperludistOO.h"
74#endif
75#ifdef HAVE_AMESOS_SUPERLU
76#include "Amesos_Superlu.h"
77#endif
78#ifdef HAVE_AMESOS_SUPERLUDIST
79#include "Amesos_Superludist.h"
80#endif
81#ifdef HAVE_AMESOS_SLUD2
82#include "Superludist2_OO.h"
83#endif
84#ifdef TEST_SPOOLES
85#include "SpoolesOO.h"
86#endif
87#include "SparseSolverResult.h"
88#include "Amesos_TestSolver.h"
89#include "CrsMatrixTranspose.h"
91
92#include <vector>
93//
94// TestMrhsSolver.cpp reads in a matrix in Harwell-Boeing format,
95// calls one of the sparse direct solvers, using multiple right hand sides
96// (one per solve) and computes the error and residual.
97//
98// TestSolver ignores the Harwell-Boeing right hand sides, creating
99// random right hand sides instead.
100//
101// TestMrhsSolver can test either A x = b or A^T x = b.
102// This can be a bit confusing because sparse direct solvers
103// use compressed column storage - the transpose of Trilinos'
104// sparse row storage.
105//
106// Matrices:
107// readA - Serial. As read from the file.
108// transposeA - Serial. The transpose of readA.
109// serialA - if (transpose) then transposeA else readA
110// distributedA - readA distributed to all processes
111// passA - if ( distributed ) then distributedA else serialA
112//
113//
114int Amesos_TestMrhsSolver( Epetra_Comm &Comm, char *matrix_file, int numsolves,
115 SparseSolverType SparseSolver, bool transpose,
116 int special, AMESOS_MatrixType matrix_type ) {
117
118
119 Comm.Barrier();
120
121 Epetra_Map * readMap;
122 Epetra_CrsMatrix * readA;
123 Epetra_Vector * readx;
124 Epetra_Vector * readb;
125 Epetra_Vector * readxexact;
126
127 std::string FileName = matrix_file ;
128 int FN_Size = FileName.size() ;
129 std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
130 std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
131 bool NonContiguousMap = false;
132
133 if ( LastFiveBytes == ".triU" ) {
134 // Call routine to read in unsymmetric Triplet matrix
135 NonContiguousMap = true;
136 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, false, Comm, readMap, readA, readx,
137 readb, readxexact, NonContiguousMap ) );
138 } else {
139 if ( LastFiveBytes == ".triS" ) {
140 NonContiguousMap = true;
141 // Call routine to read in symmetric Triplet matrix
142 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, true, Comm, readMap, readA, readx,
143 readb, readxexact) );
144 } else {
145 if ( LastFourBytes == ".mtx" ) {
146 EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( matrix_file, Comm, readMap,
147 readA, readx, readb, readxexact) );
148 } else {
149 // Call routine to read in HB problem
150 Trilinos_Util_ReadHb2Epetra( matrix_file, Comm, readMap, readA, readx,
151 readb, readxexact) ;
152 }
153 }
154 }
155
156
157 Epetra_CrsMatrix transposeA(Copy, *readMap, 0);
158 Epetra_CrsMatrix *serialA ;
159
160 if ( transpose ) {
161 assert( CrsMatrixTranspose( readA, &transposeA ) == 0 );
162 serialA = &transposeA ;
163 } else {
164 serialA = readA ;
165 }
166
167
168 // Create uniform distributed map
169 Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
170 Epetra_Map* map_;
171
172 if( NonContiguousMap ) {
173 //
174 // map gives us NumMyElements and MyFirstElement;
175 //
176 int NumGlobalElements = readMap->NumGlobalElements();
177 int NumMyElements = map.NumMyElements();
178 int MyFirstElement = map.MinMyGID();
179 std::vector<int> MapMap_( NumGlobalElements );
180 readMap->MyGlobalElements( &MapMap_[0] ) ;
181 Comm.Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ;
182 map_ = new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);
183 } else {
184 map_ = new Epetra_Map( map ) ;
185 }
186
187
188 // Create Exporter to distribute read-in matrix and vectors
189 Epetra_Export exporter(*readMap, *map_);
190 Epetra_CrsMatrix A(Copy, *map_, 0);
191
192 Epetra_RowMatrix * passA = 0;
193 Epetra_MultiVector * passx = 0;
194 Epetra_MultiVector * passb = 0;
195 Epetra_MultiVector * passxexact = 0;
196 Epetra_MultiVector * passresid = 0;
197 Epetra_MultiVector * passtmp = 0;
198
199 Epetra_MultiVector x(*map_,numsolves);
200 Epetra_MultiVector b(*map_,numsolves);
201 Epetra_MultiVector xexact(*map_,numsolves);
202 Epetra_MultiVector resid(*map_,numsolves);
203 Epetra_MultiVector tmp(*map_,numsolves);
204
205
206 Epetra_MultiVector serialx(*readMap,numsolves);
207 Epetra_MultiVector serialb(*readMap,numsolves);
208 Epetra_MultiVector serialxexact(*readMap,numsolves);
209 Epetra_MultiVector serialresid(*readMap,numsolves);
210 Epetra_MultiVector serialtmp(*readMap,numsolves);
211
212 bool distribute_matrix = ( matrix_type == AMESOS_Distributed ) ;
213 if ( distribute_matrix ) {
214 //
215 // Initialize x, b and xexact to the values read in from the file
216 //
217
218 A.Export(*serialA, exporter, Add);
219 Comm.Barrier();
220
221 assert(A.FillComplete()==0);
222 Comm.Barrier();
223
224 passA = &A;
225 passx = &x;
226 passb = &b;
227 passxexact = &xexact;
228 passresid = &resid;
229 passtmp = &tmp;
230 } else {
231 passA = serialA;
232 passx = &serialx;
233 passb = &serialb;
234 passxexact = &serialxexact;
235 passresid = &serialresid;
236 passtmp = &serialtmp;
237 }
238
239 passxexact->SetSeed(131) ;
240 passxexact->Random();
241 passx->SetSeed(11231) ;
242 passx->Random();
243
244 passb->PutScalar( 0.0 );
245 passA->Multiply( transpose, *passxexact, *passb ) ;
246
247 Epetra_MultiVector CopyB( *passb ) ;
248
249 double Anorm = passA->NormInf() ;
250 SparseDirectTimingVars::SS_Result.Set_Anorm(Anorm) ;
251
252 Epetra_LinearProblem Problem( (Epetra_RowMatrix *) passA,
253 (Epetra_MultiVector *) passx,
254 (Epetra_MultiVector *) passb );
255
256 double max_resid = 0.0;
257 for ( int j = 0 ; j < special+1 ; j++ ) {
258
259 Epetra_Time TotalTime( Comm ) ;
260 if ( false ) {
261#ifdef TEST_UMFPACK
262
263 unused code
264
265 } else if ( SparseSolver == UMFPACK ) {
266 UmfpackOO umfpack( (Epetra_RowMatrix *) passA,
267 (Epetra_MultiVector *) passx,
268 (Epetra_MultiVector *) passb ) ;
269
270 umfpack.SetTrans( transpose ) ;
271 umfpack.Solve() ;
272#endif
273#ifdef TEST_SUPERLU
274 } else if ( SparseSolver == SuperLU ) {
275 SuperluserialOO superluserial ;
276 superluserial.SetUserMatrix( (Epetra_RowMatrix *) passA) ;
277
278 superluserial.SetPermc( SuperLU_permc ) ;
279 superluserial.SetTrans( transpose ) ;
280 superluserial.SetUseDGSSV( special == 0 ) ;
281
282 for ( int i= 0 ; i < numsolves ; i++ ) {
283 // set up to sovle A X[:,i] = B[:,i]
284 Epetra_Vector *passb_i = (*passb)(i) ;
285 Epetra_Vector *passx_i = (*passx)(i) ;
286 superluserial.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
287 superluserial.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
288 // superluserial.SetRHS( (Epetra_MultiVector *) passb_i ;
289 superluserial.Solve() ;
290 if ( i == 0 ) {
291 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
292 } else {
293 if ( i < numsolves-1 )
294 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
295 else
296 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
297 }
298
299 }
300#endif
301#ifdef HAVE_AMESOS_SLUD
302 } else if ( SparseSolver == SuperLUdist ) {
303 SuperludistOO superludist( Problem ) ;
304 superludist.SetTrans( transpose ) ;
305
306 bool factor = true;
307 for ( int i= 0 ; i < numsolves ; i++ ) {
308 // set up to sovle A X[:,i] = B[:,i]
309 Epetra_Vector *passb_i = (*passb)(i) ;
310 Epetra_Vector *passx_i = (*passx)(i) ;
311 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
312 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
313 EPETRA_CHK_ERR( superludist.Solve( factor ) );
314 factor = false;
315 if ( i == 0 )
316 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
317 else {
318 if ( i < numsolves-1 )
319 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
320 else
321 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
322 }
323
324 }
325#endif
326#ifdef HAVE_AMESOS_SLUD2
327 } else if ( SparseSolver == SuperLUdist2 ) {
328 Superludist2_OO superludist2( Problem ) ;
329 superludist2.SetTrans( transpose ) ;
330
331 bool factor = true;
332 for ( int i= 0 ; i < numsolves ; i++ ) {
333 // set up to sovle A X[:,i] = B[:,i]
334 Epetra_Vector *passb_i = (*passb)(i) ;
335 Epetra_Vector *passx_i = (*passx)(i) ;
336 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
337 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
338 EPETRA_CHK_ERR( superludist2.Solve( factor ) );
339 factor = false;
340 if ( i == 0 )
341 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
342 else {
343 if ( i < numsolves-1 )
344 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
345 else
346 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
347 }
348
349 }
350#endif
351#ifdef HAVE_AMESOS_DSCPACK
352 } else if ( SparseSolver == DSCPACK ) {
353 Teuchos::ParameterList ParamList ;
354 Amesos_Dscpack dscpack( Problem ) ;
355 ParamList.set( "MaxProcs", -3 );
356 EPETRA_CHK_ERR( dscpack.SetParameters( ParamList ) );
357
358 for ( int i= 0 ; i < numsolves ; i++ ) {
359 // set up to sovle A X[:,i] = B[:,i]
360 Epetra_Vector *passb_i = (*passb)(i) ;
361 Epetra_Vector *passx_i = (*passx)(i) ;
362 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
363 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
364 EPETRA_CHK_ERR( dscpack.Solve( ) );
365 if ( i == 0 )
366 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
367 else {
368 if ( i < numsolves-1 )
369 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
370 else
371 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
372 }
373
374 }
375#endif
376#ifdef HAVE_AMESOS_UMFPACK
377 } else if ( SparseSolver == UMFPACK ) {
378 Teuchos::ParameterList ParamList ;
379 Amesos_Umfpack umfpack( Problem ) ;
380 ParamList.set( "MaxProcs", -3 );
381 EPETRA_CHK_ERR( umfpack.SetParameters( ParamList ) );
382 EPETRA_CHK_ERR( umfpack.SetUseTranspose( transpose ) );
383
384 for ( int i= 0 ; i < numsolves ; i++ ) {
385 // set up to sovle A X[:,i] = B[:,i]
386 Epetra_Vector *passb_i = (*passb)(i) ;
387 Epetra_Vector *passx_i = (*passx)(i) ;
388 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
389 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
390 EPETRA_CHK_ERR( umfpack.Solve( ) );
391 if ( i == 0 )
392 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
393 else {
394 if ( i < numsolves-1 )
395 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
396 else
397 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
398 }
399
400 }
401#endif
402#ifdef HAVE_AMESOS_SUPERLU
403 } else if ( SparseSolver == SUPERLU ) {
404 Teuchos::ParameterList ParamList ;
405 Amesos_Superlu superlu( Problem ) ;
406 ParamList.set( "MaxProcs", -3 );
407 EPETRA_CHK_ERR( superlu.SetParameters( ParamList ) );
408 EPETRA_CHK_ERR( superlu.SetUseTranspose( transpose ) );
409
410 EPETRA_CHK_ERR( superlu.SymbolicFactorization( ) );
411 EPETRA_CHK_ERR( superlu.NumericFactorization( ) );
412 for ( int i= 0 ; i < numsolves ; i++ ) {
413 // set up to sovle A X[:,i] = B[:,i]
414 Epetra_Vector *passb_i = (*passb)(i) ;
415 Epetra_Vector *passx_i = (*passx)(i) ;
416 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
417 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
418 EPETRA_CHK_ERR( superlu.Solve( ) );
419 if ( i == 0 )
420 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
421 else {
422 if ( i < numsolves-1 )
423 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
424 else
425 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
426 }
427
428 }
429#endif
430#ifdef HAVE_AMESOS_SLUS
431 } else if ( SparseSolver == SuperLU ) {
432 Epetra_SLU superluserial( &Problem ) ;
433
434 bool factor = true;
435
436 for ( int i= 0 ; i < numsolves ; i++ ) {
437 // set up to sovle A X[:,i] = B[:,i]
438 Epetra_Vector *passb_i = (*passb)(i) ;
439 Epetra_Vector *passx_i = (*passx)(i) ;
440 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
441 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
442 EPETRA_CHK_ERR( superluserial.Solve( true, false, factor, 2, -1, true, transpose ) );
443 if ( i == 0 )
444 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
445 else {
446 if ( i < numsolves-1 )
447 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
448 else
449 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
450 }
451
452 }
453#endif
454#ifdef HAVE_AMESOS_KLU
455 } else if ( SparseSolver == KLU ) {
456 Teuchos::ParameterList ParamList ;
457 // ParamList.set("OutputLevel",2);
458 Amesos_Klu klu( Problem ) ;
459 // ParamList.set ("ScaleMethod", 0) ;
460 ParamList.set( "MaxProcs", -3 );
461 EPETRA_CHK_ERR( klu.SetParameters( ParamList ) );
462 ParamList.set( "MaxProcs", -3 );
463 EPETRA_CHK_ERR( klu.SetParameters( ParamList ) );
464 EPETRA_CHK_ERR( klu.SetUseTranspose( transpose ) );
465
466 EPETRA_CHK_ERR( klu.SymbolicFactorization( ) );
467 for ( int trials = 0 ; trials <= 1 ; trials++) {
468 EPETRA_CHK_ERR( klu.NumericFactorization( ) );
469 for ( int i= 0 ; i < numsolves ; i++ ) {
470 // set up to sovle A X[:,i] = B[:,i]
471 Epetra_Vector *passb_i = (*passb)(i) ;
472 Epetra_Vector *passx_i = (*passx)(i) ;
473 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
474 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
475
476 EPETRA_CHK_ERR( klu.Solve( ) );
477 if ( i == 0 ) {
479 TotalTime.ElapsedTime() );
480 } else {
481 if ( i < numsolves-1 )
482 SparseDirectTimingVars::SS_Result.Set_Middle_Time(
483 TotalTime.ElapsedTime() );
484 else
486 TotalTime.ElapsedTime() );
487 }
488 }
489 }
490#endif
491#ifdef HAVE_AMESOS_LAPACK
492 } else if ( SparseSolver == LAPACK ) {
493 Teuchos::ParameterList ParamList ;
494 Amesos_Lapack lapack( Problem ) ;
495 EPETRA_CHK_ERR( lapack.SetUseTranspose( transpose ) );
496
497 EPETRA_CHK_ERR( lapack.SymbolicFactorization( ) );
498 EPETRA_CHK_ERR( lapack.NumericFactorization( ) );
499 for ( int i= 0 ; i < numsolves ; i++ ) {
500 // set up to sovle A X[:,i] = B[:,i]
501 Epetra_Vector *passb_i = (*passb)(i) ;
502 Epetra_Vector *passx_i = (*passx)(i) ;
503 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
504 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
505 EPETRA_CHK_ERR( lapack.Solve( ) );
506 if ( i == 0 )
507 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
508 else {
509 if ( i < numsolves-1 )
510 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
511 else
512 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
513 }
514
515 }
516#endif
517#ifdef HAVE_AMESOS_TAUCS
518 } else if ( SparseSolver == TAUCS ) {
519 Teuchos::ParameterList ParamList ;
520 Amesos_Taucs taucs( Problem ) ;
521 ParamList.set( "MaxProcs", -3 );
522 EPETRA_CHK_ERR( taucs.SetParameters( ParamList ) );
523 EPETRA_CHK_ERR( taucs.SetUseTranspose( transpose ) );
524 EPETRA_CHK_ERR( taucs.SymbolicFactorization( ) );
525 EPETRA_CHK_ERR( taucs.NumericFactorization( ) );
526
527 for ( int i= 0 ; i < numsolves ; i++ ) {
528 // set up to sovle A X[:,i] = B[:,i]
529 Epetra_Vector *passb_i = (*passb)(i) ;
530 Epetra_Vector *passx_i = (*passx)(i) ;
531 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
532 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
533 EPETRA_CHK_ERR( taucs.Solve( ) );
534 if ( i == 0 )
535 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
536 else {
537 if ( i < numsolves-1 )
538 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
539 else
540 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
541 }
542
543 }
544#endif
545#ifdef HAVE_AMESOS_PARDISO
546 } else if ( SparseSolver == PARDISO ) {
547 Teuchos::ParameterList ParamList ;
548 Amesos_Pardiso pardiso( Problem ) ;
549 ParamList.set( "MaxProcs", -3 );
550 EPETRA_CHK_ERR( pardiso.SetParameters( ParamList ) );
551 EPETRA_CHK_ERR( pardiso.SetUseTranspose( transpose ) );
552 EPETRA_CHK_ERR( pardiso.SymbolicFactorization( ) );
553 EPETRA_CHK_ERR( pardiso.NumericFactorization( ) );
554
555 for ( int i= 0 ; i < numsolves ; i++ ) {
556 // set up to sovle A X[:,i] = B[:,i]
557 Epetra_Vector *passb_i = (*passb)(i) ;
558 Epetra_Vector *passx_i = (*passx)(i) ;
559 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
560 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
561 EPETRA_CHK_ERR( pardiso.Solve( ) );
562 if ( i == 0 )
563 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
564 else {
565 if ( i < numsolves-1 )
566 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
567 else
568 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
569 }
570
571 }
572#endif
573#ifdef HAVE_AMESOS_PARAKLETE
574 } else if ( SparseSolver == PARAKLETE ) {
575 Teuchos::ParameterList ParamList ;
576 Amesos_Paraklete paraklete( Problem ) ;
577 ParamList.set( "MaxProcs", -3 );
578 EPETRA_CHK_ERR( paraklete.SetParameters( ParamList ) );
579 EPETRA_CHK_ERR( paraklete.SetUseTranspose( transpose ) );
580 EPETRA_CHK_ERR( paraklete.SymbolicFactorization( ) );
581 EPETRA_CHK_ERR( paraklete.NumericFactorization( ) );
582
583 for ( int i= 0 ; i < numsolves ; i++ ) {
584 // set up to sovle A X[:,i] = B[:,i]
585 Epetra_Vector *passb_i = (*passb)(i) ;
586 Epetra_Vector *passx_i = (*passx)(i) ;
587 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
588 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
589 EPETRA_CHK_ERR( paraklete.Solve( ) );
590 if ( i == 0 )
591 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
592 else {
593 if ( i < numsolves-1 )
594 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
595 else
596 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
597 }
598
599 }
600#endif
601#if defined(HAVE_AMESOS_MUMPS) && defined(HAVE_MPI)
602 } else if ( SparseSolver == MUMPS ) {
603 Teuchos::ParameterList ParamList ;
604 Amesos_Mumps mumps( Problem ) ;
605 ParamList.set( "MaxProcs", -3 );
606 EPETRA_CHK_ERR( mumps.SetParameters( ParamList ) );
607 EPETRA_CHK_ERR( mumps.SetUseTranspose( transpose ) );
608 EPETRA_CHK_ERR( mumps.SymbolicFactorization( ) );
609 EPETRA_CHK_ERR( mumps.NumericFactorization( ) );
610
611 for ( int i= 0 ; i < numsolves ; i++ ) {
612 // set up to sovle A X[:,i] = B[:,i]
613 Epetra_Vector *passb_i = (*passb)(i) ;
614 Epetra_Vector *passx_i = (*passx)(i) ;
615 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
616 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
617 EPETRA_CHK_ERR( mumps.Solve( ) );
618 if ( i == 0 )
619 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
620 else {
621 if ( i < numsolves-1 )
622 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
623 else
624 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
625 }
626
627 }
628#endif
629#ifdef HAVE_AMESOS_SCALAPACK
630 } else if ( SparseSolver == SCALAPACK ) {
631 Teuchos::ParameterList ParamList ;
632 Amesos_Scalapack scalapack( Problem ) ;
633 ParamList.set( "MaxProcs", -3 );
634 EPETRA_CHK_ERR( scalapack.SetParameters( ParamList ) );
635 EPETRA_CHK_ERR( scalapack.SetUseTranspose( transpose ) );
636
637 EPETRA_CHK_ERR( scalapack.SymbolicFactorization( ) );
638 EPETRA_CHK_ERR( scalapack.NumericFactorization( ) );
639 for ( int i= 0 ; i < numsolves ; i++ ) {
640 // set up to sovle A X[:,i] = B[:,i]
641 Epetra_Vector *passb_i = (*passb)(i) ;
642 Epetra_Vector *passx_i = (*passx)(i) ;
643 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
644 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
645 EPETRA_CHK_ERR( scalapack.Solve( ) );
646 if ( i == 0 )
647 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
648 else {
649 if ( i < numsolves-1 )
650 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
651 else
652 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
653 }
654
655 }
656#endif
657#ifdef HAVE_AMESOS_SUPERLUDIST
658 } else if ( SparseSolver == SUPERLUDIST ) {
659 Teuchos::ParameterList ParamList ;
660 ParamList.set( "MaxProcs", -3 );
661 Amesos_Superludist superludist( Problem ) ;
662 EPETRA_CHK_ERR( superludist.SetParameters( ParamList ) );
663 EPETRA_CHK_ERR( superludist.SetUseTranspose( transpose ) );
664 EPETRA_CHK_ERR( superludist.SymbolicFactorization( ) );
665 EPETRA_CHK_ERR( superludist.NumericFactorization( ) );
666 SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() );
667
668 for ( int i= 0 ; i < numsolves ; i++ ) {
669 // set up to sovle A X[:,i] = B[:,i]
670 Epetra_Vector *passb_i = (*passb)(i) ;
671 Epetra_Vector *passx_i = (*passx)(i) ;
672 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
673 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
674 EPETRA_CHK_ERR( superludist.Solve( ) );
675 if ( i < numsolves-1 )
676 SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() );
677 else
678 SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() );
679 }
680#endif
681#ifdef TEST_SPOOLES
682 } else if ( SparseSolver == SPOOLES ) {
683 SpoolesOO spooles( (Epetra_RowMatrix *) passA,
684 (Epetra_MultiVector *) passx,
685 (Epetra_MultiVector *) passb ) ;
686
687 spooles.SetTrans( transpose ) ;
688 spooles.Solve() ;
689#endif
690#ifdef TEST_SPOOLESSERIAL
691 } else if ( SparseSolver == SPOOLESSERIAL ) {
692 SpoolesserialOO spoolesserial( (Epetra_RowMatrix *) passA,
693 (Epetra_MultiVector *) passx,
694 (Epetra_MultiVector *) passb ) ;
695
696 spoolesserial.Solve() ;
697#endif
698 } else {
699 SparseDirectTimingVars::log_file << "Solver not implemented yet" << std::endl ;
700 std::cerr << "\n\n#################### Requested solver not available (Or not tested with multiple RHS) on this platform #####################\n" << std::endl ;
701 }
702
703 SparseDirectTimingVars::SS_Result.Set_Total_Time( TotalTime.ElapsedTime() );
704
705 //
706 // Compute the error = norm(xcomp - xexact )
707 //
708 std::vector <double> error(numsolves) ;
709 double max_error = 0.0;
710
711 passresid->Update(1.0, *passx, -1.0, *passxexact, 0.0);
712
713 passresid->Norm2(&error[0]);
714 for ( int i = 0 ; i< numsolves; i++ )
715 if ( error[i] > max_error ) max_error = error[i] ;
716 SparseDirectTimingVars::SS_Result.Set_Error(max_error) ;
717
718 // passxexact->Norm2(&error[0] ) ;
719 // passx->Norm2(&error ) ;
720
721 //
722 // Compute the residual = norm(Ax - b)
723 //
724 std::vector <double> residual(numsolves) ;
725
726 passtmp->PutScalar(0.0);
727 passA->Multiply( transpose, *passx, *passtmp);
728 passresid->Update(1.0, *passtmp, -1.0, *passb, 0.0);
729 // passresid->Update(1.0, *passtmp, -1.0, CopyB, 0.0);
730 passresid->Norm2(&residual[0]);
731
732 for ( int i = 0 ; i< numsolves; i++ )
733 if ( residual[i] > max_resid ) max_resid = residual[i] ;
734
735
736 SparseDirectTimingVars::SS_Result.Set_Residual(max_resid) ;
737
738 std::vector <double> bnorm(numsolves);
739 passb->Norm2( &bnorm[0] ) ;
740 SparseDirectTimingVars::SS_Result.Set_Bnorm(bnorm[0]) ;
741
742 std::vector <double> xnorm(numsolves);
743 passx->Norm2( &xnorm[0] ) ;
744 SparseDirectTimingVars::SS_Result.Set_Xnorm(xnorm[0]) ;
745
746 }
747 delete readA;
748 delete readx;
749 delete readb;
750 delete readxexact;
751 delete readMap;
752 delete map_;
753
754 Comm.Barrier();
755 return 0;
756}
int SuperLU_permc
int Amesos_TestMrhsSolver(Epetra_Comm &Comm, char *matrix_file, int numsolves, 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.
Class Amesos_Umfpack: An object-oriented wrapper for UMFPACK.
Epetra_SLU: An object-oriented wrapper for Xiaoye Li's serial sparse solver package: Superlu.
Definition Epetra_SLU.h:52
static std::ofstream log_file
static SparseSolverResult SS_Result
SpoolesOO: An object-oriented wrapper for Spooles.
Definition SpoolesOO.h:51
Superludist2_OO: An object-oriented wrapper for Superludist.
SuperludistOO: An object-oriented wrapper for Xiaoye Li's Superludist.