Amesos Package Browser (Single Doxygen Collection) Development
example_AmesosFactory_Tridiag.cpp
Go to the documentation of this file.
00001 
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //            Trilinos: An Object-Oriented Solver Framework
00006 //                 Copyright (2001) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 #include "Amesos_ConfigDefs.h"
00031 #ifdef HAVE_MPI
00032 #include "mpi.h"
00033 #include "Epetra_MpiComm.h"
00034 #else
00035 #include "Epetra_SerialComm.h"
00036 #endif
00037 #include "Amesos_ConfigDefs.h"
00038 #include "Amesos.h"
00039 #include "Epetra_Map.h"
00040 #include "Epetra_CrsMatrix.h"
00041 #include "Epetra_Vector.h"
00042 #include "Epetra_LinearProblem.h"
00043 #include "Teuchos_ParameterList.hpp"
00044 
00045 // ==================== //
00046 // M A I N  D R I V E R //
00047 // ==================== //
00048 //
00049 // This example will:
00050 // 1.- Create an tridiagonal matrix;
00051 // 2.- Call SymbolicFactorization();
00052 // 3.- Change the numerical values of the matrix;
00053 // 4.- Call NumericFactorization();
00054 // 5.- Set the entries of the RHS;
00055 // 6.- Call Solve().
00056 //
00057 // This example is intended to show the required data
00058 // for each phase. Phase (2) requires the matrix structure only. 
00059 // Phase (4) requires the matrix structure (supposed unchanged 
00060 // from phase (2)) and the matrix data. Phase (6) requires the 
00061 // RHS and solution vector.
00062 //
00063 // This example can be run with any number of processors.
00064 //
00065 // Author: Marzio Sala, SNL 9214
00066 // Last modified: Apr-05.
00067 
00068 int main(int argc, char *argv[]) 
00069 {
00070 
00071 #ifdef HAVE_MPI
00072   MPI_Init(&argc, &argv);
00073   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00074 #else
00075   Epetra_SerialComm Comm;
00076 #endif
00077 
00078   int NumGlobalElements = 100; // global dimension of the problem.
00079 
00080   // ======================================================= //
00081   // B E G I N N I N G   O F   M A T R I X   C R E A T I O N //
00082   // ======================================================= //
00083  
00084   // Construct a Map that puts approximatively the same number of 
00085   // equations on each processor. `0' is the index base (that is,
00086   // numbering starts from 0.
00087   Epetra_Map Map(NumGlobalElements, 0, Comm);
00088 
00089   // Create an empty EpetraCrsMatrix 
00090   Epetra_CrsMatrix A(Copy, Map, 0);
00091 
00092   // Create the structure of the matrix (tridiagonal)
00093   int NumMyElements = Map.NumMyElements();
00094 
00095   // Add  rows one-at-a-time
00096   // Need some vectors to help
00097 
00098   double Values[3];
00099   // Right now, we put zeros only in the matrix.
00100   Values[0] = 0.0;
00101   Values[1] = 0.0;
00102   Values[2] = 0.0;
00103   int Indices[3];
00104   int NumEntries;
00106   int* MyGlobalElements = Map.MyGlobalElements();
00107 
00108   // At this point we simply set the nonzero structure of A.
00109   // Actual values will be inserted later (now all zeros)
00110   for (int i = 0; i < NumMyElements; i++) 
00111   {
00112     if (MyGlobalElements[i] == 0) 
00113     {
00114       Indices[0] = 0; 
00115       Indices[1] = 1;
00116       NumEntries = 2;
00117     }
00118     else if (MyGlobalElements[i] == NumGlobalElements-1) 
00119     {
00120       Indices[0] = NumGlobalElements-1;
00121       Indices[1] = NumGlobalElements-2;
00122       NumEntries = 2;
00123     }
00124     else 
00125     {
00126       Indices[0] = MyGlobalElements[i]-1;
00127       Indices[1] = MyGlobalElements[i];
00128       Indices[2] = MyGlobalElements[i]+1;
00129       NumEntries = 3;
00130     }
00131 
00132     AMESOS_CHK_ERR(A.InsertGlobalValues(MyGlobalElements[i], 
00133           NumEntries, Values, Indices));
00134   }
00135 
00136   // Finish up.
00137   A.FillComplete();
00138 
00139   // =========================================== //
00140   // E N D   O F   M A T R I X   C R E A T I O N //
00141   // =========================================== //
00142 
00143   // Now the matrix STRUCTURE is set. We cannot add
00144   // new nonzero elements, but we can still change the
00145   // numerical values of all inserted elements (as we will
00146   // do later).
00147 
00148   // ===================================================== //
00149   // B E G I N N I N G   O F  T H E   AM E S O S   P A R T //
00150   // ===================================================== //
00151 
00152   // For comments on the commands in this section, please
00153   // see file example_AmesosFactory.cpp.
00154   
00155   Epetra_LinearProblem Problem;
00156 
00157   Problem.SetOperator(&A);
00158 
00159   // Initializes Amesos solver. Here we solve with Amesos_Klu.
00160 
00161   Amesos_BaseSolver * Solver;
00162   Amesos Factory;
00163   Solver = Factory.Create("Amesos_Klu", Problem);
00164 
00165   // Factory.Create() returns 0 if the requested solver
00166   // is not available
00167  
00168   if (Solver == 0) {
00169     std::cerr << "Selected solver is not available" << std::endl;
00170     // return ok not to break the test harness
00171 #ifdef HAVE_MPI
00172     MPI_Finalize();
00173 #endif
00174     exit(EXIT_SUCCESS);
00175   }
00176 
00177   // At this point we can perform the numeric factorization.
00178   // Note that the matrix contains 0's only.
00179 
00180   Solver->SymbolicFactorization();
00181 
00182   // Now, we repopulate the matrix with entries corresponding
00183   // to a 1D Laplacian. LHS and RHS are still untouched.
00184 
00185   for (int i = 0; i < NumMyElements; i++) 
00186   {
00187     if (MyGlobalElements[i] == 0) 
00188     {
00189       Indices[0] = 0;   
00190       Indices[1] = 1;
00191       Values[0]  = 2.0; 
00192       Values[1]  = -1.0;
00193       NumEntries = 2;
00194     }
00195     else if (MyGlobalElements[i] == NumGlobalElements-1) 
00196     {
00197       Indices[0] = NumGlobalElements - 1;
00198       Indices[1] = NumGlobalElements - 2;
00199       Values[0]  = 2.0; 
00200       Values[1]  = -1.0;
00201       NumEntries = 2;
00202     }
00203     else 
00204     {
00205       Indices[0] = MyGlobalElements[i] - 1;
00206       Indices[1] = MyGlobalElements[i];
00207       Indices[2] = MyGlobalElements[i] + 1;
00208       Values[0] = -1.0; 
00209       Values[1] = 2.0; 
00210       Values[2] = -1.0;
00211       NumEntries = 3;
00212     }
00213 
00214     AMESOS_CHK_ERR(A.ReplaceGlobalValues(MyGlobalElements[i], 
00215                                          NumEntries, Values, Indices));
00216   }
00217 
00218   // ... and we can compute the numeric factorization.
00219   Solver->NumericFactorization();
00220 
00221   // Finally, we set up the LHS and the RHS vector (Random().
00222   Epetra_Vector b(Map);
00223   b.Random();
00224   Epetra_Vector x(Map);
00225   x.PutScalar(0.0);
00226   
00227   Problem.SetLHS(&x);
00228   Problem.SetRHS(&b);
00229   
00230   Solver->Solve();
00231 
00232   // Print out the timing information and get it from the solver
00233   Solver->PrintTiming();
00234   
00235   Teuchos::ParameterList TimingsList;
00236   Solver->GetTiming( TimingsList );
00237   
00238   // You can find out how much time was spent in ...
00239   double sfact_time, nfact_time, solve_time;
00240   double mtx_conv_time, mtx_redist_time, vec_redist_time;
00241 
00242   // 1) The symbolic factorization
00243   //    (parameter doesn't always exist)
00244   sfact_time = TimingsList.get( "Total symbolic factorization time", 0.0 );
00245 
00246   // 2) The numeric factorization
00247   //    (always exists if NumericFactorization() is called)
00248   nfact_time = Teuchos::getParameter<double>( TimingsList, "Total numeric factorization time" );
00249 
00250   // 3) Solving the linear system
00251   //    (always exists if Solve() is called)
00252   solve_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
00253 
00254   // 4) Converting the matrix to the accepted format for the solver
00255   //    (always exists if SymbolicFactorization() is called)
00256   mtx_conv_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
00257 
00258   // 5) Redistributing the matrix for each solve to the accepted format for the solver
00259   mtx_redist_time = TimingsList.get( "Total matrix redistribution time", 0.0 );
00260 
00261   // 6) Redistributing the vector for each solve to the accepted format for the solver
00262   vec_redist_time = TimingsList.get( "Total vector redistribution time", 0.0 );
00263 
00264   // =========================================== //
00265   // E N D   O F   T H E   A M E S O S   P A R T //
00266   // =========================================== //
00267 
00268   // ================== //
00269   // compute ||Ax - b|| //
00270   // ================== //
00271 
00272   double residual;
00273 
00274   Epetra_Vector Ax(b.Map());
00275   A.Multiply(false, x, Ax);
00276   Ax.Update(1.0, b, -1.0);
00277   Ax.Norm2(&residual);
00278 
00279   if (!Comm.MyPID())
00280     std::cout << "After AMESOS solution, ||b-Ax||_2 = " << residual << std::endl;
00281 
00282   // delete Solver. Do this before MPI_Finalize()
00283   // as MPI calls can occur in the destructor.
00284   delete Solver;
00285     
00286   if (residual > 1e-5)
00287     return(EXIT_FAILURE);
00288 
00289 #ifdef HAVE_MPI
00290   MPI_Finalize();
00291 #endif
00292   return(EXIT_SUCCESS);
00293 
00294 } // end of main()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines