I have a matrix A which is of size (n2, n1) and I am multiplying it by a matrix, B, of size (n1, n0). I have identified this single matrix multiplication as the bottleneck in my Fortran code. Out of ~2000 lines of code, this single line takes about 77% of the runtime.

A is a double precision matrix with floating point values. B is, currently, a double precision matrix containing only values 1.0 and 0.0. I can easily make this integer, or even binary, but I was using it as real so that I could preserver precision in matmul(A,B).

What is a better way to perform this matrix multiplication to cut down on runtime?

Before anyone suggests it, I am using DGEMM and compiling with -O3 and -mavx for gfortran, and -O3 with -xhost on ifort.

The largest data I have implemented this program on so far, N = 5000, results in n2 = 1668, n1 = 1701, and n0 = 1631. This algorithm was implemented in Matlab and has shorter runtime. Matlab version is about 2.5 seconds, while this fortran program is about 7 seconds. Since this single matrix multiplication is so large, I'm thinking that Matlab is doing something interesting with the variable types.

I have compiled this with ifort using MKL and am current linking against -lblas and using -fexternal-blas, relying on matmul to perform the underlying BLAS routines. The result of ldd on my binary executable is:

linux-vdso.so.1 =>  (0x00002aaaaaacb000)
liblapack.so.3 => /usr/lib64/atlas/liblapack.so.3 (0x00002aaaaaccd000)
libblas.so.3 => /usr/lib64/libblas.so.3 (0x00002aaaab4f0000)
libgfortran.so.3 => /usr/lib64/libgfortran.so.3 (0x00002aaaab747000)
libm.so.6 => /lib64/libm.so.6 (0x00002aaaaba39000)
libgcc_s.so.1 => /lib64/libgcc_s.so.1 (0x0000003f78c00000)
libc.so.6 => /lib64/libc.so.6 (0x00002aaaabcbe000)
libf77blas.so.3 => /usr/lib64/atlas/libf77blas.so.3 (0x00002aaaac052000)
libcblas.so.3 => /usr/lib64/atlas/libcblas.so.3 (0x00002aaaac272000)
/lib64/ld-linux-x86-64.so.2 (0x00002aaaaaaab000)
libatlas.so.3 => /usr/lib64/atlas/libatlas.so.3 (0x00002aaaac492000)
libpthread.so.0 => /lib64/libpthread.so.0 (0x00002aaaacaee000)

B is structured in the way that it has zeros and ones. The lower left portion (not truly lower triangular) has ones and the upper right portion (not triangular) is zeros.

It appears that the Matlab code is treating the B matrix as logical.

  • $\begingroup$ The right way to link to the MKL these days with ifort is to use -mkl. -lblas may find a system library that is not optimized. $\endgroup$
    – Bill Barth
    Commented Jul 28, 2015 at 2:38
  • $\begingroup$ I only use -lblas when compiling with gfortran. I use specific Intel compilation flags when linking with ifort $\endgroup$
    – drjrm3
    Commented Jul 28, 2015 at 2:40
  • $\begingroup$ a) you can use those same options with gfortran to get the MKL, and b) you might try -mkl with ifort just to be sure. $\endgroup$
    – Bill Barth
    Commented Jul 28, 2015 at 2:43
  • $\begingroup$ I don't see how using real for the integer matrix would give you better precision. You can check for your BLAS using ldd on the binary you compile. $\endgroup$
    – AlexE
    Commented Jul 28, 2015 at 11:47
  • 2
    $\begingroup$ Your latest edit shows you are actually not linking against the MKL, but rather against ATLAS (which I'd personally drop any time in favour of OpenBLAS). You may try to leave out -lblas in linking. $\endgroup$
    – AlexE
    Commented Jul 28, 2015 at 13:26

1 Answer 1


What is the percentage of nonzero entries in $B$? If a high percentage of the entries are 0's, then you might well be better off treating $B$ as a sparse matrix in the multiplication.

You haven't told us the size of the matrices. What are n0 and n1? Or at least their order of magnitude?

You haven't said what implementation of the BLAS you're currently using. Your specific comments about compilation flags suggest that you might be using the reference BLAS implementation, which would be a very poor choice in comparison with optimized cache-aware and multithreaded implementations of the BLAS routines such as ATLAS, OpenBLAS, ACML, and MKL.

  • $\begingroup$ Will update ... $\endgroup$
    – drjrm3
    Commented Jul 28, 2015 at 1:59
  • $\begingroup$ Are you sure that you're running the MKL blas routines multithreaded? What is your setting of OMP_NUM_THREADS? $\endgroup$ Commented Jul 28, 2015 at 2:31
  • $\begingroup$ I am intentionally running this singlethreaded. This program is to be used for high throughput systems and is to be run millions of times, so spending theads on sub linear speed up would not be efficient. I am comparing this to singlethreaded matlab, though, so I believe it can be spedup. $\endgroup$
    – drjrm3
    Commented Jul 28, 2015 at 2:39
  • 1
    $\begingroup$ Are you sure that MATLAB is running single threaded? It's kind of difficult to turn off the multithreading in the MATLAB blas... $\endgroup$ Commented Jul 28, 2015 at 2:41
  • 1
    $\begingroup$ I'm sure. I launch it from command line with the single computation thread flag and watch the cpu usage. $\endgroup$
    – drjrm3
    Commented Jul 28, 2015 at 2:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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