LAPACK interfaces

I’m using a few LAPACK routines and to ensure that I use them properly, I’ve written interface blocks like

    pure subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
      integer, intent(in)                         :: n,nrhs,lda,ldb
      real,    intent(inout), dimension(lda,n)    :: a
      integer, intent(out),   dimension(n)        :: ipiv
      real,    intent(inout), dimension(ldb,nrhs) :: b
      integer, intent(out)                        :: info
    end subroutine dgesv

where the pure attribute is debatable because some routines can write to screen.

Theoretically, the interfaces can be auto-generated from the LAPACK Doxygen documentation, but that would require that it follows a very rigid scheme of notation.

I was wondering if such interfaces are available somewhere, @ivanpribec for example uses

  subroutine factorize_sp(this,info)
    use lapack, only: lapack_factorize => sgetrf
    type(lu_workspace(sp,*)), intent(inout) :: this
    integer, intent(out), optional :: info  
    include "lu_pdt.inc"
  end subroutine

in Purpose of kind parameters in derived types - #2 by ivanpribec.

As far as I see, LAPACK95 -- Fortran95 interface to LAPACK has a similar/the same aim but it seems to be abandoned for 20 years.

2 Likes

We are talking about ~1000 procedures for a complete interface to LAPACK this is a quite daunting task. I’ve my own set of wrappers for the subset <100 procedures I usually use, but not enough to make a good wrapper library. LAPACK95 seems like the best choice at the moment.

Maybe also interesting:

1 Like

I think you can use the Intel Fortran option

gen-interfaces
Tells the compiler to generate an interface block for each routine in a source file.

Typically, I’ve just used my own interfaces extracted manually from the doxygen documentation or the Intel MKL documentation. In most of my applications I would need only between 1 to 6 routines, so this was feasible for me.

I’ve been hesitant to adopt LAPACK95, because of

  1. the extra argument checking that takes place, which I’ve feared would reduce performance (this is probably not justified)
  2. the extra work to unpack and install the library from Netlib

The Netlib and MKL versions of LAPACK95 also slightly differ in the naming of the interfaces, adding an extra complication if you want your code to work with multiple compilers.

The low-level (external) versions on the other hand only require adjustment to the linking sequence and selection of an optimized BLAS library.

There is also the mfi project by @14ngp which uses the fypp preprocessor to automatically generate both the original BLAS/LAPACK interfaces as well as some new higher-level modern interfaces. It’s not yet complete, but it seems like a promising project!

(This repo would benefit from a similar CI workflow as stdlib to perform the preprocessing and make available the full set of f90 files for fpm or direct usage.)

1 Like

Maybe worth mentioning, that there are several libraries which create and expose abstractions for matrices and usually come with a more high-level wrapper for BLAS, LAPACK and ScaLAPACK like operations. An example would be GitHub - lanl/bml: The Basic Matrix Library (bml), which offers multiple sparse and dense storage formats.

1 Like

Many thanks for all the suggestions, it’s good to see that there is some work going on.

The LAPACK docstrings are very nicely formatted and can be parsed to generate interfaces automatically. I’ve written a script that gets the information for most of the routines. Once it is fully operational I will post it here, hopefully it can be used in one of the projects or maybe even included in the official LAPACK release. It works without changing a line of code, I have just unified the style of the comments in a few cases (GitHub - MarDiehl/lapack at standardized-documentation).

Does anyone knows why all array arguments are always assumed size (dimension(*)) even if the shape is explicitly given? LAPACK seem to to put the cart before the horse: The documentation contains the information about dimension and intent and the user needs to make sure that the routine is called correctly. It would be much easier to put that information in the function signature and auto-generate the documentation for the user.

Partially for compatibility with C, it makes it the caller’s responsibility to provide a contiguous array section.

The second and main reason I believe is that Lapack originated in F77, before dynamic memory allocation was widely available upon the arrival of Fortran 90.

Instead you might have just declared a big work-array:

integer, parameter :: n = 3
real :: work(10,11)    ! Sufficient for systems of size 10 × 10
integer :: iwork(10), info
external :: sgetrf, sgetrs

! Fill the top-left n × n block of work
! ...

call sgetrf(n,n,work(1,1),10,iwork,info)

! Right-hand side is stored in the last column
work(1,11) =  1.0
work(2,11) =  2.0
work(3,11) = -3.0

call sgetrs('N',n,1,work(1,1),10,iwork,work(1,11),10,info)
end 

If you add the intent attributes now and try to compile such a legacy code, it will be riddled with compiler warnings and errors…

Edit: here was a semi-related example of a BLAS call to banded matrix vector multiplication. With an explicit interface, it turns out you need to use pointer remapping to match the compiler’s expectations. With an external procedure, you just pass the corresponding first array element.

I got the impression that the development at GitHub - Reference-LAPACK/lapack: LAPACK development repository is still on-going, including the transition to free format source for recently touched routines.

Taking the extra mile and adding intent declarations to all procedures could have good chances to make it into the mainline (maybe ask first in an issue). It might even be worth to contribute the interfaces we craft for our projects again and again as include files there and use them to replace the with external declared procedures.

We have tried to do that here: fortran-utils/lapack.f90 at master · certik/fortran-utils · GitHub, as well as providing “modern Fortran” interface to Lapack built on top: fortran-utils/linalg.f90 at master · certik/fortran-utils · GitHub

Because many of the Lapack routines were originally designed for Fortran 77, there are cases where the notion of “intent” is difficult to interpret and apply.

Consider, as an example, the Lapack95 subroutine call

call gels(A,b,info=info)

The matrix A is m X n, b is m X 1, and info is a scalar. When m > n, this routine returns the least squares solution to the problem A.x = b, but it also returns a compact vector v with the same 2-norm as the residual vector r = A.x - b . The second argument contains the vector b on input. However, on output b is replaced by a m-vector, the first n rows of which contain x, and the remaining m-n rows contain v. Conceptually, b has intent IN, and x and v have intent OUT. Programmatically, however, we have to regard b as having intent INOUT, with content b on input and [x;v] on output.

Now, try to answer the related question of the assumed sizes/shapes of the subroutine arguments, and you will see a number of problems. Similarly, if a wrapper routine in modern Fortran is written, instead of getting by with just two required arguments for GELS, we may need three or four: A, b, x, v. In the current Lapack95 design, the second argument, b, is replaced by the concatenation of x and v.

I don’t think we can have intent(out) in any LAPACK / BLAS routine without changing the behaviour of the API, so this is not an option. In the case of ?gels we would have intent(inout) for both, similar for other cases where it is ambiguous. I don’t really see a fundamental issue in adding intent and explicit interface declaration.

According to the documentation, there are cases where variables have out arguments, for example LAPACK: cgesc2.

I have now a small python script that generates interfaces. It is successful for 1715 of the 2063 files with *.f extension in SRC.

With some fine tuning and using dimension(*) as fallback for the cases where the dimension is not fixed, an even higher ratio should be possible.

Attaches is the interface file, it compiles with gfortran if long lines are allowed, unfortunately it is not standard conforming yet.

lapack.f90 (845.2 KB)

1 Like