Amesos Package Browser (Single Doxygen Collection) Development
Test_Singular/cxx_main.cpp
Go to the documentation of this file.
00001 #include "Amesos_ConfigDefs.h"
00002 
00003 #ifdef HAVE_MPI
00004 #include "mpi.h"
00005 #include "Epetra_MpiComm.h"
00006 #else
00007 #include "Epetra_SerialComm.h"
00008 #endif
00009 #include "Epetra_Util.h"
00010 #include "Epetra_Map.h"
00011 #include "Epetra_Vector.h"
00012 #include "Epetra_CrsMatrix.h"
00013 #include "Epetra_LinearProblem.h"
00014 #include "Epetra_Export.h"
00015 #include "Amesos.h"
00016 #include "Teuchos_ParameterList.hpp"
00017 #include "Teuchos_RefCountPtr.hpp"
00018 #include "Galeri_ReadHB.h"
00019 
00020 //============ //
00021 // main driver //
00022 //============ //
00023 
00024 int main(int argc, char *argv[]) 
00025 {
00026 #ifdef HAVE_MPI
00027   MPI_Init(&argc, &argv);
00028   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00029 #else
00030   Epetra_SerialComm Comm;
00031 #endif
00032 
00033 // ================= //
00034   // reading HB matrix //
00035   // ================= //
00036 
00037   // HB files are for serial matrices. Hence, only
00038   // process 0 reads this matrix (and if present
00039   // solution and RHS). Then, this matrix will be redistributed
00040   // using epetra capabilities.
00041   // All variables that begin with "read" refer to the
00042   // HB matrix read by process 0.
00043   Epetra_Map* readMap;
00044   Epetra_CrsMatrix* readA;
00045   Epetra_Vector* readx;
00046   Epetra_Vector* readb;
00047   Epetra_Vector* readxexact;
00048 
00049   // Name of input file with a numerical singularity
00050   std::string matrix_file="bcsstm05_ns.rua";
00051   try
00052   {
00053     Galeri::ReadHB(matrix_file.c_str(), Comm, readMap,
00054                    readA, readx, readb, readxexact);
00055   }
00056   catch(...)
00057   {
00058     std::cout << "Caught exception, maybe file name is incorrect" << std::endl;
00059 #ifdef HAVE_MPI
00060     MPI_Finalize();
00061 #else
00062     // not to break test harness
00063     exit(EXIT_SUCCESS);
00064 #endif
00065   }
00066 
00067   // Create uniform distributed map.
00068   // Note that linear map are used for simplicity only!
00069   // Amesos (through Epetra) can support *any* map.
00070 
00071   Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
00072 
00073   // Create the distributed matrix, based on Map.
00074   Epetra_CrsMatrix A(Copy, map, 0);
00075 
00076   const Epetra_Map &OriginalMap = readA->RowMatrixRowMap() ;
00077   assert (OriginalMap.SameAs(*readMap));
00078   Epetra_Export exporter(OriginalMap, map);
00079 
00080   Epetra_Vector x(map);          // distributed solution
00081   Epetra_Vector b(map);          // distributed rhs
00082   Epetra_Vector xexact(map);     // distributed exact solution
00083 
00084   // Exports from data defined on processor 0 to distributed.
00085   x.Export(*readx, exporter, Add);
00086   b.Export(*readb, exporter, Add);
00087   xexact.Export(*readxexact, exporter, Add);
00088   A.Export(*readA, exporter, Add);
00089   A.FillComplete();
00090 
00091   // deletes memory
00092   delete readMap;
00093   delete readA;
00094   delete readx;
00095   delete readb;
00096   delete readxexact;
00097 
00098   // Creates an epetra linear problem, contaning matrix
00099   // A, solution x and rhs b.
00100   Epetra_LinearProblem problem(&A,&x,&b);
00101 
00102   // =========== //
00103   // AMESOS PART //
00104   // =========== //
00105 
00106   // Create the factory
00107   Amesos factory;
00108 
00109   // Create the solver
00110   std::string solverType = "Klu";
00111   Teuchos::RCP<Amesos_BaseSolver> solver = Teuchos::rcp( factory.Create(solverType, problem) );
00112 
00113   // Perform symbolic factorization
00114   int symRet = solver->SymbolicFactorization();
00115   if (symRet != 0) {
00116     std::cout << "Processor "<< map.Comm().MyPID() << " : Symbolic factorization did not complete!" << std::endl;
00117   }
00118 
00119   // Perform numeric factorization
00120   int numRet = solver->NumericFactorization();
00121   if (numRet != 0) {
00122     std::cout << "Processor "<< map.Comm().MyPID() << " : Numeric factorization did not complete!" << std::endl;
00123   } 
00124 
00125   // Check that all processors returned error -22 (Numerically Singular)!
00126   int minRet = 0;
00127   Comm.MinAll( &numRet, &minRet, 1 );
00128 
00129   if (minRet != NumericallySingularMatrixError) {
00130     if (Comm.MyPID()==0)
00131       std::cout << std::endl << "End Result: TEST FAILED" << std::endl;
00132   }
00133   else {
00134     if (Comm.MyPID()==0)
00135       std::cout << std::endl << "End Result: TEST PASSED" << std::endl;
00136   }
00137  
00138 #ifdef HAVE_MPI
00139   MPI_Finalize();
00140 #endif
00141 
00142   return(0);
00143 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines