4 ms·

I don't get the whole premise of this article. If they wanted to minimize the use of "antiquated" Fortran code, that's fine I guess. But why didn't they just u

I don't get the whole premise of this article. If they wanted to minimize the use of "antiquated" Fortran code, that's fine I guess. But why didn't they just use ZGEEV for both the real and the complex case? Instead, they went on great lengths to transform complex matrices into real matrices, and recover the eigenvalues of the original complex matrix somehow...

ogogmad 4y ago
The eigenvalues and eigenvectors are going to be complex some of the time anyway. So I agree, why not just copy-and-paste the real algorithm (or use operator overloading) while remembering to insert complex-conjugates? It's educational, sure. But it also seems unnecessary.
- reikonomusha 4y ago
  I think it's a very non-trivial task to "just" copy and paste a three-decade-old program that has tens of thousands of lines of dependencies, making sure to remember to add conjugates as necessary. :) Eigenvalues and eigenvectors are notoriously difficult programs to get right for general matrices with no known properties a priori.
reikonomusha 4y ago
Getting DGEEV to work requires 50,000 lines of mechanically translated code to work from HOMPACK, BLAS, and LAPACK. At no point is the "COMPLEX*16" data type used in the public interface to DGEEV because the routine demands real and imaginary parts are split into separate columns of the eigenvalue array and eigenvector matrix, i.e., actual complex numbers are never constructed. So things are "easy". As such, ZGEEV would require (1) getting COMPLEX*16 values and arrays to translate appropriately into Common Lisp (where the representation is as a boxed, heap-allocated object, typically), and (2) introducing approximately double the amount of code (totaling well over 100,000 lines of mechanically translated code... just to get one function!). To (1), I think F2CL supports this, but it's evident from the longest-maintained distribution of the software that it hasn't been as thoroughly exercised as the real-type code. So, if LOC doesn't cost anything, this might be a fine starting point. Lastly, "antiquated" specifically refers to a 15+-year-old LAPACK distribution written in FORTRAN 77, for which F2CL works—sort of. F2CL is actually more strict than F77, which makes translation of real-world F77 code difficult. For instance, in F77, you can pass an entire string to a function that declares itself as accepting only a single character. However, after translation, Lisp will raise a type error when that happens. So, our choice is to fix 30-year-old F77-to-Lisp translation code, or manually fix those Fortran type errors in the output. The latter is what was done to make the real code work. Reference LAPACK is now written in Fortran 90+, with many numerical bug fixes in each release. As such, this cannot be translated with this tool, unless somebody writes a full-blown F90 parser and refactors the translator. If one wants to use modern LAPACK directly in MAGICL, there's already a MAGICL/EXT-LAPACK system one can load, but it requires your system to have the binary library installed.
- amelius 4y ago
  But the lapack functions have been thoroughly tested, over decades. Don't roll your own lapack functions is generally good advice.
  - reikonomusha 4y ago
    LAPACK has been written in Fortran 90+ since 2008 in version 3.2. That's nearly 15 years ago. The last release of the LAPACK reference was 3.10.1 released April 2022. If your constrained to using Common Lisp, you can - load and distribute LAPACK binaries (already well supported in principle, though Windows users typically have difficulties) - fix the F77 translator and produce millions of lines of Common Lisp code for all LAPACK types (or, at worst, >100k for ZGEEV and dependencies) - write a new F90+ translator and produce millions of lines of Common Lisp code for all LAPACK types - write a numerically stable and accurate implementation from scratch that works for all manner of conditioned matrices (really hard!) - translate a small subset of reliable F77 routines that clock in at a few 10k LOC and prove mathematically they can do more than what's written on the tin. The first and last options are what have been chosen, the first allowing the most modern implementation of numerical functions for all data types, and the last allowing expedient independence from foreign binaries (i.e., it's pure Common Lisp). Which option are you suggesting is optimal if you're seeking to relieve dependence—as a fallback for when it's unavailable—on a few shared libraries (/ foreign code)?
    - amelius 4y ago
      Use a Fortran compiler to compile the specific parts that you need, and link statically if that is what you want?
      - reikonomusha 4y ago
        This is essentially the first option I listed, which is already supported by MAGICL by loading MAGICL/EXT-LAPACK [1]. Since MAGICL is a library and not a compiled application, using LAPACK in the way you suggest is still a foreign dependency (now a C compiler, Fortran compiler, Lisp compiler, and the transitive dependencies of the *GEEV Fortran functions), even if at the end of the day it's basically packaged up in an end-user's application binary. (It's probably out of the scope of this thread to get into the weeds of how Common Lisp is compiled as a standard language, and how specific compiler implementations like SBCL work. It suffices to say though that it's almost always not a simple matter of "get .o files and link them to your compiled Lisp", in no small part due to the popularity of image-based development.) It's perfectly reasonable to consider different deployment strategies, but "statically link Fortran code (essentially into the Lisp runtime)", to a user of MAGICL, is not an ordinary thing to do. [1] https://github.com/quil-lang/magicl#extensions https://github.com/quil-lang/magicl#extensions

ogogmad 4y ago

The eigenvalues and eigenvectors are going to be complex some of the time anyway. So I agree, why not just copy-and-paste the real algorithm (or use operator overloading) while remembering to insert complex-conjugates? It's educational, sure. But it also seems unnecessary.

reikonomusha 4y ago
I think it's a very non-trivial task to "just" copy and paste a three-decade-old program that has tens of thousands of lines of dependencies, making sure to remember to add conjugates as necessary. :) Eigenvalues and eigenvectors are notoriously difficult programs to get right for general matrices with no known properties a priori.

reikonomusha 4y ago

Getting DGEEV to work requires 50,000 lines of mechanically translated code to work from HOMPACK, BLAS, and LAPACK. At no point is the "COMPLEX*16" data type used in the public interface to DGEEV because the routine demands real and imaginary parts are split into separate columns of the eigenvalue array and eigenvector matrix, i.e., actual complex numbers are never constructed. So things are "easy". As such, ZGEEV would require (1) getting COMPLEX*16 values and arrays to translate appropriately into Common Lisp (where the representation is as a boxed, heap-allocated object, typically), and (2) introducing approximately double the amount of code (totaling well over 100,000 lines of mechanically translated code... just to get one function!). To (1), I think F2CL supports this, but it's evident from the longest-maintained distribution of the software that it hasn't been as thoroughly exercised as the real-type code. So, if LOC doesn't cost anything, this might be a fine starting point. Lastly, "antiquated" specifically refers to a 15+-year-old LAPACK distribution written in FORTRAN 77, for which F2CL works—sort of. F2CL is actually more strict than F77, which makes translation of real-world F77 code difficult. For instance, in F77, you can pass an entire string to a function that declares itself as accepting only a single character. However, after translation, Lisp will raise a type error when that happens. So, our choice is to fix 30-year-old F77-to-Lisp translation code, or manually fix those Fortran type errors in the output. The latter is what was done to make the real code work. Reference LAPACK is now written in Fortran 90+, with many numerical bug fixes in each release. As such, this cannot be translated with this tool, unless somebody writes a full-blown F90 parser and refactors the translator. If one wants to use modern LAPACK directly in MAGICL, there's already a MAGICL/EXT-LAPACK system one can load, but it requires your system to have the binary library installed.

amelius 4y ago
But the lapack functions have been thoroughly tested, over decades. Don't roll your own lapack functions is generally good advice.
- reikonomusha 4y ago
  LAPACK has been written in Fortran 90+ since 2008 in version 3.2. That's nearly 15 years ago. The last release of the LAPACK reference was 3.10.1 released April 2022. If your constrained to using Common Lisp, you can - load and distribute LAPACK binaries (already well supported in principle, though Windows users typically have difficulties) - fix the F77 translator and produce millions of lines of Common Lisp code for all LAPACK types (or, at worst, >100k for ZGEEV and dependencies) - write a new F90+ translator and produce millions of lines of Common Lisp code for all LAPACK types - write a numerically stable and accurate implementation from scratch that works for all manner of conditioned matrices (really hard!) - translate a small subset of reliable F77 routines that clock in at a few 10k LOC and prove mathematically they can do more than what's written on the tin. The first and last options are what have been chosen, the first allowing the most modern implementation of numerical functions for all data types, and the last allowing expedient independence from foreign binaries (i.e., it's pure Common Lisp). Which option are you suggesting is optimal if you're seeking to relieve dependence—as a fallback for when it's unavailable—on a few shared libraries (/ foreign code)?
  - amelius 4y ago
    Use a Fortran compiler to compile the specific parts that you need, and link statically if that is what you want?
    - reikonomusha 4y ago
      This is essentially the first option I listed, which is already supported by MAGICL by loading MAGICL/EXT-LAPACK [1]. Since MAGICL is a library and not a compiled application, using LAPACK in the way you suggest is still a foreign dependency (now a C compiler, Fortran compiler, Lisp compiler, and the transitive dependencies of the *GEEV Fortran functions), even if at the end of the day it's basically packaged up in an end-user's application binary. (It's probably out of the scope of this thread to get into the weeds of how Common Lisp is compiled as a standard language, and how specific compiler implementations like SBCL work. It suffices to say though that it's almost always not a simple matter of "get .o files and link them to your compiled Lisp", in no small part due to the popularity of image-based development.) It's perfectly reasonable to consider different deployment strategies, but "statically link Fortran code (essentially into the Lisp runtime)", to a user of MAGICL, is not an ordinary thing to do. [1] https://github.com/quil-lang/magicl#extensions https://github.com/quil-lang/magicl#extensions