Performance optimization or tuning possible for Scalapack Gemm?

Question

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:

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;
}

Answer 1

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.

Answer 2

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).

Answer 3

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.