BLAS Pitfalls with Mac OSX Accelerate Framework
Today I'm going to talk about bugs you can encounter when using the Apple (i.e. Mac OSX) Accelerate framework and BLAS functions like dot product, norm, and sum. I hope that describing these bugs here saves you a headache in the future. This will also be a fun chance to use f2py and the fortran_magic Jupyter extension.
%load_ext fortranmagic
Let's start with a very simple routine written in Fortran to see how
fortran_magic works.
%%fortran
subroutine f(x, y, z)
real(4), intent(in) :: x, y
real(4), intent(out) :: z
z = x + y
end subroutine f
print(f(1, 2))
3.0
We can also work with arrays.
%%fortran
subroutine fa(x, y, z, N)
integer, intent(in) :: N
real(4), dimension(n), intent(in) :: x, y
real(4), dimension(n), intent(out) :: z
z = x + y
end subroutine fa
from numpy import array
a = array([1, 2, 3])
b = array([4, 5, 6])
print(fa(a, b))
[5. 7. 9.]
Now we can try rewriting the above function to call BLAS.
%%fortran
subroutine fb(x, y, z, n)
integer, intent(in) :: n
real(4), dimension(n), intent(in) :: x, y
real(4), dimension(n), intent(out) :: z
real(4), parameter :: fac = 1.0
call scopy(n, x, 1, z, 1)
call saxpy(n, fac, y, 1, z, 1)
end subroutine fb
print(fb(a, b))
[5. 7. 9.]
BLAS has functions, in addition to subroutines. For example, to compute the dot product.
%%fortran
subroutine fc(x, y, res, n)
integer, intent(in) :: n
real(4), dimension(n), intent(in) :: x, y
real(4), intent(out) :: res
res = sdot(n, x, 1, y, 1)
end subroutine fc
print(fc(a, b))
0.0
It seems like something has gone wrong. Did I make a mistake in my fortran routine? Maybe we should debug a little.
%%fortran
subroutine fc(x, y, res, n)
integer, intent(in) :: n
real(4), dimension(n), intent(in) :: x, y
real(4), intent(out) :: res
write(*, *) x
write(*, *) y
write(*, *) dot_product(x, y)
res = sdot(n, x, 1, y, 1)
write(*, *) res
end subroutine fc
print(fc(a, b))
0.0 1.00000000 2.00000000 3.00000000
4.00000000 5.00000000 6.00000000
32.0000000
0.00000000
Nothing seems wrong with the code. What is worse, the intrinsic dot_product works, but BLAS doesn't! There is a trick to fix this though.
%%fortran
subroutine fc(x, y, res, n)
real(8), external :: sdot
integer, intent(in) :: n
real(4), dimension(n), intent(in) :: x, y
real(8), intent(out) :: res
res = sdot(n, x, 1, y, 1)
end subroutine fc
print(fc(a, b))
32.0
It seems I had to coach my Fortran routine to believe that the returned value is double precision. Yet, it's sdot, and if we look at the documentation it should definitely return a single precision floating point number. Let's try some other routines.
%%fortran
subroutine fd(x, res, n)
real(8), external :: snrm2
integer, intent(in) :: n
real(4), dimension(n), intent(in) :: x
real(8), intent(out) :: res
res = snrm2(n, x, 1)
end subroutine fd
print(fd(a))
3.7416574954986572
Indeed, for any of the functions in BLAS LEVEL-1, I think you will suffer through this. As I mentioned before though, this is a problem specific to Apple's Accelerate framework. Let's try linking against openblas instead.
%%fortran --extra '/Users/wddawson/miniconda3/envs/fmagic/lib/libopenblas.0.dylib'
subroutine fe(x, res, n)
integer, intent(in) :: n
real(4), dimension(n), intent(in) :: x
real(4), intent(out) :: res
res = snrm2(n, x, 1)
end subroutine fe
print(fe(a))
3.7416574954986572
If you find yourself calling BLAS directly, you might want to take some care to provide a backup plan if your code is being compiled and executed on a Mac.