Amesos Package Browser (Single Doxygen Collection) Development
RunParaklete.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 int MyCreateCrsMatrix( char *in_filename, const Epetra_Comm &Comm, 
00051          Epetra_Map *& readMap,
00052          const bool transpose, const bool distribute, 
00053          bool& symmetric, Epetra_CrsMatrix *& Matrix ) {
00054 
00055   Epetra_CrsMatrix * readA = 0; 
00056   Epetra_Vector * readx = 0; 
00057   Epetra_Vector * readb = 0;
00058   Epetra_Vector * readxexact = 0;
00059 
00060   //
00061   //  This hack allows TestOptions to be run from either the test/TestOptions/ directory or from 
00062   //  the test/ directory (as it is in nightly testing and in make "run-tests")
00063   //
00064   FILE *in_file = fopen( in_filename, "r");
00065 
00066   char *filename;
00067   if (in_file == NULL ) 
00068     filename = &in_filename[1] ; //  Strip off ithe "." from
00069          //  "../" and try again 
00070   else {
00071     filename = in_filename ;
00072     fclose( in_file );
00073   }
00074 
00075   symmetric = false ; 
00076   std::string FileName = filename ;
00077 
00078   int FN_Size = FileName.size() ; 
00079   std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
00080   std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
00081 
00082   if ( LastFiveBytes == ".triU" ) { 
00083     // Call routine to read in unsymmetric Triplet matrix
00084     EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx, 
00085                   readb, readxexact) );
00086     symmetric = false; 
00087   } else {
00088     if ( LastFiveBytes == ".triS" ) { 
00089       // Call routine to read in symmetric Triplet matrix
00090       EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, true, Comm, readMap, readA, readx, 
00091               readb, readxexact) );
00092       symmetric = true; 
00093     } else {
00094       if (  LastFourBytes == ".mtx" ) { 
00095   EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( filename, Comm, readMap, 
00096                      readA, readx, readb, readxexact) );   
00097   FILE* in_file = fopen( filename, "r");
00098   assert (in_file != NULL) ;  // Checked in Trilinos_Util_CountMatrixMarket() 
00099   const int BUFSIZE = 800 ; 
00100   char buffer[BUFSIZE] ; 
00101   fgets( buffer, BUFSIZE, in_file ) ;  // Pick symmetry info off of this string 
00102   std::string headerline1 = buffer;
00103 #ifdef TFLOP
00104   if ( headerline1.find("symmetric") < BUFSIZE ) symmetric = true;
00105 #else
00106   if ( headerline1.find("symmetric") != std::string::npos) symmetric = true; 
00107 
00108 #endif
00109   fclose(in_file);
00110 
00111       } else {
00112   // Call routine to read in HB problem
00113   Trilinos_Util_ReadHb2Epetra( filename, Comm, readMap, readA, readx, 
00114                  readb, readxexact) ;
00115   if (  LastFourBytes == ".rsa" ) symmetric = true ; 
00116       }
00117     }
00118   }
00119 
00120 
00121   if ( readb )  delete readb;
00122   if ( readx ) delete readx;
00123   if ( readxexact ) delete readxexact;
00124 
00125   Epetra_CrsMatrix *serialA ; 
00126   Epetra_CrsMatrix *transposeA;
00127 
00128   if ( transpose ) {
00129     transposeA = new Epetra_CrsMatrix( Copy, *readMap, 0 );
00130     assert( CrsMatrixTranspose( readA, transposeA ) == 0 ); 
00131     serialA = transposeA ; 
00132     delete readA;
00133     readA = 0 ; 
00134   } else {
00135     serialA = readA ; 
00136   }
00137 
00138   assert( (void *) &serialA->Graph() ) ;
00139   assert( (void *) &serialA->RowMap() ) ;
00140   assert( serialA->RowMap().SameAs(*readMap) ) ; 
00141 
00142   if ( distribute ) { 
00143     // Create uniform distributed map
00144     Epetra_Map DistMap(readMap->NumGlobalElements(), 0, Comm);
00145 
00146     // Create Exporter to distribute read-in matrix and vectors
00147     Epetra_Export exporter( *readMap, DistMap );
00148     
00149     Epetra_CrsMatrix *Amat = new Epetra_CrsMatrix( Copy, DistMap, 0 );
00150     Amat->Export(*serialA, exporter, Add);
00151     assert(Amat->FillComplete()==0);    
00152     
00153     Matrix = Amat; 
00154     //
00155     //  Make sure that deleting Amat->RowMap() will delete map 
00156     //
00157     //  Bug:  We can't manage to delete map his way anyway,
00158     //        and this fails on tranposes, so for now I just accept
00159     //        the memory loss.
00160     //    assert( &(Amat->RowMap()) == map ) ; 
00161     delete readMap; 
00162     readMap = 0 ; 
00163     delete serialA; 
00164   } else { 
00165 
00166     Matrix = serialA; 
00167   }
00168 
00169 
00170   return 0;
00171 }
00172 
00173 
00174 
00175 
00176 
00177 
00178 // ===================== //
00179 // M A I N   D R I V E R //
00180 // ===================== //
00181 //
00182 // This example compares all the available Amesos solvers
00183 // for the solution of the same linear system. 
00184 //
00185 // The example can be run in serial and in parallel.
00186 //
00187 // Author: Marzio Sala, SNL 9214
00188 // Last modified: Oct-05.
00189 
00190 int main(int argc, char *argv[]) 
00191 {
00192 #ifdef HAVE_MPI
00193   MPI_Init(&argc, &argv);
00194   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00195 #else
00196   Epetra_SerialComm Comm;
00197 #endif
00198 
00199   bool verbose = (Comm.MyPID() == 0);
00200   double TotalResidual = 0.0;
00201 
00202   // Create the Map, defined as a grid, of size nx x ny x nz,
00203   // subdivided into mx x my x mz cubes, each assigned to a 
00204   // different processor.
00205 
00206 #if 0
00207   ParameterList GaleriList;
00208   GaleriList.set("nx", 4);
00209   GaleriList.set("ny", 4);
00210   GaleriList.set("nz", 4 * Comm.NumProc());
00211   GaleriList.set("mx", 1);
00212   GaleriList.set("my", 1);
00213   GaleriList.set("mz", Comm.NumProc());
00214   Epetra_Map* Map = CreateMap("Cartesian3D", Comm, GaleriList);
00215   
00216   // Create a matrix, in this case corresponding to a 3D Laplacian
00217   // discretized using a classical 7-point stencil. Please refer to 
00218   // the Galeri documentation for an overview of available matrices.
00219   // 
00220   // NOTE: matrix must be symmetric if DSCPACK is used.
00221 
00222   Epetra_CrsMatrix* Matrix = CreateCrsMatrix("Laplace3D", Map, GaleriList);
00223 #else
00224   bool transpose = false ; 
00225   bool distribute = false ; 
00226   bool symmetric ; 
00227   Epetra_CrsMatrix *Matrix = 0 ;
00228   Epetra_Map *Map = 0 ;
00229   CreateCrsMatrix( "ibm.triU", Comm, Map, transpose, distribute, &symmetric, Matrix ) ;
00230 
00231 
00232 
00233 
00234 #endif
00235 
00236   // build vectors, in this case with 1 vector
00237   Epetra_MultiVector LHS(*Map, 1); 
00238   Epetra_MultiVector RHS(*Map, 1); 
00239     
00240   // create a linear problem object
00241   Epetra_LinearProblem Problem(Matrix, &LHS, &RHS);
00242 
00243   // use this list to set up parameters, now it is required
00244   // to use all the available processes (if supported by the
00245   // underlying solver). Uncomment the following two lines
00246   // to let Amesos print out some timing and status information.
00247   ParameterList List;
00248   List.set("PrintTiming",true);
00249   List.set("PrintStatus",true);
00250   List.set("MaxProcs",Comm.NumProc());
00251 
00252   std::vector<std::string> SolverType;
00253   SolverType.push_back("Amesos_Paraklete");
00254 
00255   Epetra_Time Time(Comm);
00256   
00257   // this is the Amesos factory object that will create 
00258   // a specific Amesos solver.
00259   Amesos Factory;
00260 
00261   // Cycle over all solvers.
00262   // Only installed solvers will be tested.
00263   for (unsigned int i = 0 ; i < SolverType.size() ; ++i) 
00264   {
00265     // Check whether the solver is available or not
00266     if (Factory.Query(SolverType[i])) 
00267     {
00268       // 1.- set exact solution (constant vector)
00269       LHS.PutScalar(1.0);
00270  
00271       // 2.- create corresponding rhs
00272       Matrix->Multiply(false, LHS, RHS);
00273  
00274       // 3.- randomize solution vector
00275       LHS.Random();
00276  
00277       // 4.- create the amesos solver object
00278       Amesos_BaseSolver* Solver = Factory.Create(SolverType[i], Problem);
00279       assert (Solver != 0);
00280 
00281       Solver->SetParameters(List);
00282 
00283       // 5.- factorize and solve
00284       
00285       Time.ResetStartTime();
00286       AMESOS_CHK_ERR(Solver->SymbolicFactorization());
00287 
00288       // 7.- delete the object
00289       delete Solver;
00290 
00291     }
00292   }
00293 
00294   delete Matrix;
00295   delete Map;
00296 
00297   if (TotalResidual > 1e-9) 
00298     exit(EXIT_FAILURE);
00299 
00300 #ifdef HAVE_MPI
00301   MPI_Finalize();
00302 #endif
00303 
00304   return(EXIT_SUCCESS);
00305 } // end of main()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines