LAPACK interfaces

1

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
2

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
3

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.

4

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.

5

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
6

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
7

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.

8

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.

9

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.

10

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

11

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.

12

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.

13

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

14

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
15

The Numerical Algorithms Group appears to have a large collection of interface_blocks, as well as examples and a test suite.
https://github.com/numericalalgorithmsgroup/LAPACK_examples

1 Like
16

Looks great! But it seems that they are provided for double precision only?

17

One possible reason is that LAPACK was developed in the f77 era. In f77, I think it was not allowed to have a zero dimension for an explicit shape dummy argument, and some compilers had runtime tests for such violations. Linear algebra codes, including LAPACK, sometimes need to allow an N=0 argument as a way to signify that the array should be ignored or that certain calculations within the routine should be skipped. So the f77 workaround was to allow the N=0 argument, but to use an assumed size dummy argument declaration. Most of the LAPACK matrices have a separate leading dimension argument, so the dummy argument is declared as A(LDA,*) for an array that would otherwise have been declared as A(LDA,N). In those cases, LDA is required to be positive even when the N=0 argument results in ignoring the array. That LDA requirement is also due to the f77 dummy argument restriction. Another common convention in LAPACK is that a work array is passed to computational routines that is declared as WORK(*) and effectively treated as WORK(LWORK) when LWORK is positive. This is because LWORK=-1 is treated as a special case of a workspace query; in this case, only WORK(1) is referenced within the routine.

That f77 dummy argument restriction was removed in f90, where zero-length arrays (and character strings) became legitimate. Other than the shortlived LAPACK95 interface, the LAPACK routines themselves have never been updated beyond the f77 fortran interface. There are many practical reasons why this would have been a good idea, including the use of optional arguments for work arrays, the use of derived types to simplify the argument lists, explicit interfaces, and generic interfaces to eliminate floating point KIND dependencies in the user code.

18

LAPACK is often called by languages other than Fortran, such as C, C++, Python/SciPy, Matlab, and R, in which case the Fortran 77 interface is ok.

I have not seen much use of LAPACK95, as you say, even by Fortran programmers. Are there technical reasons for that? Fortran programmers who want a better interface for LAPACK typically write their own wrappers.

19

This is one of the first lines at the netlib site:

# NOTE:  This directory was updated in November, 2000!

Nowadays, someone would see that and assume it is outdated code. Even the name LAPACK95 suggests that it is outdated. Originally, 2000 was about the time that people were first developing f90/f95 code, so many people were just getting started using modern features of the language when LAPACK95 development ended. For some context, g95 development started in 2000, and gfortran was forked in 2003.

20

I’ve used it a few times via the MKL interface. But I’ve always avoided the Netlib distribution because of the extra work to unpack the tarball, adjust the Makefile, and link it in my own project build system. There is also the burden of extra logic to choose between MKL and the Netlib version in your own build system.