Amesos Package Browser (Single Doxygen Collection) Development
Amesos_TestMultiSolver.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                Amesos: Direct Sparse Solver Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Amesos_ConfigDefs.h"
00030 #include "Teuchos_ParameterList.hpp"
00031 //#include "Trilinos_Util_ReadTriples2Epetra.h"
00032 //#include "Trilinos_Util_ReadMatrixMarket2Epetra.h"
00033 #include "Trilinos_Util.h"
00034 #include "Epetra_LocalMap.h"
00035 #include "Epetra_Map.h"
00036 #include "Epetra_Vector.h"
00037 #include "Epetra_Export.h"
00038 #include "Epetra_CrsMatrix.h"
00039 #include "Epetra_LinearProblem.h"
00040 #include "Epetra_Time.h"
00041 #ifdef HAVE_AMESOS_SLUD
00042 #include "SuperludistOO.h"
00043 #endif
00044 #ifdef HAVE_AMESOS_SLUS
00045 #include "Epetra_SLU.h"
00046 #endif
00047 #ifdef HAVE_AMESOS_SLUD2
00048 #include "Superludist2_OO.h"
00049 #endif
00050 #ifdef HAVE_AMESOS_DSCPACK
00051 #include "Amesos_Dscpack.h"
00052 #endif
00053 #ifdef HAVE_AMESOS_UMFPACK
00054 #include "Amesos_Umfpack.h"
00055 #endif
00056 #ifdef HAVE_AMESOS_KLU
00057 #include "Amesos_Klu.h"
00058 #endif
00059 #ifdef HAVE_AMESOS_LAPACK
00060 #include "Amesos_Lapack.h"
00061 #endif
00062 #ifdef HAVE_AMESOS_TAUCS
00063 #include "Amesos_Taucs.h"
00064 #endif
00065 #ifdef HAVE_AMESOS_PARDISO
00066 #include "Amesos_Pardiso.h"
00067 #endif
00068 #ifdef HAVE_AMESOS_PARAKLETE
00069 #include "Amesos_Paraklete.h"
00070 #endif
00071 #ifdef HAVE_AMESOS_MUMPS
00072 #include "Amesos_Mumps.h"
00073 #endif
00074 #ifdef HAVE_AMESOS_SCALAPACK
00075 #include "Amesos_Scalapack.h"
00076 #endif
00077 #ifdef HAVE_AMESOS_SUPERLUDIST
00078 #include "Amesos_Superludist.h"
00079 #endif
00080 #ifdef HAVE_AMESOS_SUPERLU
00081 #include "Amesos_Superlu.h"
00082 #endif
00083 #ifdef TEST_SPOOLES
00084 #include "SpoolesOO.h"
00085 #endif
00086 
00087 #include "Amesos_TestSolver.h"
00088 #include "CrsMatrixTranspose.h"
00089 #include "SparseDirectTimingVars.h"
00090 
00091 //
00092 //  Amesos_TestMultiSolver.cpp reads in a matrix in Harwell-Boeing format, 
00093 //  calls one of the sparse direct solvers, using blocked right hand sides
00094 //  and computes the error and residual.  
00095 //
00096 //  TestSolver ignores the Harwell-Boeing right hand sides, creating
00097 //  random right hand sides instead.  
00098 //
00099 //  Amesos_TestMultiSolver can test either A x = b or A^T x = b.
00100 //  This can be a bit confusing because sparse direct solvers 
00101 //  use compressed column storage - the transpose of Trilinos'
00102 //  sparse row storage.
00103 //
00104 //  Matrices:
00105 //    readA - Serial.  As read from the file.
00106 //    transposeA - Serial.  The transpose of readA.
00107 //    serialA - if (transpose) then transposeA else readA 
00108 //    distributedA - readA distributed to all processes
00109 //    passA - if ( distributed ) then distributedA else serialA
00110 //
00111 //
00112 int Amesos_TestMultiSolver( Epetra_Comm &Comm, char *matrix_file, int numsolves, 
00113           SparseSolverType SparseSolver, bool transpose,
00114           int special, AMESOS_MatrixType matrix_type ) {
00115 
00116 
00117   int iam = Comm.MyPID() ;
00118 
00119   
00120   //  int hatever;
00121   //  if ( iam == 0 )  std::cin >> hatever ; 
00122   Comm.Barrier();
00123 
00124 
00125   Epetra_Map * readMap;
00126   Epetra_CrsMatrix * readA; 
00127   Epetra_Vector * readx; 
00128   Epetra_Vector * readb;
00129   Epetra_Vector * readxexact;
00130    
00131   std::string FileName = matrix_file ;
00132   int FN_Size = FileName.size() ; 
00133   std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
00134   std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
00135   bool NonContiguousMap = false; 
00136 
00137   if ( LastFiveBytes == ".triU" ) { 
00138     NonContiguousMap = true; 
00139     // Call routine to read in unsymmetric Triplet matrix
00140     EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, false, Comm, readMap, readA, readx, 
00141                   readb, readxexact, NonContiguousMap ) );
00142   } else {
00143     if ( LastFiveBytes == ".triS" ) { 
00144       NonContiguousMap = true; 
00145       // Call routine to read in symmetric Triplet matrix
00146       EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, true, Comm, 
00147               readMap, readA, readx, 
00148               readb, readxexact, NonContiguousMap ) );
00149     } else {
00150       if (  LastFourBytes == ".mtx" ) { 
00151   EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( matrix_file, Comm, readMap, 
00152                      readA, readx, readb, readxexact) );
00153       } else {
00154   // Call routine to read in HB problem
00155   Trilinos_Util_ReadHb2Epetra( matrix_file, Comm, readMap, readA, readx, 
00156                  readb, readxexact) ;
00157       }
00158     }
00159   }
00160 
00161   Epetra_CrsMatrix transposeA(Copy, *readMap, 0);
00162   Epetra_CrsMatrix *serialA ; 
00163 
00164   if ( transpose ) {
00165     assert( CrsMatrixTranspose( readA, &transposeA ) == 0 ); 
00166     serialA = &transposeA ; 
00167   } else {
00168     serialA = readA ; 
00169   }
00170 
00171   // Create uniform distributed map
00172   Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
00173   Epetra_Map* map_;
00174 
00175   if( NonContiguousMap ) {
00176     //
00177     //  map gives us NumMyElements and MyFirstElement;
00178     //
00179     int NumGlobalElements =  readMap->NumGlobalElements();
00180     int NumMyElements = map.NumMyElements();
00181     int MyFirstElement = map.MinMyGID();
00182     std::vector<int> MapMap_( NumGlobalElements );
00183     readMap->MyGlobalElements( &MapMap_[0] ) ;
00184     Comm.Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ; 
00185     map_ = new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);
00186   } else {
00187     map_ = new Epetra_Map( map ) ; 
00188   }
00189 
00190 
00191   // Create Exporter to distribute read-in matrix and vectors
00192   Epetra_Export exporter(*readMap, *map_);
00193   Epetra_CrsMatrix A(Copy, *map_, 0);
00194 
00195   Epetra_RowMatrix * passA = 0; 
00196   Epetra_MultiVector * passx = 0; 
00197   Epetra_MultiVector * passb = 0;
00198   Epetra_MultiVector * passxexact = 0;
00199   Epetra_MultiVector * passresid = 0;
00200   Epetra_MultiVector * passtmp = 0;
00201 
00202   Epetra_MultiVector x(*map_,numsolves);
00203   Epetra_MultiVector b(*map_,numsolves);
00204   Epetra_MultiVector xexact(*map_,numsolves);
00205   Epetra_MultiVector resid(*map_,numsolves);
00206   Epetra_MultiVector tmp(*map_,numsolves);
00207 
00208   Epetra_MultiVector serialx(*readMap,numsolves);
00209   Epetra_MultiVector serialb(*readMap,numsolves);
00210   Epetra_MultiVector serialxexact(*readMap,numsolves);
00211   Epetra_MultiVector serialresid(*readMap,numsolves);
00212   Epetra_MultiVector serialtmp(*readMap,numsolves);
00213 
00214   bool distribute_matrix = ( matrix_type == AMESOS_Distributed ) ; 
00215   if ( distribute_matrix ) { 
00216     //
00217     //  Initialize x, b and xexact to the values read in from the file
00218     //
00219     
00220     A.Export(*serialA, exporter, Add);
00221     Comm.Barrier();
00222 
00223     assert(A.FillComplete()==0);    
00224     Comm.Barrier();
00225 
00226     passA = &A; 
00227     passx = &x; 
00228     passb = &b;
00229     passxexact = &xexact;
00230     passresid = &resid;
00231     passtmp = &tmp;
00232   } else { 
00233     passA = serialA; 
00234     passx = &serialx; 
00235     passb = &serialb;
00236     passxexact = &serialxexact;
00237     passresid = &serialresid;
00238     passtmp = &serialtmp;
00239   }
00240 
00241   passxexact->SetSeed(131) ; 
00242   passxexact->Random();
00243   passx->SetSeed(11231) ; 
00244   passx->Random();
00245 
00246   passb->PutScalar( 0.0 );
00247   passA->Multiply( transpose, *passxexact, *passb ) ; 
00248 
00249   Epetra_MultiVector CopyB( *passb ) ;
00250 
00251   double Anorm = passA->NormInf() ; 
00252   SparseDirectTimingVars::SS_Result.Set_Anorm(Anorm) ;
00253 
00254   Epetra_LinearProblem Problem(  (Epetra_RowMatrix *) passA, 
00255          (Epetra_MultiVector *) passx, 
00256          (Epetra_MultiVector *) passb );
00257 
00258   double max_resid = 0.0;
00259   for ( int j = 0 ; j < special+1 ; j++ ) { 
00260     
00261     Epetra_Time TotalTime( Comm ) ; 
00262     if ( false ) { 
00263 #ifdef TEST_UMFPACK
00264 
00265       unused code
00266 
00267     } else if ( SparseSolver == UMFPACK ) { 
00268       UmfpackOO umfpack( (Epetra_RowMatrix *) passA, 
00269        (Epetra_MultiVector *) passx, 
00270        (Epetra_MultiVector *) passb ) ; 
00271     
00272       umfpack.SetTrans( transpose ) ; 
00273       umfpack.Solve() ; 
00274 #endif
00275 #ifdef TEST_SUPERLU
00276     } else if ( SparseSolver == SuperLU ) { 
00277       SuperluserialOO superluserial( (Epetra_RowMatrix *) passA, 
00278              (Epetra_MultiVector *) passx, 
00279              (Epetra_MultiVector *) passb ) ; 
00280 
00281       superluserial.SetPermc( SuperLU_permc ) ; 
00282       superluserial.SetTrans( transpose ) ; 
00283       superluserial.SetUseDGSSV( special == 0 ) ; 
00284       superluserial.Solve() ; 
00285 #endif
00286 #ifdef HAVE_AMESOS_SLUD
00287     } else if ( SparseSolver == SuperLUdist ) { 
00288       SuperludistOO superludist( Problem ) ; 
00289       superludist.SetTrans( transpose ) ; 
00290       EPETRA_CHK_ERR( superludist.Solve( true ) ) ;
00291 #endif 
00292 #ifdef HAVE_AMESOS_SLUD2
00293     } else if ( SparseSolver == SuperLUdist2 ) { 
00294       Superludist2_OO superludist2( Problem ) ; 
00295       superludist2.SetTrans( transpose ) ; 
00296       EPETRA_CHK_ERR( superludist2.Solve( true ) ) ;
00297 #endif 
00298 #ifdef TEST_SPOOLES
00299     } else if ( SparseSolver == SPOOLES ) { 
00300       SpoolesOO spooles( (Epetra_RowMatrix *) passA, 
00301        (Epetra_MultiVector *) passx, 
00302        (Epetra_MultiVector *) passb ) ; 
00303     
00304       spooles.SetTrans( transpose ) ; 
00305       spooles.Solve() ; 
00306 #endif
00307 #ifdef HAVE_AMESOS_DSCPACK
00308     } else if ( SparseSolver == DSCPACK ) { 
00309       Teuchos::ParameterList ParamList ;
00310       Amesos_Dscpack dscpack( Problem ) ; 
00311       ParamList.set( "MaxProcs", -3 );
00312       EPETRA_CHK_ERR( dscpack.SetParameters( ParamList ) ); 
00313     
00314       EPETRA_CHK_ERR( dscpack.Solve( ) ); 
00315 #endif
00316 #ifdef HAVE_AMESOS_UMFPACK
00317     } else if ( SparseSolver == UMFPACK ) { 
00318       Teuchos::ParameterList ParamList ;
00319       Amesos_Umfpack umfpack( Problem ) ; 
00320       ParamList.set( "MaxProcs", -3 );
00321       EPETRA_CHK_ERR( umfpack.SetParameters( ParamList ) ); 
00322       EPETRA_CHK_ERR( umfpack.SetUseTranspose( transpose ) ); 
00323     
00324       EPETRA_CHK_ERR( umfpack.Solve( ) ); 
00325 #endif
00326 #ifdef HAVE_AMESOS_KLU
00327     } else if ( SparseSolver == KLU ) { 
00328       Teuchos::ParameterList ParamList ;
00329       Amesos_Klu klu( Problem ) ; 
00330       ParamList.set( "MaxProcs", -3 );
00331       EPETRA_CHK_ERR( klu.SetParameters( ParamList ) ); 
00332       EPETRA_CHK_ERR( klu.SetUseTranspose( transpose ) ); 
00333     
00334       EPETRA_CHK_ERR( klu.SymbolicFactorization(  ) ); 
00335       EPETRA_CHK_ERR( klu.NumericFactorization(  ) ); 
00336       EPETRA_CHK_ERR( klu.Solve( ) ); 
00337 #endif
00338 #ifdef HAVE_AMESOS_PARAKLETE
00339     } else if ( SparseSolver == PARAKLETE ) { 
00340       Teuchos::ParameterList ParamList ;
00341       Amesos_Paraklete paraklete( Problem ) ; 
00342       ParamList.set( "MaxProcs", -3 );
00343       EPETRA_CHK_ERR( paraklete.SetParameters( ParamList ) ); 
00344       EPETRA_CHK_ERR( paraklete.SetUseTranspose( transpose ) ); 
00345     
00346       EPETRA_CHK_ERR( paraklete.SymbolicFactorization(  ) ); 
00347       EPETRA_CHK_ERR( paraklete.NumericFactorization(  ) ); 
00348       EPETRA_CHK_ERR( paraklete.Solve( ) ); 
00349 #endif
00350 #ifdef HAVE_AMESOS_SLUS
00351     } else if ( SparseSolver == SuperLU ) { 
00352       Epetra_SLU superluserial( &Problem ) ; 
00353       EPETRA_CHK_ERR( superluserial.SetUseTranspose( transpose ) ); 
00354     
00355       EPETRA_CHK_ERR( superluserial.SymbolicFactorization(  ) ); 
00356       EPETRA_CHK_ERR( superluserial.NumericFactorization(  ) ); 
00357 
00358       EPETRA_CHK_ERR( superluserial.Solve( ) ); 
00359 #endif
00360 #ifdef HAVE_AMESOS_LAPACK
00361     } else if ( SparseSolver == LAPACK ) { 
00362       Teuchos::ParameterList ParamList ;
00363       ParamList.set( "MaxProcs", -3 );
00364       Amesos_Lapack lapack( Problem ) ; 
00365       EPETRA_CHK_ERR( lapack.SetUseTranspose( transpose ) ); 
00366     
00367       EPETRA_CHK_ERR( lapack.SymbolicFactorization(  ) ); 
00368       EPETRA_CHK_ERR( lapack.NumericFactorization(  ) ); 
00369       EPETRA_CHK_ERR( lapack.Solve( ) ); 
00370 #endif
00371 #ifdef HAVE_AMESOS_TAUCS
00372     } else if ( SparseSolver == TAUCS ) { 
00373       Teuchos::ParameterList ParamList ;
00374       Amesos_Taucs taucs( Problem ) ; 
00375       ParamList.set( "MaxProcs", -3 );
00376       EPETRA_CHK_ERR( taucs.SetParameters( ParamList ) ); 
00377       EPETRA_CHK_ERR( taucs.SetUseTranspose( transpose ) ); 
00378     
00379       EPETRA_CHK_ERR( taucs.SymbolicFactorization( ) ); 
00380       EPETRA_CHK_ERR( taucs.NumericFactorization( ) ); 
00381       EPETRA_CHK_ERR( taucs.Solve( ) ); 
00382 #endif
00383 #ifdef HAVE_AMESOS_PARDISO
00384     } else if ( SparseSolver == PARDISO ) { 
00385       Teuchos::ParameterList ParamList ;
00386       Amesos_Pardiso pardiso( Problem ) ; 
00387       ParamList.set( "MaxProcs", -3 );
00388       EPETRA_CHK_ERR( pardiso.SetParameters( ParamList ) ); 
00389       EPETRA_CHK_ERR( pardiso.SetUseTranspose( transpose ) ); 
00390     
00391       EPETRA_CHK_ERR( pardiso.SymbolicFactorization( ) ); 
00392       EPETRA_CHK_ERR( pardiso.NumericFactorization( ) ); 
00393       EPETRA_CHK_ERR( pardiso.Solve( ) ); 
00394 #endif
00395 #ifdef HAVE_AMESOS_PARKLETE
00396     } else if ( SparseSolver == PARKLETE ) { 
00397       Teuchos::ParameterList ParamList ;
00398       Amesos_Parklete parklete( Problem ) ; 
00399       ParamList.set( "MaxProcs", -3 );
00400       EPETRA_CHK_ERR( parklete.SetParameters( ParamList ) ); 
00401       EPETRA_CHK_ERR( parklete.SetUseTranspose( transpose ) ); 
00402     
00403       EPETRA_CHK_ERR( parklete.SymbolicFactorization( ) ); 
00404       EPETRA_CHK_ERR( parklete.NumericFactorization( ) ); 
00405       EPETRA_CHK_ERR( parklete.Solve( ) ); 
00406 #endif
00407 #ifdef HAVE_AMESOS_MUMPS
00408     } else if ( SparseSolver == MUMPS ) { 
00409       Teuchos::ParameterList ParamList ;
00410       Amesos_Mumps mumps( Problem ) ; 
00411       ParamList.set( "MaxProcs", -3 );
00412       EPETRA_CHK_ERR( mumps.SetParameters( ParamList ) ); 
00413       EPETRA_CHK_ERR( mumps.SetUseTranspose( transpose ) ); 
00414     
00415       EPETRA_CHK_ERR( mumps.SymbolicFactorization( ) ); 
00416       EPETRA_CHK_ERR( mumps.NumericFactorization( ) ); 
00417       EPETRA_CHK_ERR( mumps.Solve( ) ); 
00418 #endif
00419 #ifdef HAVE_AMESOS_SCALAPACK
00420     } else if ( SparseSolver == SCALAPACK ) { 
00421       Teuchos::ParameterList ParamList ;
00422       Amesos_Scalapack scalapack( Problem ) ; 
00423       ParamList.set( "MaxProcs", -3 );
00424       EPETRA_CHK_ERR( scalapack.SetParameters( ParamList ) ); 
00425       EPETRA_CHK_ERR( scalapack.SetUseTranspose( transpose ) ); 
00426     
00427       EPETRA_CHK_ERR( scalapack.SymbolicFactorization( ) ); 
00428       EPETRA_CHK_ERR( scalapack.NumericFactorization( ) ); 
00429       EPETRA_CHK_ERR( scalapack.Solve( ) ); 
00430 #endif
00431 #ifdef HAVE_AMESOS_SUPERLUDIST
00432     } else if ( SparseSolver == SUPERLUDIST ) { 
00433       Teuchos::ParameterList ParamList ;
00434       Amesos_Superludist superludist( Problem ) ; 
00435       ParamList.set( "MaxProcs", -3 );
00436       EPETRA_CHK_ERR( superludist.SetParameters( ParamList ) ); 
00437 
00438       EPETRA_CHK_ERR( superludist.SetUseTranspose( transpose ) ); 
00439     
00440       EPETRA_CHK_ERR( superludist.SymbolicFactorization(  ) ); 
00441       EPETRA_CHK_ERR( superludist.NumericFactorization(  ) ); 
00442       EPETRA_CHK_ERR( superludist.Solve( ) ); 
00443 #endif
00444 #ifdef HAVE_AMESOS_SUPERLU
00445     } else if ( SparseSolver == SUPERLU ) { 
00446       Teuchos::ParameterList ParamList ;
00447       Amesos_Superlu superlu( Problem ) ; 
00448       ParamList.set( "MaxProcs", -3 );
00449       EPETRA_CHK_ERR( superlu.SetParameters( ParamList ) ); 
00450       EPETRA_CHK_ERR( superlu.SetUseTranspose( transpose ) ); 
00451     
00452       EPETRA_CHK_ERR( superlu.SymbolicFactorization(  ) ); 
00453       EPETRA_CHK_ERR( superlu.NumericFactorization(  ) ); 
00454       EPETRA_CHK_ERR( superlu.Solve( ) ); 
00455 #endif
00456 #ifdef TEST_SPOOLESSERIAL 
00457     } else if ( SparseSolver == SPOOLESSERIAL ) { 
00458       SpoolesserialOO spoolesserial( (Epetra_RowMatrix *) passA, 
00459              (Epetra_MultiVector *) passx, 
00460              (Epetra_MultiVector *) passb ) ; 
00461     
00462       spoolesserial.Solve() ;
00463 #endif
00464     } else { 
00465       SparseDirectTimingVars::log_file << "Solver not implemented yet" << std::endl ;
00466       std::cerr << "\n\n####################  Requested solver not available (Or not tested with blocked RHS) on this platform #####################\n" << std::endl ;
00467     }
00468 
00469     SparseDirectTimingVars::SS_Result.Set_Total_Time( TotalTime.ElapsedTime() ); 
00470     //    SparseDirectTimingVars::SS_Result.Set_First_Time( 0.0 ); 
00471     //    SparseDirectTimingVars::SS_Result.Set_Middle_Time( 0.0 ); 
00472     //    SparseDirectTimingVars::SS_Result.Set_Last_Time( 0.0 ); 
00473 
00474     //
00475     //  Compute the error = norm(xcomp - xexact )
00476     //
00477     std::vector <double> error(numsolves) ; 
00478     double max_error = 0.0;
00479   
00480     passresid->Update(1.0, *passx, -1.0, *passxexact, 0.0);
00481 
00482     passresid->Norm2(&error[0]);
00483     for ( int i = 0 ; i< numsolves; i++ ) 
00484       if ( error[i] > max_error ) max_error = error[i] ; 
00485     SparseDirectTimingVars::SS_Result.Set_Error(max_error) ;
00486 
00487     //  passxexact->Norm2(&error[0] ) ; 
00488     //  passx->Norm2(&error ) ; 
00489 
00490     //
00491     //  Compute the residual = norm(Ax - b)
00492     //
00493     std::vector <double> residual(numsolves) ; 
00494   
00495     passtmp->PutScalar(0.0);
00496     passA->Multiply( transpose, *passx, *passtmp);
00497     passresid->Update(1.0, *passtmp, -1.0, *passb, 0.0); 
00498     //    passresid->Update(1.0, *passtmp, -1.0, CopyB, 0.0); 
00499     passresid->Norm2(&residual[0]);
00500 
00501     for ( int i = 0 ; i< numsolves; i++ ) 
00502       if ( residual[i] > max_resid ) max_resid = residual[i] ; 
00503 
00504 
00505     SparseDirectTimingVars::SS_Result.Set_Residual(max_resid) ;
00506     
00507     std::vector <double> bnorm(numsolves); 
00508     passb->Norm2( &bnorm[0] ) ; 
00509     SparseDirectTimingVars::SS_Result.Set_Bnorm(bnorm[0]) ;
00510 
00511     std::vector <double> xnorm(numsolves); 
00512     passx->Norm2( &xnorm[0] ) ; 
00513     SparseDirectTimingVars::SS_Result.Set_Xnorm(xnorm[0]) ;
00514 
00515 
00516     if ( false && iam == 0 ) { 
00517 
00518       std::cout << " Amesos_TestMutliSolver.cpp " << std::endl ; 
00519       for ( int i = 0 ; i< numsolves && i < 10 ; i++ ) {
00520   std::cout << "i=" << i 
00521        << " error = " << error[i] 
00522        << " xnorm = " << xnorm[i] 
00523        << " residual = " << residual[i] 
00524        << " bnorm = " << bnorm[i] 
00525        << std::endl ; 
00526       
00527       }
00528     
00529       std::cout << std::endl << " max_resid = " << max_resid ; 
00530       std::cout << " max_error = " << max_error << std::endl ; 
00531       std::cout << " Get_residual() again = " << SparseDirectTimingVars::SS_Result.Get_Residual() << std::endl ;
00532 
00533     }
00534   }
00535   delete readA;
00536   delete readx;
00537   delete readb;
00538   delete readxexact;
00539   delete readMap;
00540   delete map_;
00541   
00542   Comm.Barrier();
00543 
00544 return 0 ;
00545 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines