Adventures with Fortran

If you do scientific and/or numerical computing, you probably use or at least know about BLAS, and LAPACK. If you've done much programming with them, you may have passed in data that it cannot handle, and then you'll probably know about xerbla_.

xerbla_ is the error-handler in those two libraries, and by default it terminates the whole process. To do otherwise you have to provide a replacement routine, and have the linker for your program insert that replacement, which the library will then use instead. This is increasingly hard or impossible to do, including on macOS, Windows, and static libraries. There is now an open pull-request in the reference versions of those libraries that will permanently solve this, by allowing client programs to just pass in a new handler, which will then be used. No linker stuff involved.

Implementing this involved figuring out how to achieve this in Fortran. As a test-bed, I made a directory with this simple, no-nonsense Makefile:

FC = gfortran
FFLAGS = -g

.f90.o:
        $(FC) $(FFLAGS) -c $< -o $@
.f90:
        $(FC) $(FFLAGS) $< -o $@ $(LDFLAGS)
.SUFFIXES:
.SUFFIXES: .f90 .o

all: callbacks

The first try was using somewhat modern Fortran, with a module:

module xerbla_callbacks
  implicit none
  private
  public :: active_callback, xerbla_interface
  procedure(xerbla_interface), pointer :: active_callback => null()
  interface
    subroutine xerbla_interface(srname, info)
      character*(*), intent(in) :: srname
      integer, intent(in) :: info
    end subroutine
  end interface
end module xerbla_callbacks

subroutine set_xerbla(cb)
  use xerbla_callbacks
  implicit none
  procedure() :: cb
  active_callback => cb
end subroutine set_xerbla

subroutine get_xerbla(cb_ret)
  use xerbla_callbacks
  procedure(xerbla_interface), pointer :: cb_ret
  cb_ret => active_callback
end subroutine get_xerbla

subroutine xerbla(srname, info)
  use xerbla_callbacks
  character*(*), intent(in) :: srname
  integer, intent(in) :: info
  if (associated(active_callback)) then
    call active_callback(srname, info)
  else
    print *, 'I am the main xerbla'
  end if
end subroutine

subroutine xerbla_replacement(srname, info)
  character*(*), intent(in) :: srname
  integer, intent(in) :: info
  print *, 'I am the replacement xerbla'
end subroutine

program hello
  use xerbla_callbacks
  implicit none
  external xerbla_replacement, set_xerbla, xerbla
  procedure(), pointer :: already_cb
  interface
    subroutine get_xerbla(cb_ret)
      procedure(), pointer :: cb_ret
    end subroutine
  end interface
  call xerbla('Hello, World!', 5)
  call get_xerbla(already_cb)
  call set_xerbla(xerbla_replacement)
  call xerbla('Hello, World!', 5)
  call set_xerbla(already_cb)
  call xerbla('Hello, World!', 5)
end program hello

This didn't work in LAPACK because it has an automatic translation facility to turn its source files into versions that suffix routines with _64 that cannot handle modules well. Using an older idiom that allows stored, persistent data was necessary, using the ENTRY keyword.

subroutine xerbla(srname, info)
  character*(*), intent(in) :: srname
  integer, intent(in) :: info
  procedure(xerbla_interface), pointer :: active_callback => null(), cb_ret
  procedure(xerbla_interface) :: cb
  abstract interface
    subroutine xerbla_interface(srname, info)
      character*(*), intent(in) :: srname
      integer, intent(in) :: info
    end subroutine
  end interface
  if (associated(active_callback)) then
    call active_callback(srname, info)
  else
    print *, 'I am the main xerbla'
  end if
  return
  entry set_xerbla(cb)
    active_callback => cb
  return
  entry get_xerbla(cb_ret)
    cb_ret => active_callback
  return
end subroutine

subroutine xerbla_replacement(srname, info)
  character*(*), intent(in) :: srname
  integer, intent(in) :: info
  print *, 'I am the replacement xerbla'
end subroutine

program hello
  implicit none
  external xerbla_replacement, set_xerbla, xerbla
  procedure(), pointer :: already_cb
  interface
    subroutine get_xerbla(cb_ret)
      procedure(), pointer :: cb_ret
    end subroutine
  end interface
  call xerbla('Hello, World!', 5)
  call get_xerbla(already_cb)
  call set_xerbla(xerbla_replacement)
  call xerbla('Hello, World!', 5)
  call set_xerbla(already_cb)
  call xerbla('Hello, World!', 5)
end program hello

Amazingly, this actually all works, and when the C interfaces for these things are figured out, this facility will become available. That will solve a lot of tricky problems, especially for those who make dynamic-language bindings for BLAS etc, but also for those using it on macOS, Windows, etc.

Leave a comment

About Mohawk

user-pic I blog about Perl.