I have thoughts below, but my overall questions are:
What is the easiest way to get eigenvectors/eigenvalues of a generic (diagonalizable) matrix in Fortran for the purposes of prototyping? Additionally, for performant code in large codebases that have a large need for this functionality, what is the best way (with good programming practices) to do this?
I recently found, from prototyping an algorithm in Fortran, something to dislike, compared to say Matlab or Python:
The difficulty in finding eigenvalues and eigenvectors of dense matrices. I’ll share my thoughts below, and pose a question.
Context: I was prepping to give a Fortran tutorial (a second part) given to a particular community. Its focus was on scientific computing. An extended example I had in mind included an eigenvalue/eigenvector solve in the algorithm. An online search seemed to tell me LAPACK was the easiest way to go for this. Downloading and compiling LAPACK is fine, linking them a little annoying but fine, but actually using them is quite annoying compared to
Es, Vs = np.linalg.eig(A). I found (in my limited time) there was no good way to introduce this functionality to the people I was giving the tutorial for, so that they could meaningfully follow along themselves during the active tutorial.
The annoyance/difficulty in calling LAPACK routines is a big hamper to my desire to prototype algorithms in Fortran, or teach people anything involving spectra of matrices with Fortran. Fortran’s in-built dot_product and matmul routines are great; it would be great to have an in-built eigen() routine, or an easily importable library with that functionality. It doesn’t have to even be nearly as performant as LAPACK, but it would make prototyping a lot quicker. As it is, eigenvalue calculations are an extremely common part of science, and if I put on my ‘instructor hat’, I would never have students use Fortran to do this when Matlab or Python can do it much more easily. It’s kind of sad when, for eigenvalue/eigenvectors, Matlab/Python really can do “Formula Translation” unlike Fortran, especially when they’re usually calling LAPACK themselves!
To have speed, using LAPACK or equivalent is inevitable. Using Google, it seemed for dense matrices it was best to use LAPACK. For optimal performance, one should perform the workspace query for LAPACK routines, and re-use that information. The workspace also means carrying around an extra two vectors for SSYEV (or more depending on the routine, like SSYEVR), plus several extra scalars. For performance, this is necessary but unwieldy.
For maintenance and readability of large codebases
Many of Fortran’s unique features as a language are useful for increasing readability of large codebases. Declarations of intent(in/out), separation of functions and subroutines, existence of the ‘pure’ and ‘elemental’ keywords, ability to name constructs and exit nested ones at any level of nesting, etc. It almost seems like Fortran is built for HPC with large codebases.
This is amazing, as large codebases need a lot of help in code readability for future modification and maintenance. LAPACK routines, however, are an obscure language of their own. For readability, I’d need to wrap the Fortran routines I need in a nicer wrapper that users/future programmers don’t need to worry about. For performance, however, I need to keep all the extra workspace-related baggage that LAPACK has. I’d wrap the variables into a derived type, probably, and instantiate different objects if I’m using LAPACK routines on different sized matrices, but re-using the workspace variables when I can. The names would be chosen for readability. However, writing my own wrappers and types every time seems silly. At least high-performant libraries like FFTW have you create a ‘plan’ containing the equivalent information for this. Creating a ‘compute plan’ is also now a large part of the CUDA framework, including recent features like CUDA Graphs. Perhaps this could be explored.
I’m not sure what the Fortran community can do about LAPACK in particular, but if people are thinking of adding eigenvalue/eigenvector solvers to the Fortran Package Manager (or standard library), or at least making modern wrappers to LAPACK easily available in fpm, the above may be something to think about.
Again, my questions are about: the easiest way to get eigenvalues/eigenvectors for prototyping in Fortran (is it really LAPACK?), and the best way to do it for large codebases that need good and robust performance.