Amesos Package Browser (Single Doxygen Collection) Development
example_AmesosFactory.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 
00031 // This example needs Galeri to generate the linear system.
00032 // You must have configured Trilinos with --enable-galeri
00033 // in order to compile this example
00034 
00035 #ifdef HAVE_MPI
00036 #include "mpi.h"
00037 #include "Epetra_MpiComm.h"
00038 #else
00039 #include "Epetra_SerialComm.h"
00040 #endif
00041 #include "Amesos.h"
00042 #include "Epetra_RowMatrix.h"
00043 #include "Epetra_MultiVector.h"
00044 #include "Epetra_LinearProblem.h"
00045 // following header file and namespace declaration
00046 // are  required by this example to generate the linear system,
00047 // not by Amesos itself.
00048 #include "Galeri_Maps.h"
00049 #include "Galeri_CrsMatrices.h"
00050 using namespace Teuchos;
00051 using namespace Galeri;
00052 
00053 // ==================== //
00054 // M A I N  D R I V E R //
00055 // ==================== //
00056 //
00057 // This example will:
00058 // 1.- create a linear system, stored as an
00059 //     Epetra_LinearProblem. The matrix corresponds
00060 //     to a 5pt Laplacian (2D on Cartesian grid).
00061 //     The user can change the global size of the problem 
00062 //     by modifying variables nx and ny.
00063 // 2.- The linear system matrix, solution and rhs
00064 //     are distributed among the available processors,
00065 //     using a linear distribution. This is for 
00066 //     simplicity only! Amesos can support any Epetra_Map.
00067 // 3.- Once the linear problem is created, we
00068 //     create an Amesos Factory object.
00069 // 4.- Using the Factory, we create the required Amesos_BaseSolver
00070 //     solver. Any supported (and compiled) Amesos
00071 //     solver can be used. If the selected solver
00072 //     is not available (that is, if Amesos has *not*
00073 //     been configured with support for this solver),
00074 //     the factory returns 0. Usually, Amesos_Klu
00075 //     is always available.
00076 // 5.- At this point we can factorize the matrix,
00077 //     and solve the linear system. Only three methods
00078 //     should be used for any Amesos_BaseSolver object:
00079 //     1.- NumericFactorization();
00080 //     2.- SymbolicFactorization();
00081 //     3.- Solve();
00082 // 6.- We note that the header files of Amesos-supported
00083 //     libraries are *not* required in this file. They are
00084 //     actually needed to compile the Amesos library only.
00085 //
00086 // NOTE: this example can be run with any number of processors.
00087 //
00088 // Author: Marzio Sala, SNL 9214
00089 // Last modified: Oct-05.
00090 
00091 int main(int argc, char *argv[]) 
00092 {
00093 
00094 #ifdef HAVE_MPI
00095   MPI_Init(&argc, &argv);
00096   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00097 #else
00098   Epetra_SerialComm Comm;
00099 #endif
00100 
00101   int nx = 100;                  // number of grid points in the x direction
00102   int ny = 100 * Comm.NumProc(); // number of grid points in the y direction
00103   int NumVectors = 1;        // number of rhs's. Amesos
00104                              // supports single or
00105            // multiple RHS.
00106 
00107   // Initializes an Gallery object.
00108   // NOTE: this example uses the Trilinos package Galeri
00109   // to define in an easy way the linear system matrix.
00110   // The user can easily change the matrix type; consult the 
00111   // Galeri documentation for mode details.
00112   //
00113   // Here the problem has size nx x ny, and the 2D Cartesian
00114   // grid is divided into mx x my subdomains.
00115   ParameterList GaleriList;
00116   GaleriList.set("nx", nx);
00117   GaleriList.set("ny", ny);
00118   GaleriList.set("mx", 1);
00119   GaleriList.set("my", Comm.NumProc());
00120 
00121   Epetra_Map* Map = CreateMap("Cartesian2D", Comm, GaleriList);
00122   Epetra_CrsMatrix* Matrix = CreateCrsMatrix("Laplace2D", Map, GaleriList);
00123 
00124   // Creates vectors for right-hand side and solution, and the
00125   // linear problem container.
00126 
00127   Epetra_MultiVector LHS(*Map, NumVectors); LHS.PutScalar(0.0); // zero solution
00128   Epetra_MultiVector RHS(*Map, NumVectors); RHS.Random();       // random rhs
00129   Epetra_LinearProblem Problem(Matrix, &LHS, &RHS);
00130 
00131   // ===================================================== //
00132   // B E G I N N I N G   O F  T H E   AM E S O S   P A R T //
00133   // ===================================================== //
00134 
00135   // Initializes the Amesos solver. This is the base class for
00136   // Amesos. It is a pure virtual class (hence objects of this
00137   // class cannot be allocated, and can exist only as pointers 
00138   // or references).
00139   //
00140   Amesos_BaseSolver* Solver;
00141 
00142   // Initializes the Factory. Factory is a function class (a
00143   // class that contains methods only, no data). Factory
00144   // will be used to create Amesos_BaseSolver derived objects.
00145   //
00146   Amesos Factory;
00147 
00148   // Specifies the solver. String ``SolverType'' can assume one 
00149   // of the following values:
00150   // - Lapack
00151   // - Klu
00152   // - Umfpack
00153   // - Pardiso
00154   // - Taucs
00155   // - Superlu
00156   // - Superludist
00157   // - Mumps
00158   // - Dscpack
00159   // 
00160   std::string SolverType = "Klu";
00161   Solver = Factory.Create(SolverType, Problem);
00162 
00163   // Factory.Create() returns 0 if the requested solver
00164   // is not available
00165   //
00166   if (Solver == 0) {
00167     std::cerr << "Specified solver is not available" << std::endl;
00168     // return ok not to break test harness even if
00169     // the solver is not available
00170 #ifdef HAVE_MPI
00171     MPI_Finalize();
00172 #endif
00173     return(EXIT_SUCCESS);
00174   }
00175 
00176   // Parameters for all Amesos solvers are set through
00177   // a call to SetParameters(List). List is a Teuchos
00178   // parameter list (Amesos requires Teuchos to compile).
00179   // In most cases, users can proceed without calling
00180   // SetParameters(). Please refer to the Amesos guide
00181   // for more details.
00182   // NOTE: you can skip this call; then the solver will
00183   // use default parameters.
00184   //
00185   // Parameters in the list are set using 
00186   // List.set("parameter-name", ParameterValue);
00187   // In this example, we specify that we want more output.
00188   //
00189   Teuchos::ParameterList List;
00190   List.set("PrintTiming", true);
00191   List.set("PrintStatus", true);
00192   
00193   Solver->SetParameters(List);
00194   
00195   // Now we are ready to solve. Generally, users will
00196   // call SymbolicFactorization(), then NumericFactorization(),
00197   // and finally Solve(). Note that:
00198   // - the numerical values of the linear system matrix
00199   //   are *not* required before NumericFactorization();
00200   // - solution and rhs are *not* required before calling
00201   //   Solve().
00202   if (Comm.MyPID() == 0)
00203     std::cout << "Starting symbolic factorization..." << std::endl;
00204   Solver->SymbolicFactorization();
00205   
00206   // you can change the matrix values here
00207   if (Comm.MyPID() == 0)
00208     std::cout << "Starting numeric factorization..." << std::endl;
00209   Solver->NumericFactorization();
00210   
00211   // you can change LHS and RHS here
00212   if (Comm.MyPID() == 0)
00213     std::cout << "Starting solution phase..." << std::endl;
00214   Solver->Solve();
00215   
00216   // you can get the timings here
00217   Teuchos::ParameterList TimingsList;
00218   Solver->GetTiming( TimingsList );
00219   
00220   // you can find out how much time was spent in ...
00221   double sfact_time, nfact_time, solve_time;
00222   double mtx_conv_time, mtx_redist_time, vec_redist_time;
00223 
00224   // 1) The symbolic factorization 
00225   //    (parameter doesn't always exist)
00226   sfact_time = TimingsList.get( "Total symbolic factorization time", 0.0 );
00227 
00228   // 2) The numeric factorization 
00229   //    (always exists if NumericFactorization() is called)
00230   nfact_time = Teuchos::getParameter<double>( TimingsList, "Total numeric factorization time" );
00231 
00232   // 3) Solving the linear system 
00233   //    (always exists if Solve() is called)
00234   solve_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
00235 
00236   // 4) Converting the matrix to the accepted format for the solver
00237   //    (always exists if SymbolicFactorization() is called)
00238   mtx_conv_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
00239 
00240   // 5) Redistributing the matrix for each solve to the accepted format for the solver
00241   mtx_redist_time = TimingsList.get( "Total matrix redistribution time", 0.0 );
00242 
00243   // 6) Redistributing the vector for each solve to the accepted format for the solver
00244   vec_redist_time = TimingsList.get( "Total vector redistribution time", 0.0 );
00245 
00246   // =========================================== //
00247   // E N D   O F   T H E   A M E S O S   P A R T //
00248   // =========================================== //
00249 
00250   // delete Solver. MPI calls can occur.
00251   delete Solver;
00252     
00253   // delete the objects created by Galeri
00254   delete Matrix;
00255   delete Map;
00256 
00257 #ifdef HAVE_MPI
00258   MPI_Finalize();
00259 #endif
00260 
00261   return(EXIT_SUCCESS);
00262 
00263 } // end of main()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines