Re: linear algebra & LAPACK - questions and thoughts

[Dr . Jürgen Sauermann]

Hi Rowan,

i would like to share some thoughts of mine with you.

First of all some history. As you have correctly noticed, LAPACK was a build
requirement for GNU APL up to including version 1.4. Looking at the
configure.ac in 1.4:

AC_CHECK_LIB([lapack],    [dgelsy_])
AC_CHECK_LIB([blas],      [dcopy_])
#AC_CHECK_LIB([gfortran], [_gfortran_transpose])

So, I only needed in single function dgelsy (in order to implement ⌹).
dgelsy was kindly provided by liblapack, liblapack in turn was dependent
on libblas, and libblas was dependent on libgfortran. At that time I
wasn't sure if these libraries would be available on all platforms. My
suspicion at that time was that libgfortran might limit the use of GNU APL
to GNU-only platforms. Also, installing these libraries was rather tedious
at that time.

In order to simplify the GNU APL build process, I decided to manually translate
the FORTRAN version of dgelsy_ (and all functions that it called, and all
functions that they called...) to C++. Before doing that (which was a major
task) I looked at existing ports of LAPACK to C or C++, but none of them
really convinced me. You can still see the result of the manual translation
in file src/LApack.hh of GNU APL.

So, as a summary, reintroducing a dependency on liblapack is not an idea
that would make me happy. The situations here is different though because
unlike ⌹ which is a built-in primitive of every APL interpreter, while
LAPACK support would be an optional feature that could be dropped on all
platforms that cause problems. Similar to the emacs and SQL functions in
GNU APL today.

A second thought is the following. As you may have noticed, I added a
new primitive in GNU APL which computes the QR factorization of a real
or complex matrix. Before implementing that primitive I browsed through
quite a number of linear algebra publications on the web, most of them
concerning the householder transformation which is the core algorithm of
both ⌹ and of QR factorization. I should at this point admit, that I never
fully understood what my own code in LApack.hh was doing, Rather I
translated every LAPACK function needed from FORTRAN to C++ in a 1:1
fashion, and the result seemed to be OK. For that reason I was never too
happy with the ⌹ implementation in GNU APL even though it apparently does
what it should.

Then, while I was looking into QR factorization, Jay Foad pointed us to
the now available ACM papers. I was formerly looking for some of these
papers because they were cited in the ISO APL standard, but at that time
I could not get them without paying for them (I am not an ACM member).
On the other hand the papers were not that important either. While
browsing through the now freely available papers, I stumbled upon a
paper "THE HOUSEHOLDER ALGORITHM AND APPLICATIONS" written by Garry Helzer
in the 1970s. That paper contains in its Appendix a 12-line APL
function of the algorithm that I was trying to implement. Suddenly,
simply by looking at those 12 lines of APL code, the matter became so
crystal clear to me that I could implement it right away.

The point that I am trying to make here is the following. Suppose we would
create a good translation from (parts of) LApack to APL. The resulting APL
code would then be fairly valuable not only for those who want to solve a
particular problem (like eigenvalues) but also for those who want to
understand the mathematics behind the problem. Such a translation would
not be limited to GNU APL but would be useful for all APL programmers.

Of course a translation of LApack to APL will have a worse performance than
a direct implementation in, say, C/C++. But that is the only drawback that
I can currently see. Given the many advantages that a translation of LApack
to APL could have, that would be my preferred option (= the last option in
your "way forward" list below).


Best Regards,
Jürgen

On 4/12/20 6:08 PM, Rowan Cannaday wrote:

I've been mulling over methods of bringing linear algebra methods into gnu-apl, specifically eigenvalues & eigenvectors (leaving open the possibility for more).

The canonical methods for this are from LAPACK (which was formerly a compilation dependency).
Specifically:
    dgeev
    dsyev
    zheev
    zgeev

There are a few paths forward I can think of:
    - include LAPACK as an optional dependency and wrap the functions (similar to how ⌹ was implemented before it was ported to C++). Possibly with a ⎕FOO wrapper.
    - Port (or find) C/C++ versions of the functions and use the ⎕FX interface
    - Implement the underlying algorithms as APL libraries

Mostly looking for advice on how to proceed, that would make this:
    1. Simple to integrate
    2. 'okay' performance
    3. Possibility to expand to more LAPACK functionality
    4. Not introducing unnecessary compilation challenges

Since this is something I suspect has been considered, I figured I'd reach out to the broader list before attempting something.

Thanks,

- Rowan