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
.
-mkl
.-lblas
may find a system library that is not optimized. $\endgroup$gfortran
to get the MKL, and b) you might try -mkl withifort
just to be sure. $\endgroup$-lblas
in linking. $\endgroup$