Amesos Package Browser (Single Doxygen Collection) Development
compare_solvers.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                Amesos: An Interface to Direct Solvers
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 // This example needs Galeri to generate the linear system.
00031 #ifdef HAVE_MPI
00032 #include "mpi.h"
00033 #include "Epetra_MpiComm.h"
00034 #else
00035 #include "Epetra_SerialComm.h"
00036 #endif
00037 #include "Epetra_Vector.h"
00038 #include "Epetra_Time.h"
00039 #include "Epetra_RowMatrix.h"
00040 #include "Epetra_CrsMatrix.h"
00041 #include "Amesos.h"
00042 #include "Amesos_BaseSolver.h"
00043 #include "Teuchos_ParameterList.hpp"
00044 #include "Galeri_Maps.h"
00045 #include "Galeri_CrsMatrices.h"
00046 
00047 using namespace Teuchos;
00048 using namespace Galeri;
00049 
00050 #include "Trilinos_Util.h"
00051 #include "Trilinos_Util_ReadMatrixMarket2Epetra.h"
00052 #include "CrsMatrixTranspose.h"
00053 #include "Teuchos_RCP.hpp"
00054 #include "Epetra_Export.h"
00055 
00056 int MyCreateCrsMatrix( char *in_filename, const Epetra_Comm &Comm, 
00057          Epetra_Map *& readMap,
00058          const bool transpose, const bool distribute, 
00059          bool& symmetric, Epetra_CrsMatrix *& Matrix ) {
00060 
00061   Epetra_CrsMatrix * readA = 0; 
00062   Epetra_Vector * readx = 0; 
00063   Epetra_Vector * readb = 0;
00064   Epetra_Vector * readxexact = 0;
00065 
00066   //
00067   //  This hack allows TestOptions to be run from either the test/TestOptions/ directory or from 
00068   //  the test/ directory (as it is in nightly testing and in make "run-tests")
00069   //
00070   FILE *in_file = fopen( in_filename, "r");
00071 
00072   char *filename;
00073   if (in_file == NULL ) 
00074     filename = &in_filename[1] ; //  Strip off ithe "." from
00075          //  "../" and try again 
00076   else {
00077     filename = in_filename ;
00078     fclose( in_file );
00079   }
00080 
00081   symmetric = false ; 
00082   std::string FileName = filename ;
00083 
00084   int FN_Size = FileName.size() ; 
00085   std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
00086   std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
00087 
00088   if ( LastFiveBytes == ".TimD" ) { 
00089     // Call routine to read in a file in a Zero Based File in Tim Davis format 
00090     EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx, 
00091                   readb, readxexact, false, true, true ) );
00092     symmetric = false; 
00093   } else {
00094     if ( LastFiveBytes == ".triU" ) { 
00095       // Call routine to read in unsymmetric Triplet matrix
00096       EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx, 
00097               readb, readxexact) );
00098       symmetric = false; 
00099     } else {
00100       if ( LastFiveBytes == ".triS" ) { 
00101   // Call routine to read in symmetric Triplet matrix
00102   EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, true, Comm, readMap, readA, readx, 
00103                 readb, readxexact) );
00104   symmetric = true; 
00105       } else {
00106   if (  LastFourBytes == ".mtx" ) { 
00107     EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( filename, Comm, readMap, 
00108                  readA, readx, readb, readxexact) );   
00109     FILE* in_file = fopen( filename, "r");
00110     assert (in_file != NULL) ;  // Checked in Trilinos_Util_CountMatrixMarket() 
00111     const int BUFSIZE = 800 ; 
00112     char buffer[BUFSIZE] ; 
00113     fgets( buffer, BUFSIZE, in_file ) ;  // Pick symmetry info off of this string 
00114     std::string headerline1 = buffer;
00115 #ifdef TFLOP
00116     if ( headerline1.find("symmetric") < BUFSIZE ) symmetric = true;
00117 #else
00118     if ( headerline1.find("symmetric") != std::string::npos) symmetric = true; 
00119     
00120 #endif
00121     fclose(in_file);
00122     
00123   } else {
00124     // Call routine to read in HB problem
00125     Trilinos_Util_ReadHb2Epetra( filename, Comm, readMap, readA, readx, 
00126                readb, readxexact) ;
00127     if (  LastFourBytes == ".rsa" ) symmetric = true ; 
00128   }
00129       }
00130     }
00131   }
00132 
00133 
00134   if ( readb )  delete readb;
00135   if ( readx ) delete readx;
00136   if ( readxexact ) delete readxexact;
00137 
00138   Epetra_CrsMatrix *serialA ; 
00139   Epetra_CrsMatrix *transposeA;
00140 
00141   if ( transpose ) {
00142     transposeA = new Epetra_CrsMatrix( Copy, *readMap, 0 );
00143     assert( CrsMatrixTranspose( readA, transposeA ) == 0 ); 
00144     serialA = transposeA ; 
00145     delete readA;
00146     readA = 0 ; 
00147   } else {
00148     serialA = readA ; 
00149   }
00150 
00151   assert( (void *) &serialA->Graph() ) ;
00152   assert( (void *) &serialA->RowMap() ) ;
00153   assert( serialA->RowMap().SameAs(*readMap) ) ; 
00154 
00155   if ( distribute ) { 
00156     // Create uniform distributed map
00157     Epetra_Map DistMap(readMap->NumGlobalElements(), 0, Comm);
00158 
00159     // Create Exporter to distribute read-in matrix and vectors
00160     Epetra_Export exporter( *readMap, DistMap );
00161     
00162     Epetra_CrsMatrix *Amat = new Epetra_CrsMatrix( Copy, DistMap, 0 );
00163     Amat->Export(*serialA, exporter, Add);
00164     assert(Amat->FillComplete()==0);    
00165     
00166     Matrix = Amat; 
00167     //
00168     //  Make sure that deleting Amat->RowMap() will delete map 
00169     //
00170     //  Bug:  We can't manage to delete map his way anyway,
00171     //        and this fails on tranposes, so for now I just accept
00172     //        the memory loss.
00173     //    assert( &(Amat->RowMap()) == map ) ; 
00174     delete readMap; 
00175     readMap = 0 ; 
00176     delete serialA; 
00177   } else { 
00178 
00179     Matrix = serialA; 
00180   }
00181 
00182 
00183   return 0;
00184 }
00185 
00186 
00187 
00188 
00189 
00190 // ===================== //
00191 // M A I N   D R I V E R //
00192 // ===================== //
00193 //
00194 // This example compares all the available Amesos solvers
00195 // for the solution of the same linear system. 
00196 //
00197 // The example can be run in serial and in parallel.
00198 //
00199 // Author: Marzio Sala, SNL 9214
00200 // Last modified: Oct-05.
00201 
00202 int main(int argc, char *argv[]) 
00203 {
00204 #ifdef HAVE_MPI
00205   MPI_Init(&argc, &argv);
00206   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00207 #else
00208   Epetra_SerialComm Comm;
00209 #endif
00210 
00211   bool verbose = (Comm.MyPID() == 0);
00212   double TotalResidual = 0.0;
00213 
00214   // Create the Map, defined as a grid, of size nx x ny x nz,
00215   // subdivided into mx x my x mz cubes, each assigned to a 
00216   // different processor.
00217 
00218 #ifndef FILENAME_SPECIFIED_ON_COMMAND_LINE
00219   ParameterList GaleriList;
00220   GaleriList.set("nx", 4);
00221   GaleriList.set("ny", 4);
00222   GaleriList.set("nz", 4 * Comm.NumProc());
00223   GaleriList.set("mx", 1);
00224   GaleriList.set("my", 1);
00225   GaleriList.set("mz", Comm.NumProc());
00226   Epetra_Map* Map = CreateMap("Cartesian3D", Comm, GaleriList);
00227   
00228   // Create a matrix, in this case corresponding to a 3D Laplacian
00229   // discretized using a classical 7-point stencil. Please refer to 
00230   // the Galeri documentation for an overview of available matrices.
00231   // 
00232   // NOTE: matrix must be symmetric if DSCPACK is used.
00233 
00234   Epetra_CrsMatrix* Matrix = CreateCrsMatrix("Laplace3D", Map, GaleriList);
00235 #else
00236   bool transpose = false ; 
00237   bool distribute = false ; 
00238   bool symmetric ; 
00239   Epetra_CrsMatrix *Matrix = 0 ;
00240   Epetra_Map *Map = 0 ;
00241   MyCreateCrsMatrix( argv[1], Comm, Map, transpose, distribute, symmetric, Matrix ) ;
00242 
00243 #endif
00244 
00245   // build vectors, in this case with 1 vector
00246   Epetra_MultiVector LHS(*Map, 1); 
00247   Epetra_MultiVector RHS(*Map, 1); 
00248     
00249   // create a linear problem object
00250   Epetra_LinearProblem Problem(Matrix, &LHS, &RHS);
00251 
00252   // use this list to set up parameters, now it is required
00253   // to use all the available processes (if supported by the
00254   // underlying solver). Uncomment the following two lines
00255   // to let Amesos print out some timing and status information.
00256   ParameterList List;
00257   List.set("PrintTiming",true);
00258   List.set("PrintStatus",true);
00259   List.set("MaxProcs",Comm.NumProc());
00260 
00261   std::vector<std::string> SolverType;
00262   SolverType.push_back("Amesos_Paraklete");
00263   SolverType.push_back("Amesos_Klu");
00264   Comm.Barrier() ; 
00265 #if 1
00266   SolverType.push_back("Amesos_Lapack");
00267   SolverType.push_back("Amesos_Umfpack");
00268   SolverType.push_back("Amesos_Pardiso");
00269   SolverType.push_back("Amesos_Taucs");
00270   SolverType.push_back("Amesos_Superlu");
00271   SolverType.push_back("Amesos_Superludist");
00272   SolverType.push_back("Amesos_Mumps");
00273   SolverType.push_back("Amesos_Dscpack");
00274   SolverType.push_back("Amesos_Scalapack");
00275 #endif
00276   Epetra_Time Time(Comm);
00277   
00278   // this is the Amesos factory object that will create 
00279   // a specific Amesos solver.
00280   Amesos Factory;
00281 
00282   // Cycle over all solvers.
00283   // Only installed solvers will be tested.
00284   for (unsigned int i = 0 ; i < SolverType.size() ; ++i) 
00285   {
00286     // Check whether the solver is available or not
00287     if (Factory.Query(SolverType[i])) 
00288     {
00289       // 1.- set exact solution (constant vector)
00290       LHS.PutScalar(1.0);
00291  
00292       // 2.- create corresponding rhs
00293       Matrix->Multiply(false, LHS, RHS);
00294  
00295       // 3.- randomize solution vector
00296       LHS.Random();
00297  
00298       // 4.- create the amesos solver object
00299       Amesos_BaseSolver* Solver = Factory.Create(SolverType[i], Problem);
00300       assert (Solver != 0);
00301 
00302       Solver->SetParameters(List);
00303       Solver->SetUseTranspose( true) ; 
00304 
00305       // 5.- factorize and solve
00306       
00307 
00308       Comm.Barrier() ; 
00309       if (verbose) 
00310         std::cout << std::endl
00311              << "Solver " << SolverType[i] 
00312              << ", verbose = " << verbose << std::endl ; 
00313       Comm.Barrier() ; 
00314 
00315 
00316       Time.ResetStartTime();
00317       AMESOS_CHK_ERR(Solver->SymbolicFactorization());
00318       if (verbose) 
00319         std::cout << std::endl
00320              << "Solver " << SolverType[i] 
00321              << ", symbolic factorization time = " 
00322              << Time.ElapsedTime() << std::endl;
00323       Comm.Barrier() ; 
00324 
00325       AMESOS_CHK_ERR(Solver->NumericFactorization());
00326       if (verbose) 
00327         std::cout << "Solver " << SolverType[i] 
00328              << ", numeric factorization time = " 
00329              << Time.ElapsedTime() << std::endl;
00330       Comm.Barrier() ; 
00331 
00332       AMESOS_CHK_ERR(Solver->Solve());
00333       if (verbose) 
00334         std::cout << "Solver " << SolverType[i] 
00335              << ", solve time = " 
00336              << Time.ElapsedTime() << std::endl;
00337       Comm.Barrier() ; 
00338   
00339       // 6.- compute difference between exact solution and Amesos one
00340       //     (there are other ways of doing this in Epetra, but let's
00341       //     keep it simple)
00342       double d = 0.0, d_tot = 0.0;
00343       for (int j = 0 ; j< LHS.Map().NumMyElements() ; ++j)
00344         d += (LHS[0][j] - 1.0) * (LHS[0][j] - 1.0);
00345 
00346       Comm.SumAll(&d,&d_tot,1);
00347       if (verbose)
00348         std::cout << "Solver " << SolverType[i] << ", ||x - x_exact||_2 = " 
00349              << sqrt(d_tot) << std::endl;
00350 
00351       // 7.- delete the object
00352       delete Solver;
00353 
00354       TotalResidual += d_tot;
00355     }
00356   }
00357 
00358   delete Matrix;
00359   delete Map;
00360 
00361   if (TotalResidual > 1e-9) 
00362     exit(EXIT_FAILURE);
00363 
00364 #ifdef HAVE_MPI
00365   MPI_Finalize();
00366 #endif
00367 
00368   return(EXIT_SUCCESS);
00369 } // end of main()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines