RFC: Naming of functions in an Lapack module

[dmbates]

@ViralBShah has suggested that all the code that ends up calling

ccall(dlsym(_jl_liblpack, ...) ...)

which includes much of the current base/linalg_lapack.jl and base/linalg_factorizations.jl be moved to an Lapack module. See the discussion of pull request #1243.

After doing so the _jl_lapack prefixes on the names of functions could be removed. I can start work on this but I wanted to request comments on some function names first:

1. Remove all _jl_lapack prefixes
2. Perhaps leave all symbol names unexported. This way these functions will only be called by a user sufficiently knowledgeable to reference an unexported symbol name.
3. Add ! to the names of those functions that modify an input array.

I would also suggest continuing the convention of passing StridedMatrix or VectorOrMatrix templated types and not passing arguments like lda, m, n, etc. that can be derived from the other arguments. (The underlying Lapack routines just get a pointer to the contents of the vector or matrix and must have dimensions, leading dimensions, etc. passed separately.)

I don't know how much cost there is in adding many small functions to the namespace of a module. If this cost is negligible then I would suggest creating a more-or-less complete set of interfaces to the Lapack subroutines. For example, in addition to the QR decompositions allow for the LQ, RQ and QL decompositions.

Finally, @timholy has created a Tridiagonal type in pull request #1243. We may want to consider other specialized types of matrix structures reflecting the specialized Lapack types, such as Triangular, Symmetric, GeneralBanded, Hessenberg, Bidiagonal, Schur, etc.

[timholy]

I'm fine with the idea of a module in principle. To me, matrix algebra is part of the core of Julia, and I would not want to have to type

julia> load("lapack.jl")

julia> import LapackMod.*

at the beginning of every Julia session (unless there is some significant cost to Julia's potential for mainstream programming by having Lapack available at startup). @ViralBShah , care to comment further on the reasons for a separate module?

Re the naming, those plans seem good. I think we could do most of that independently of the idea for a separate module, since Base is already module and those are not exported.

While working on the last set of changes, I came to the conclusion that the low-level Lapack wrappers should be moved from factorizations.jl to linalg_lapack.jl. There are already wrappers for some Lapack factorization routines in linalg_lapack.jl, this would just make that file the go-to place for all Lapack wrappers.

[ViralBShah]

That is right. I am suggesting to put the LAPACK wrappers in a module to start with.

Linear Algebra could be in a module too - but I agree that it would be really annoying to load a module to do lu(), and it should always be loaded by default. If such a module is going to be loaded by default, then it might as well just be in Base.

-viral

On 03-Sep-2012, at 8:55 PM, dmbates wrote:

@ViralBShah has suggested that all the code that ends up calling

ccall(dlsym(_jl_liblpack, ...) ...)
which includes much of the current base/linalg_lapack.jl and base/linalg_factorizations.jl be moved to an Lapack module. See the discussion of pull req[https://github.com//pull/1243]

—
Reply to this email directly or view it on GitHub.

[dmbates]

@timholy The idea is that the user would not typically access functions in the Lapack or BLAS modules. The linear algebra functions would stay essentially as they are - it is just that behind the scenes they would call, e.g.

Lapack.getrf!(...)

instead of the current _jl_lapack_getrf

I agree about the low-level Lapack wrappers in factorizations.jl being moved to linalg_lapack.jl. That was part of the plan although I might now have expressed that well.

[StefanKarpinski]

I wholeheartedly support this plan. The non-surface BLAS and LAPACK functions should all be squirreled away in modules for sure.

[dmbates]

During the implementation I encountered order-dependency problems and ended up modifying the functions to interface with Lapack tridiagonal subroutines. I wanted to assure @timholy that it was not just on a whim that I changed his interfaces.

I have added a file Lapack.jl defining a module of the Lapack interface functions. Tim had Lapack interface functions defined in linalg_specialized.jl that had methods for Tridiagonal objects. If these are to be moved to the Lapack.jl file then in sysimg.jl it must follow the file containing the declaration of the Tridiagonal type or contain the declaration itself. That becomes awkward. I imposed the rule that direct calls to the Lapack interface functions use only scalars, vectors and matrices as arguments. When a specialized matrix type is defined, methods for eig, solve, \, lud, qrd, etc. access fields and pass them to the Lapack interface routines.

I also changed the declaration of info in Tim's code. IIRC, to allow the info value to be written by the Lapack subroutine it must be an Array(Int32, 1), not an Int32. If it is an Int32 and passed as &info the Lapack subroutine can change the value but the change will not be visible in the calling function.

[StefanKarpinski]

These sound like good things to me. Can we see these changes somewhere?

[timholy]

If it is an Int32 and passed as &info the Lapack subroutine can change the value but the change will not be visible in the calling function.

Didn't know that, thanks for explaining.

As for the other changes you describe, I agree that it's probably best to separate the Tridiagonal type from the most basic Lapack wrapper. Assuming you're fixing both the core functions and the higher-level wrappers, I'm all in favor. Thanks for the cleanup!

[dmbates]

@StefanKarpinski Changes will be committed to git@github.com:dmbates/julia in a LapackModule branch once I get them building again. As you know it gets tricky creating modules in files that are included in sysimg.jl

[StefanKarpinski]

Oh, yes I do. Good reason to move more things out of sysimg.jl in my opinion. Maybe we need to have a third bootstrapping stage. Viral's been struggling with the same sort of things for his refactoring of RNG stuff.

[dmbates]

@StefanKarpinski The LapackModule branch is now on my fork of the julia repository. Still much cleaning up to do.

[dmbates]

A further clean up, though perhaps not today, will be to replace all the two-state flag arguments to Lapack.* (like diag, uplo, ...) that currently take LapackChar arguments with Bool arguments so we don't need to check for legitimate values.

Longer term: Should wrappers for BLAS be in a separate module or should these modules be combined?

[StefanKarpinski]

Is there a clear logical distinction between what would go into the BLAS module versus the LAPACK module? E.g. are the calling conventions significantly different? I know there's a clear historical distinction, but that doesn't necessarily mean there's also a clear logical distinction. OpenBLAS, for example, includes LAPACK functionality from what I can tell.

[dmbates]

There is no clear distinction between BLAS and Lapack subroutines. The naming conventions are the same but the names themselves don't overlap.

I don't think there would be a problem in putting both sets of interfaces in one module.

[StefanKarpinski]

Sounds good. Should probably be called LAPACK then and include BLAS routines as a lower level.

[dmbates]

I have found that Jeff's allowing a pointer argument in ccall to be constructed as &foo is very helpful in calling Fortran subroutines. I have gone even further and condensed sequences like

   ldb = stride(B,2)
   ccall(dlsym(...), Void,
          (Ptr{Int32}), &ldb)

to

   ccall(dlsym(...), Void,
         (Ptr{Int32}), &stride(B,2))

and that does seem to work. (Clever fellow that @JeffBezanson !)

I just wanted to check that I should be able to expect it to continue to work as desired.