Amesos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
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"
90 #include "SparseDirectTimingVars.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 //
114 int 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() ;
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 ) {
292  } else {
293  if ( i < numsolves-1 )
295  else
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 )
317  else {
318  if ( i < numsolves-1 )
320  else
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 )
342  else {
343  if ( i < numsolves-1 )
345  else
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 )
367  else {
368  if ( i < numsolves-1 )
370  else
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 )
393  else {
394  if ( i < numsolves-1 )
396  else
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 )
421  else {
422  if ( i < numsolves-1 )
424  else
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 )
445  else {
446  if ( i < numsolves-1 )
448  else
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 
467  for ( int trials = 0 ; trials <= 1 ; trials++) {
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 )
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 
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 )
508  else {
509  if ( i < numsolves-1 )
511  else
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 ) );
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 )
536  else {
537  if ( i < numsolves-1 )
539  else
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 )
564  else {
565  if ( i < numsolves-1 )
567  else
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 )
592  else {
593  if ( i < numsolves-1 )
595  else
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 ) );
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 )
620  else {
621  if ( i < numsolves-1 )
623  else
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 )
648  else {
649  if ( i < numsolves-1 )
651  else
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( ) );
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 )
677  else
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 
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] ;
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 
737 
738  std::vector <double> bnorm(numsolves);
739  passb->Norm2( &bnorm[0] ) ;
741 
742  std::vector <double> xnorm(numsolves);
743  passx->Norm2( &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 SetUseTranspose(bool UseTranspose)
SetUseTranpose(true) is more efficient in Amesos_Scalapack.
int NumGlobalElements() const
Amesos_Klu: A serial, unblocked code ideal for getting started and for very sparse matrices...
Definition: Amesos_Klu.h:117
void SetLHS(Epetra_MultiVector *X)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetUseTranspose(bool UseTranspose_in)
If set true, X will be set to the solution of AT X = B (not A X = B)
Epetra_SLU: An object-oriented wrapper for Xiaoye Li&#39;s serial sparse solver package: Superlu...
Definition: Epetra_SLU.h:57
int Solve()
Solves A X = B (or AT x = B)
int Solve()
Solves A X = B (or AT x = B)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void SetTrans(bool trans)
Definition: SpoolesOO.h:80
SuperludistOO: An object-oriented wrapper for Xiaoye Li&#39;s Superludist.
Amesos_Paraklete: A serial, unblocked code ideal for getting started and for very sparse matrices...
Amesos_Superludist: An object-oriented wrapper for Superludist.
int Solve()
Solves A X = B (or AT x = B)
Amesos_Mumps: An object-oriented wrapper for the double precision version of MUMPS.
Definition: Amesos_Mumps.h:118
int Solve()
Solves A X = B (or AT x = B)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int MyGlobalElements(int *MyGlobalElementList) const
int Solve()
Solves A X = B (or AT x = B)
Definition: Amesos_Klu.cpp:729
void Set_Total_Time(double Total_Time_in)
double ElapsedTime(void) const
Amesos_Dscpack: An object-oriented wrapper for Dscpack.
int Amesos_TestMrhsSolver(Epetra_Comm &Comm, char *matrix_file, int numsolves, SparseSolverType SparseSolver, bool transpose, int special, AMESOS_MatrixType matrix_type)
void Set_Anorm(double anorm_in)
int SetUseTranspose(bool UseTranspose)
SetUseTranpose()
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
SparseSolverType
void SetTrans(bool trans)
Setting the transpose flag to true causes Solve() to compute A^t x = b.
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
virtual void Barrier() const =0
Amesos_Pardiso: Interface to the PARDISO package.
Amesos_Superlu: Amesos interface to Xioye Li&#39;s SuperLU 3.0 serial code.
int Solve()
Solves A X = B (or AT X = B)
int FillComplete(bool OptimizeDataStorage=true)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void Set_Error(double error_in)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int Solve()
Definition: SpoolesOO.cpp:128
int SetUseTranspose(bool UseTranspose_in)
SetUseTranpose()
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetUseTranspose(bool UseTranspose_in)
If set true, X will be set to the solution of AT X = B (not A X = B)
Definition: Amesos_Mumps.h:151
int SetUseTranspose(bool UseTranspose_in)
If set true, X will be set to the solution of AT X = B (not A X = B)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Set parameters from the input parameters list, returns 0 if successful.
int SuperLU_permc
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
AMESOS_MatrixType
virtual double NormInf() const =0
void Set_Middle_Time(double Middle_Time_in)
int Solve(bool Factor)
All computation is performed during the call to Solve()
int NumericFactorization()
Performs NumericFactorization on the matrix A.
Definition: Amesos_Klu.cpp:681
#define EPETRA_CHK_ERR(xxx)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
std::string error
void Set_Xnorm(double xnorm_in)
int SetUseTranspose(bool useTheTranspose)
If set true, X will be set to the solution of AT X = B (not A X = B)
int Solve()
Solves A X = B (or AT x = B)
virtual int Broadcast(double *MyVals, int Count, int Root) const =0
int SetUseTranspose(bool UseTranspose_in)
SetUseTranpose(true) is more efficient in Amesos_Klu.
Definition: Amesos_Klu.h:166
int SetUseTranspose(bool UseTranspose)
Amesos_Taucs supports only symmetric matrices, hence transpose is irrelevant, but harmless...
Definition: Amesos_Taucs.h:118
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
Amesos_Lapack: an interface to LAPACK.
Definition: Amesos_Lapack.h:71
void Set_Last_Time(double Last_Time_in)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
Definition: Amesos_Klu.cpp:617
void Set_First_Time(double First_Time_in)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
void SetRHS(Epetra_MultiVector *B)
int SetUseTranspose(bool UseTranspose)
Amesos_Superludist does not support transpose at this time.
int Solve()
Solves A X = B (or AT x = B)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int Solve(bool Verbose=false, bool Equil=true, bool Factor=true, int perm_type=2, double pivot_thresh=-1, bool Refact=true, bool Trans=false)
All computation is performed during the call to Solve()
Definition: Epetra_SLU.cpp:202
SpoolesOO: An object-oriented wrapper for Spooles.
Definition: SpoolesOO.h:57
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void Set_Residual(double residual_in)
int Solve(bool Factor)
All computation is performed during the call to Solve()
int Solve()
Solves A X = B (or AT x = B)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
Superludist2_OO: An object-oriented wrapper for Superludist.
static std::ofstream log_file
int NumericFactorization()
Performs NumericFactorization on the matrix A.
static SparseSolverResult SS_Result
int Solve()
Solves A X = B (or AT X = B)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void SetTrans(bool trans)
Setting the transpose flag to true causes Solve() to compute A^t x = b.
Class Amesos_Umfpack: An object-oriented wrapper for UMFPACK.
int Solve()
Solves A X = B (or AT x = B)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
Definition: Amesos_Klu.cpp:449
void Set_Bnorm(double bnorm_in)
Amesos_Scalapack: A serial and parallel dense solver. For now, we implement only the unsymmetric ScaL...
Amesos_Taucs: An interface to the TAUCS package.
Definition: Amesos_Taucs.h:79
#define EPETRA_MAX(x, y)
int NumericFactorization()
Performs NumericFactorization on the matrix A.