Amesos Package Browser (Single Doxygen Collection) Development
Test_MUMPS/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_Map.h"
00010 #include "Epetra_Vector.h"
00011 #include "Epetra_Time.h"
00012 #include "Epetra_CrsMatrix.h"
00013 #include "Epetra_Util.h"
00014 #include "Amesos_Mumps.h"
00015 #include "Amesos_TestRowMatrix.h"
00016 #include "Teuchos_ParameterList.hpp"
00017 //  using namespace Trilinos_Util;   commented out to resolve bug #1886
00018 
00019 //=============================================================================
00020 bool CheckError(const Epetra_RowMatrix& A,
00021     const Epetra_MultiVector& x,
00022     const Epetra_MultiVector& b,
00023     const Epetra_MultiVector& x_exact)
00024 {
00025   std::vector<double> Norm;
00026   int NumVectors = x.NumVectors();
00027   Norm.resize(NumVectors);
00028   Epetra_MultiVector Ax(x);
00029   A.Multiply(false,x,Ax);
00030   Ax.Update(1.0,b,-1.0);
00031   Ax.Norm2(&Norm[0]);
00032   bool TestPassed = false;
00033   double TotalNorm = 0.0;
00034   for (int i = 0 ; i < NumVectors ; ++i) {
00035     TotalNorm += Norm[i];
00036   }
00037   if (A.Comm().MyPID() == 0)
00038     std::cout << "||Ax - b||  = " << TotalNorm << std::endl;
00039   if (TotalNorm < 1e-5 )
00040     TestPassed = true;
00041   else
00042     TestPassed = false;
00043 
00044   Ax.Update (1.0,x,-1.0,x_exact,0.0);
00045   Ax.Norm2(&Norm[0]);
00046   for (int i = 0 ; i < NumVectors ; ++i) {
00047     TotalNorm += Norm[i];
00048   }
00049   if (A.Comm().MyPID() == 0)
00050     std::cout << "||Ax - b||  = " << TotalNorm << std::endl;
00051 #ifdef HAVE_AMESOS_SMUMPS
00052   if (TotalNorm < 1e-2 )
00053 #else
00054   if (TotalNorm < 1e-5 )
00055 #endif
00056     TestPassed = true;
00057   else
00058     TestPassed = false;
00059 
00060   return(TestPassed);
00061 }
00062 
00063 //=============================================================================
00064 int main(int argc, char *argv[]) {
00065 
00066 #ifdef HAVE_MPI
00067   MPI_Init(&argc, &argv);
00068   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00069 #else
00070   Epetra_SerialComm Comm;
00071 #endif
00072 
00073   int NumGlobalElements = 100;
00074   int NumVectors = 7;
00075   
00076   // =================== //
00077   // create a random map //
00078   // =================== //
00079 
00080   int* part = new int[NumGlobalElements];
00081 
00082   if (Comm.MyPID() == 0) {
00083     Epetra_Util Util;
00084 
00085     for( int i=0 ; i<NumGlobalElements ; ++i ) {
00086       unsigned int r = Util.RandomInt();
00087       part[i] = r%(Comm.NumProc());
00088     }
00089   }
00090 
00091   Comm.Broadcast(part,NumGlobalElements,0);
00092 
00093   // count the elements assigned to this proc
00094   int NumMyElements = 0;
00095   for (int i = 0 ; i < NumGlobalElements ; ++i) {
00096     if (part[i] == Comm.MyPID())
00097       NumMyElements++;
00098   }
00099 
00100   // get the loc2global list
00101   int* MyGlobalElements = new int[NumMyElements];
00102   int count = 0;
00103   for (int i = 0 ; i < NumGlobalElements ; ++i) {
00104     if (part[i] == Comm.MyPID() )
00105       MyGlobalElements[count++] = i;
00106   }
00107 
00108   Epetra_Map Map(NumGlobalElements,NumMyElements,MyGlobalElements,
00109                  0,Comm);
00110 
00111   delete [] part;
00112   
00113   // ===================== //
00114   // Create a dense matrix //
00115   // ===================== //
00116  
00117   Epetra_CrsMatrix Matrix(Copy,Map,NumGlobalElements);
00118 
00119   int* Indices = new int[NumGlobalElements];
00120   double* Values = new double[NumGlobalElements];
00121 
00122   for (int i = 0 ; i < NumGlobalElements ; ++i) 
00123     Indices[i] = i;
00124 
00125   for (int i = 0 ; i < NumMyElements ; ++i) {
00126     int iGlobal = MyGlobalElements[i];
00127     for (int jGlobal = 0 ; jGlobal < NumGlobalElements ; ++jGlobal) {
00128       if (iGlobal >= jGlobal) 
00129   Values[jGlobal] = 1.0 * (jGlobal + 1);
00130       else
00131   Values[jGlobal] = 1.0 * (iGlobal + 1);
00132     }
00133     Matrix.InsertGlobalValues(MyGlobalElements[i],
00134                               NumGlobalElements, Values, Indices);
00135 
00136   }
00137 
00138   Matrix.FillComplete();
00139 
00140   delete [] Indices;
00141   delete [] Values;
00142  
00143   // ======================== //
00144   // other data for this test //
00145   // ======================== //
00146 
00147   Amesos_TestRowMatrix A(&Matrix);
00148   Epetra_MultiVector x(Map,NumVectors);
00149   Epetra_MultiVector x_exact(Map,NumVectors);
00150   Epetra_MultiVector b(Map,NumVectors);
00151   x_exact.Random();
00152   A.Multiply(false,x_exact,b);
00153 
00154   // =========== //
00155   // AMESOS PART //
00156   // =========== //
00157 
00158   Epetra_LinearProblem Problem;
00159   Amesos_Mumps* Solver = new Amesos_Mumps(Problem);
00160 
00161   Problem.SetOperator(&A);
00162   Problem.SetLHS(&x);
00163   Problem.SetRHS(&b);
00164 
00165   Teuchos::ParameterList List;
00166   List.set("MaxProcs",2);
00167   AMESOS_CHK_ERR(Solver->SetParameters(List));
00168 
00169   AMESOS_CHK_ERR(Solver->SymbolicFactorization());
00170   AMESOS_CHK_ERR(Solver->NumericFactorization());
00171   AMESOS_CHK_ERR(Solver->Solve());
00172 
00173   bool TestPassed = true;
00174 
00175   TestPassed = TestPassed &&
00176     CheckError(A,x,b,x_exact);
00177 
00178   delete Solver;
00179 
00180   AMESOS_CHK_ERR( ! TestPassed ) ; 
00181 
00182 #ifdef HAVE_MPI
00183   MPI_Finalize();
00184 #endif
00185 
00186   return(EXIT_SUCCESS);
00187 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines