5
$\begingroup$

I'm comparing the performance of distributed gemm, using Scalapack over OpenBLAS, with threaded gemm, using OpenBLAS. It seems quite hard for me to get scalapack to give better results than multithreaded BLAS. Is there some tweaking/optimization/configuration I need to do, or something I need to do differently somehow?

My hardware:

  • 4 blades in a chassis, which I will call 'nodes' henceforth
  • connected by gigabit ethernet, not the fastest I know, maybe the cause of the poor results here?
  • each node has 12 physical cores
  • software: mpich2 2.3, Ubuntu 12.04 64-bit, gfortran

I'm comparing three software configurations, for a GEMM of two square matrices, of dimension n:

  • multi-threaded BLAS, using OpenBLAS, with 12 threads. 1 single process, on 1 single node
  • scalapack over OpenBLAS, 48 mpi processes, 1 thread per mpi process, on 4 nodes
  • scalapack over OpenBLAS, 4 mpi processes, 12 threads per mpi process, on 4 nodes

Here is a graph of some results:

Scalapack versus multithreaded BLAS

This graph seems to show:

  • using pure undistributed BLAS is actually super fast
  • using distributed scalapack, with 1 mpi process per node gives the best performance, but is not dramatically faster than pure undistributed BLAS, even for really large matrices (the largest n in the graph is 30,000)
  • using distributed scalapack with 1 thread per mpi process gives the worst results, always worse than simple undistributed BLAS

Some questions:

  • how typical are these results?
  • is there something I'm doing terribly wrong, or could do better, that might improve the results of the distributed scalapack results?
  • is it normal that the 1 thread per mpi process configuration gives worse results than the other configurations tested?

Here is the test code used:

For GEMM:

#include <cstdio>
#include <iostream>
#include <cmath>
using namespace std;

#include "mycblas.h"
#include "utils/NanoTimer.h"
#include "utils/stringhelper.h"
#include "args.h"

extern "C" {
    void openblas_set_num_threads(int num_threads);
}

int main( int argc, char *argv[] ) {
    int N, its, threads;
    Args( argc, argv ).arg("N", &N).arg("its", &its).arg("threads",&threads).go();

    openblas_set_num_threads( threads );    

    NanoTimer timer;
    double *A = (double*)malloc(sizeof(double)*N*N);
    double *B = (double*)malloc(sizeof(double)*N*N);
    int linsize = N * N;
    for( int i = 0; i < linsize; i++ ) {
        A[i] = i + 3;
        B[i] = i * 2;
    }
    int m = N;
    int n = N;
    int k = N;
    double alpha = 1;
    double beta = 0;
    double *C = (double*)malloc(sizeof(double)*N*N);
    timer.toc("setup input matrices");
    for( int it = 0; it < its; it++ ) {
        dgemm(false,false,N,N,N, 1, A, N, B, N, 0, C, N );
        timer.toc("it " + toString(it) );
    }
    int sum = 0;
    for( int mult = 0; mult < log(N)/log(10); mult++ ) {
        int offset = pow(10,mult);
        sum += C[offset];
    }
    cout << "sum, to prevent short-cut optimization " << sum << endl;
    return 0;
}

For scalapack:

#include <iostream>
#include <stdexcept>
#include <cstring>
#include <cmath>
using namespace std;

#include "mpi.h"

#include "utils/NanoTimer.h"
#include "utils/stringhelper.h"
#include "args.h"
#include "scalapack.h"

extern "C" {
    void openblas_set_num_threads(int num_threads);
}

int getRootFactor( int n ) {
    for( int t = sqrt(n); t > 0; t-- ) {
        if( n % t == 0 ) {
            return t;
        }
    }
    return 1;
}

// conventions:
// M_ by N_ matrix block-partitioned into MB_ by NB_ blocks, then
// distributed according to 2d block-cyclic scheme

// based on http://acts.nersc.gov/scalapack/hands-on/exercise3/pspblasdriver.f.html

int main( int argc, char *argv[] ) {
    int p, P;
    blacs_pinfo( &p, &P );
//    mpi_print( toString(p) + " / " + toString(P) );

    int n;
    int numthreads;
    int its;
    int blocksize;
    Args( argc, argv ).arg("N", &n ).arg("num iterations", &its ).arg("numthreads", &numthreads ).arg("blocksize", &blocksize).go();
    openblas_set_num_threads( numthreads );

    int nprows = getRootFactor(P);
    int npcols = P / nprows;
    if( p == 0 ) cout << "grid: " << nprows << " x " << npcols << endl;

    int system = blacs_get( -1, 0 );
    int grid = blacs_gridinit( system, true, nprows, npcols );
    if( p == 0 ) cout << "system context " << system << " grid context: " << grid << endl;

    int myrow, mycol;
    blacs_gridinfo( grid, nprows, npcols, &myrow, &mycol );
//    mpi_print("grid, me: " + toString(myrow) + ", " + toString(mycol) );

    if( myrow >= nprows || mycol >= npcols ) {
//        mpi_print("not needed, exiting");
        blacs_gridexit( grid );
        blacs_exit(0);
        exit(0);
    }

    // A     B       C
    // m x k k x n = m x n
    // nb: blocksize

    // nprows: process grid, number rows
    // npcols: process grid, number cols
    // myrow: process grid, our row
    // mycol: process grid, our col
    int m = n;
    int k = n;
//    int nb = min(n,128); // nb is column block size for A, and row blocks size for B
    int nb=min(n/P,128);

    int mp = numroc( m, nb, myrow, 0, nprows ); // mp number rows A owned by this process
    int kp = numroc( k, nb, myrow, 0, nprows ); // kp number rows B owned by this process
    int kq = numroc( k, nb, mycol, 0, npcols ); // kq number cols A owned by this process
    int nq = numroc( n, nb, mycol, 0, npcols ); // nq number cols B owned by this process
//    mpi_print( "mp " + toString(mp) + " kp " + toString(kp) + " kq " + toString(kq) + " nq " + toString(nq) );

    struct DESC desca, descb, descc;
    descinit( (&desca), m, k, nb, nb, 0, 0, grid, max(1, mp) );
    descinit( (&descb), k, n, nb, nb, 0, 0, grid, max(1, kp) );
    descinit( (&descc), m, n, nb, nb, 0, 0, grid, max(1, mp) );
//    mpi_print( "desca.LLD_ " + toString(desca.LLD_) + " kq " + toString(kq) );
    double *ipa = new double[desca.LLD_ * kq];
    double *ipb = new double[descb.LLD_ * nq];
    double *ipc = new double[descc.LLD_ * nq];

    for( int i = 0; i < desca.LLD_ * kq; i++ ) {
        ipa[i] = p;
    }
    for( int i = 0; i < descb.LLD_ * nq; i++ ) {
        ipb[i] = p;
    }

    if( p == 0 ) cout << "created matrices" << endl;
    double *work = new double[nb];
    if( n <=5 ) {
        pdlaprnt( n, n, ipa, 1, 1, &desca, 0, 0, "A", 6, work );
        pdlaprnt( n, n, ipb, 1, 1, &descb, 0, 0, "B", 6, work );
    }

    NanoTimer timer;
    for( int it = 0; it < its; it++ ) {
        pdgemm( false, false, m, n, k, 1,
                      ipa, 1, 1, &desca, ipb, 1, 1, &descb,
                      1, ipc, 1, 1, &descc );
        MPI_Barrier( MPI_COMM_WORLD );
        if( p == 0 ) timer.toc("it " + toString(it) + " pdgemm");
    }

    blacs_gridexit( grid );
    blacs_exit(0);

    return 0;
}

Other files:

mycblas.h:

extern "C" {
#define ADD_
    // blas:
    #include <cblas_f77.h>

    // lapack:
    void dpotrf_( char *uplo, int *n, double *A, int *lda, int *info );
    void dtrtrs_( char *uplo, char *trans, char *diag, int *n, int *nrhs, double *A, int *lda,
         double *B, int *ldb, int *info );
    void dtrsm_( char *side, char *uplo, char *transA, char *diag, const int *m, const int *n, 
                 const double *alpha, const double *A, const int *lda, double *B, const int *ldb );
}

char boolToChar( bool value ) {
    return value ? 't' : 'n';
}

// double, general matrix multiply
void dgemm( bool transa, bool transb, int m, int n, int k, double alpha, double *A,
    int lda, double *B, int ldb, double beta, double *C, int ldc ) {
    char transachar = boolToChar( transa );
    char transbchar = boolToChar( transb );
    dgemm_(&transachar,&transbchar,&m,&n,&k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc );    
}

// double, triangular, solve matrix
void dtrsm( bool XA, bool isUpper, bool transA, bool isUnitTriangular, int m, int n,
       double alpha, double *A, int lda, double *B, int ldb ) {
    char sideChar = XA ? 'R' : 'L';
    char isUpperChar = isUpper ? 'U' : 'L';
    char transAChar = transA ? 'T' : 'N';
    char isUnitTriangularChar = isUnitTriangular ? 'U' : 'N';
    dtrsm_( &sideChar, &isUpperChar, &transAChar, &isUnitTriangularChar, &m, &n,
        &alpha, A, &lda, B, &ldb );
}

// double, symmetric positive definite, triangular factorization (=cholesky)
int dpotrf( bool isUpper, int N, double *A, int lda ) {
    int info;
    char uplo = isUpper ? 'U' : 'L';
    dpotrf_( &uplo, &N, A, &lda, &info );
    return info;
}

// double, triangular, triangular solve
int dtrtrs( bool isUpper, bool transA, bool isUnitTriangular, int n, int nrhs, double *A, int lda,
             double *B, int ldb ) {
    int info;
    char isUpperChar = isUpper ? 'U' : 'L';
    char transChar = transA ? 'T' : 'N';
    char isUnitTriangularChar = isUnitTriangular ? 'U' : 'N';
    dtrtrs_( &isUpperChar, &transChar, &isUnitTriangularChar, &n, &nrhs, A, &lda, B, &ldb, &info );
    return info;
}

scalapack.h:

#pragma once

extern "C" {
    struct DESC{
        int DTYPE_;
        int CTXT_;
        int M_;
        int N_;
        int MB_;
        int NB_;
        int RSRC_;
        int CSRC_;
        int LLD_;
    } ;

    void blacs_pinfo_( int *iam, int *nprocs );
    void blacs_get_( int *icontxt, int *what, int *val );
    void blacs_gridinit_( int *icontxt, char *order, int *nprow, int *npcol );
    void blacs_gridinfo_( int *context, int *nprow, int *npcol, int *myrow, int *mycol );
    void blacs_gridexit_( int *context );
    void blacs_exit_( int *code );

    int numroc_( int *n, int *nb, int *iproc, int *isrcproc, int *nprocs );
    void descinit_( struct DESC *desc, int *m, int *n, int *mb, int *nb, int *irsrc, int *icsrc, int *ictxt, int *lld, int *info );
    void pdlaprnt_( int *m, int *n, double *a, int *ia, int *ja, struct DESC *desca, int *irprnt,
        int *icprnt, const char *cmatnm, int *nout, double *work, int cmtnmlen );
    void pdgemm_( char *transa, char *transb, int *m, int *n, int *k, double *alpha,
         double *a, int *ia, int *ja, struct DESC *desca, double *b, int *ib, int *jb,
        struct DESC *descb, double *beta, double *c, int *ic, int *jc, struct DESC *descc );
}

void blacs_pinfo( int *p, int *P ) {
    blacs_pinfo_( p, P );
}

int blacs_get( int icontxt, int what ) {
    int val;
    blacs_get_( &icontxt, &what, &val );
    return val;
}

int blacs_gridinit( int icontxt, bool isColumnMajor, int nprow, int npcol ) {
    int newcontext = icontxt;
    char order = isColumnMajor ? 'C' : 'R';
    blacs_gridinit_( &newcontext, &order, &nprow, &npcol );
    return newcontext;
}

void blacs_gridinfo( int context, int nprow, int npcol, int *myrow, int *mycol ) {
    blacs_gridinfo_( &context, &nprow, &npcol, myrow, mycol );
}

void blacs_gridexit( int context ) {
    blacs_gridexit_( &context );
}

void blacs_exit( int code ) {
    blacs_exit_( &code );
}

int numroc( int n, int nb, int iproc, int isrcproc, int nprocs ) {
    return numroc_( &n, &nb, &iproc, &isrcproc, &nprocs );
}

void descinit( struct DESC *desc, int m, int n, int mb, int nb, int irsrc, int icsrc, int ictxt, int lld ) {
    int info;
    descinit_( desc, &m, &n, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info );
    if( info != 0 ) {
        throw runtime_error( "non zero info: " + toString( info ) );
    }
//    return info;
}

void pdlaprnt( int m, int n, double *A, int ia, int ja, struct DESC *desc, int irprnt,
    int icprnt, const char *cmatnm, int nout, double *work ) {
    int cmatnmlen = strlen(cmatnm);
    pdlaprnt_( &m, &n, A, &ia, &ja, desc, &irprnt, &icprnt, cmatnm, &nout, work, cmatnmlen );
}

void pdgemm( bool isTransA, bool isTransB, int m, int n, int k, double alpha,
     double *a, int ia, int ja, struct DESC *desca, double *b, int ib, int jb,
    struct DESC *descb, double beta, double *c, int ic, int jc, struct DESC *descc ) {
    char transa = isTransA ? 'T' : 'N';
    char transb = isTransB ? 'T' : 'N';
    pdgemm_( &transa, &transb, &m, &n, &k, &alpha, a, &ia, &ja, desca, b, &ib, &jb,
        descb, &beta, c, &ic, &jc, descc );
}

utils/NanoTimer.h:

#pragma once
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <sys/time.h>

class NanoTimer {
public:
   struct timespec start;

   NanoTimer() {
      clock_gettime(CLOCK_MONOTONIC,  &start);

   }
   double elapsedSeconds() {
      struct timespec now;
      clock_gettime(CLOCK_MONOTONIC,  &now);
      double time = (now.tv_sec - start.tv_sec) + (double) (now.tv_nsec - start.tv_nsec) * 1e-9;
      start = now;
      return time;
   }
    void toc(string label) {
        double elapsed = elapsedSeconds();
        cout << label << ": " << elapsed << "s" << endl;        
    }
};

utils/stringhelper.h:

#pragma once

#include <vector>
#include <string>
#include <sstream>
#include <iostream>
#include <cstdlib>

template<typename T>
std::string toString(T val ) { // not terribly efficient, but works...
   std::ostringstream myostringstream;
   myostringstream << val;
   return myostringstream.str();
}

args.h:

#pragma once

// usage:
// int N, its;
// arg( "N", &N );
// arg( "its", &its );
// args( argc, argv );

class Arg {
public:
    virtual void assign( const char *argvalue ) = 0;
    virtual void print( ostream &os ) const = 0;
};
ostream &operator<<( ostream &os, const Arg &arg ) {
    arg.print( os );
    return os;
}

class IntArg : public Arg {
public:
    int *argptr;
    IntArg( int *_argptr ) : argptr(_argptr ) {
    }
    void assign( const char *argvalue ) {
        *argptr = atoi( argvalue );
    }
    void print( ostream &os ) const {
        os << (*argptr );
    }
};

vector<string> argnames;
vector<Arg *> argptrs;

void arg_usage(string cmd) {
    cout << "Usage: " << cmd;
    for( int i = 0; i < argnames.size(); i++ ) {
        cout << " [" << argnames[i] << "]";
    }
    cout << endl;
    exit(1);
}

void arg( string name, int *p_value ) {
    argnames.push_back(name);
    argptrs.push_back( new IntArg( p_value ) );
}

void args( int argc, char *argv[] ) {
    if( argc - 1 != argnames.size() ) {
        arg_usage(argv[0]);
    }
    for( int i = 0; i < argnames.size(); i++ ) {
        argptrs[i]->assign( argv[i+1] );
        cout << argnames[i] << ": " << (*argptrs[i]) << endl;
    }
}

class Args {
public:
    int argc;
    char **argv;
    Args( int _argc, char *_argv[] ) : argc(_argc), argv(_argv) {
    }
    void go() {
        args( argc, argv );
    }
    Args &_( string name, int *pvalue ) {
        ::arg( name, pvalue );
        return *this;
    }
    Args &arg( string name, int *pvalue ) {
        ::arg( name, pvalue );
        return *this;
    }
};

Edit: in answer to Jeff's question about my Elemental code, here is my Elemental code for GEMM:

#include "mpi.h"

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <sys/time.h>
using namespace std;

#include "elemental.hpp"
using namespace elem;

extern "C" {
    void openblas_set_num_threads(int num_threads);
}

class NanoTimer {
public:
   struct timespec start;

   NanoTimer() {
      clock_gettime(CLOCK_MONOTONIC,  &start);

   }
   double elapsedSeconds() {
      struct timespec now;
      clock_gettime(CLOCK_MONOTONIC,  &now);
      double time = (now.tv_sec - start.tv_sec) + (double) (now.tv_nsec - start.tv_nsec) * 1e-9;
      start = now;
      return time;
   }
    void toc(string label) {
        double elapsed = elapsedSeconds();
        cout << label << ": " << elapsed << endl;        
    }
};

int sum = 0;
void readMatrix( Matrix<double> &A ) {
    for( int i = 1; i < A.Width(); i *= 10 ) {
        sum += A.Get(i,i);
    }
}
void readMatrix( DistMatrix<double,CIRC,CIRC> &A ) {
    for( int i = 1; i < A.Width(); i *= 10 ) {
        sum += A.Get(i,i);
    }
}

int main( int argc, char *argv[] ) {
    elem::Initialize( argc, argv );
    int p = mpi::CommRank(mpi::COMM_WORLD);
    int P = mpi::CommSize(mpi::COMM_WORLD);

    if( argc < 3 ) {
        if( p == 0 ) cout << "Usage: " << argv[0] << " [N] [multithreaded: 1|0]" << endl;
        return -1;
    }
    int n = atoi( argv[1] );
    int multithreaded = atoi( argv[2] );

    Matrix<double> A(n,n);
    Matrix<double> B(n,n);
    for( int i = 0; i < n; i++ ) {
        for( int j = 0; j < n; j++ ) {
            A.Set(i,j, i*j + 2 );
            B.Set(i,j, i*j + 4 );
        }
    }

    Matrix<double> C(n,n);
    MPI_Barrier(MPI_COMM_WORLD);
    NanoTimer timer;
    if( p == 0 ) {
        Gemm<double>(NORMAL, NORMAL, 1, A, B, 0, C );  
        // read some values, to prevent being optimized out :-P
        timer.toc("blas multithreaded");
        readMatrix(C);

        if( !multithreaded ) {
            openblas_set_num_threads(1);
            Gemm<double>(NORMAL, NORMAL, 1, A, B, 0, C );    
            timer.toc("blas singlethreaded");
            readMatrix(C);
        }
    }
    if( !multithreaded ) {
        openblas_set_num_threads(1);
    }

    Grid g;
    DistMatrix<double,CIRC,CIRC> Aroot(n,n,g);
    DistMatrix<double,CIRC,CIRC> Broot(n,n,g);
    DistMatrix<double,CIRC,CIRC> Croot(n,n,g);
    Aroot.SetRoot(0);
    Broot.SetRoot(0);
    Croot.SetRoot(0);

    if( p == 0 ) {
        for( int i = 0; i < n; i++ ) {
            for( int j = 0; j < n; j++ ) {
                Aroot.Set(i,j,i*j+2);
                Broot.Set(i,j,i*j+2);
            }
        }
    }
    if( p == 0 ) timer.toc("populate root node");

//    DistMatrix<double,MC,STAR> Adist( Aroot );
//    DistMatrix<double,STAR,MR> Bdist( Broot );
//    DistMatrix<double,MC,MR> Cdist(n,n,g);
    DistMatrix<double> Adist( Aroot );
    DistMatrix<double> Bdist( Broot );
    DistMatrix<double> Cdist(n,n,g);
    MPI_Barrier(MPI_COMM_WORLD);
    if( p == 0 ) timer.toc("distributed to slaves");

    Gemm<double>(NORMAL, NORMAL, 1, Adist, Bdist, 0, Cdist );    
    MPI_Barrier(MPI_COMM_WORLD);
    if( p == 0 ) timer.toc("distmatrix gemm");

    Croot = Cdist;
    MPI_Barrier(MPI_COMM_WORLD);
    if( p == 0 ) timer.toc("gathered to master");
    readMatrix(Croot);

    if( p == 0 ) cout << "sum, to prevent optimization out: " << sum << endl;

    elem::Finalize();
    return 0;
}
$\endgroup$
5
$\begingroup$

There's nothing surprising about these results. Matrix multiplication is well known to be communication intensive, and you've got a relatively slow communications network between your four nodes.

Using MPI between two processes on the same node is certainly faster than using MPI between processes on different nodes because you don't have the bandwidth limitations of the ethernet. However, you still have to pay a price in calling the MPI library. In comparison, the multithreaded BLAS gets very fast communication between threads because the memory is fast and because you don't pay the overhead of calling the MPI library routines and using the OS to pass messages between threads as you would with MPI.

Ideally, for large $n$, you should see the 4 nodes, 1 process per node, 12 threads per process method coming in at four times as fast as 12 threads on a single node. The two curves are diverging towards that ratio of 4, but not very fast (you're 3 times as fast at n=30000)- this is an indication that the bandwidth of the gigabit ethernet is really holding you back.

$\endgroup$
  • $\begingroup$ Matrix multiplication is not communication intensive for any reasonable definition of communication intensive. It's completely compute-bound in a reasonable implementation. Yes, it involves communication but the surface-to-volume ratio, i.e. bytes communicated vs. flops, is very favorable. Maybe I misunderstand your definition of communication intensive. $\endgroup$ – Jeff Jul 4 '13 at 1:56
  • $\begingroup$ "Compute bound" vs. "Communication bound" depends very much on how man flops you've got and how many GB/s of bandwidth you have. It's clear that communication (broadly including the block size choices made by ScaLapack and the efficiency of the MPI implementation as well as the gigabit ethernet) is keeping the four node speed at well under four times the one node speed. All of the implementations should be be performing the same number of arithmetic operations so the difference comes down to communication costs. $\endgroup$ – Brian Borchers Jul 4 '13 at 3:36
  • $\begingroup$ @cjordan1 BLAS1 operations are memory-bound. Suggesting in any way that they are compute bound is indication of either an extreme misunderstanding of computing or misuse of relevant terminology. $\endgroup$ – Jeff Jul 6 '13 at 19:54
  • $\begingroup$ @BrianBorchers Yes, load-store is faster than ethernet, but that's not what I was objecting to. $\endgroup$ – Jeff Jul 6 '13 at 19:55
2
$\begingroup$

These results seem normal. Your network is terrible for distributed BLAS. I can't say about the mpich2 shared memory implementation, but it might also not be up to snuff and limiting your 1 thread per process/12 processes per node results. You might try OpenMPI or MVAPICH2 (though I've never run the latter with just Ethernet).

$\endgroup$
  • 1
    $\begingroup$ MVAPICH is based upon MPICH and will not have a superior shared memory implementation, nor is the implementation for TCP/IP any better. In all performance comparisons for shared memory and TCP/IP that I have seen, MPICH is faster than OpenMPI. $\endgroup$ – Jeff Jul 4 '13 at 1:58
  • $\begingroup$ On the shared memory side, I was thinking of LiMIC kernel module that the MVPAICH folks developed. I guess that's technically not shared memory, but a separate mechanism for doing intra-node communications that saves a copy. I'm 99.9% certain that MPICH hasn't adopted this approach, and it would be helpful to Hugh for the intra-node part of his problem. Also, when I said Ethernet, I was thinking of the RoCE and iWARP work that the MVAPICH folks have done and not TCP/IP. I can't tell from the question whether or not that will help Hugh since I don't know if his hardware supports it. $\endgroup$ – Bill Barth Jul 4 '13 at 16:48
  • $\begingroup$ runtime.bordeaux.inria.fr/knem fits nicely into 0.1% :-) $\endgroup$ – Jeff Jul 6 '13 at 20:00
  • $\begingroup$ Interesting. Both MVAPICH and MPICH use KNEM, but only MVAPICH uses LiMIC. Now I'll have to find out what the performance differences are. $\endgroup$ – Bill Barth Jul 7 '13 at 0:19
  • $\begingroup$ I would also be curious how OpenMPI performs with KNem and XPMEM if you're going to go down this path. XPMEM is built into SGI and Cray Linux but I don't know how hard it is to enable in generic systems. $\endgroup$ – Jeff Jul 10 '13 at 20:00
2
$\begingroup$

The algorithm in ScaLAPACK is not very good. Look at SUMMA (http://www.cs.utexas.edu/ftp/techreports/tr95-13.pdf) and its implementation in Elemental (code.google.com/p/elemental/) or CTF (http://ctf.eecs.berkeley.edu/) instead. Papers by Solomonik and coworkers show substantial performance improvements over ScaLAPACK.

A colleague of mine - a physicist, not a computer scientist - implemented SUMMA himself and found it was faster than ScaLAPACK, so it's not that difficult to beat the latter.

$\endgroup$
  • $\begingroup$ I tried Elemental, but all my results were better with Scalapack than Elemental, on both Gigabit ethernet, and Infiniband QR. $\endgroup$ – Hugh Perkins Jul 4 '13 at 7:57
  • $\begingroup$ Without seeing the benchmark code, I can't comment on that, but that's the opposite of what we found for QDR IB and Blue Gene/P. $\endgroup$ – Jeff Jul 6 '13 at 19:56
  • $\begingroup$ added my Elemental code to the question description, at the bottom. $\endgroup$ – Hugh Perkins Jul 7 '13 at 2:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.