[Trilinos-Users] performance of trilinos

[chenliang wang]

Hi,
I use the code listed at the end of this mail to test the speedup using Trilinos (10.10.2 compiled with mkl) in our DELL cluster.To my surprise, there is no significant speedup no matter how many processors I use.(almost 1.5x~1.2x)
Is that possible? Is there any more rough benchmarks to compare the speedup ?

Chen-Liang Wang

***********************************************************************************************************************************************************************************************************************************************************************
//the example code:
// Solve a 2D Laplacian problem
// This example builds the matrix and solves it with AztecOO.
#include "Epetra_ConfigDefs.h"
#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "AztecOO.h"

// external function
void get_neighbours( const int i, const int nx, const int ny,
              int & left, int & right,
              int & lower, int & upper);

// =========== //
// main driver //
// =========== //

int main(int argc, char *argv[])
{

#ifdef HAVE_MPI
MPI_Init(&argc, &argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif

// number of nodes along the x- and y-axis
int nx = 4806;
int ny = 4046;
int NumGlobalElements = nx * ny;

// create a linear map
Epetra_Map Map(NumGlobalElements,0,Comm);

// local number of rows
int NumMyElements = Map.NumMyElements();
// get update list
int * MyGlobalElements = new int [NumMyElements];
Map.MyGlobalElements( MyGlobalElements );

// Create a Epetra_Matrix with 5 nonzero per rows

Epetra_CrsMatrix A(Copy,Map,5);

// Add rows one-at-a-time
// Need some vectors to help

double Values[4];
int Indices[4];
int NumEntries;
int left, right, lower, upper;
double diag = 4.0;

for( int i=0 ; i<NumMyElements; ++i ) {
    int NumEntries=0;
    get_neighbours( MyGlobalElements, nx, ny,
             left, right, lower, upper);
    if( left != -1 ) {
    Indices[NumEntries] = left;
    Values[NumEntries] = -1.0;
    ++NumEntries;
    }
    if( right != -1 ) {
      Indices[NumEntries] = right;
      Values[NumEntries] = -1.0;
      ++NumEntries;
    }
    if( lower != -1 ) {
      Indices[NumEntries] = lower;
      Values[NumEntries] = -1.0;
      ++NumEntries;
    }
    if( upper != -1 ) {
      Indices[NumEntries] = upper;
      Values[NumEntries] = -1.0;
      ++NumEntries;
    }
    // put the off-diagonal entries
    A.InsertGlobalValues(MyGlobalElements, NumEntries, Values, Indices);
    // Put in the diagonal entry
    A.InsertGlobalValues(MyGlobalElements, 1, &diag, MyGlobalElements+i);
}

// Finish up
A.FillComplete();

// create x and b vectors
Epetra_Vector x(Map);
Epetra_Vector b(Map);
//b.PutScalar(100.0);
b.Random();
x.Random();
// ==================== AZTECOO INTERFACE ======================

// create linear problem  Ax=b
Epetra_LinearProblem Problem(&A,&x,&b);
// create AztecOO instance 
AztecOO Solver(Problem);

Solver.SetAztecOption( AZ_precond, AZ_none );//
//Solver.SetAztecOption( AZ_precond, AZ_Jacobi );//Jacobi
Solver.SetAztecOption( AZ_solver, AZ_cg );//CG
//Solver.SetAztecOption( AZ_solver, AZ_gmres );//
Solver.SetAztecOption(AZ_output,AZ_all);//
Solver.Iterate(100,1E-4);//

// ==================== END OF AZTECOO INTERFACE ================

if( Comm.MyPID() == 0 ) {
    cout << "Solver performed " << Solver.NumIters()
    << "iterations.\n";
    cout << "Norm of the true residual = " << Solver.TrueResidual() << endl;
}

#ifdef HAVE_MPI
MPI_Finalize();
#endif

return(EXIT_SUCCESS);

}

/****************************************************************************/
/****************************************************************************/
/****************************************************************************/

void get_neighbours( const int i, const int nx, const int ny,
              int & left, int & right,
              int & lower, int & upper)
{

int ix, iy;
ix = i%nx;
iy = (i - ix)/nx;

if( ix == 0 )
    left = -1;
else
    left = i-1;
if( ix == nx-1 )
    right = -1;
else
    right = i+1;
if( iy == 0 )
    lower = -1;
else
    lower = i-nx;
if( iy == ny-1 )
    upper = -1;
else
    upper = i+nx;

return;

}