Does anybody know of a good open source implementation of the singular value decomposition written in modern Fortran? I may need to modify an SVD implementation in order to support a custom matrix type that includes forward-mode automatic differentiation via overloaded operators (following the example from Arjen Markus’s book), and it would be very helpful to have an implementation that was close to useable right out of the gate.
Don’t know of any SVD implementations written in Modern Fortran
but you might be better off just finding an old F77 implementation
that has some F90 mods and then wrapping that routine with whatever
Modern Fortran logic you need. Linpack and/or Lapack would be your
best bet (check www.netlib.org). A routine I have used in the past is
a slightly modified version of the SVD routine from the Forsythe, Malcolm,
and Moler book (the modified version is from EISPACK, I think), Computer Methods for Mathematical Computations. The version I used ican be found as part of the FMM
software on Ralph Carmichael’s Public Domain Programs for the Aeronautical Engineer (PDAS)
web site (https://www.pdas.com). Go to Computer Methods for Mathematical Computations. In addition to the SVD routine there are several other useful programs from the FMM book such as zeroin,
quant8 etc. The download file contains the original F77 files plus some moderate refactoring
to F90 by Carmichael. (Note the FMM routine is basically the same one found in Numerical Recipes)
The term “modern Fortran” is not precisely defined, although I use it too. Alan Miller’s dsvdc.f90 uses free source form and has argument INTENTs. From John Burkardt’s site one can compile gfortran -std=legacy linpack_d.f90 lapack_d.f90 svd_test.f90.
@Beliavsky, completely agree about the use (and maybe misuse) of the term “modern Fortran”.
I personally consider anything from F90 on as “modern”. My pet name for all of Fortran (sorry
FORTRAN) before F90 is “Jurrasic Fortran”
When it comes to matrix operations, it is usually more efficient to use analytical matrix derivatives rather than attempt to differentiate the source code. The following reference is very useful in this regard: http://people.maths.ox.ac.uk/~gilesm/files/NA-08-01.pdf.
With that said, the SVD derivatives look quite involved for the orthogonal factors and I have not used them before.
Since you are looking for forward-mode differentiation you may want to consider simply using the complex step method on something like zgesvd. While requiring slightly more operations than a dual-number approach, it may end up being more efficient due to the built-in support for complex numbers in Fortran.
I think it will be tricky to modify a SVD subroutine to make it works with automatic differentiation.
I did something similar to a diagonalisation procedure, it is not trivial. Furthermore, they are a way to use the results (eigenvalues, eigenvectors) of the matrix to get the eigenvalue and eigenvector derivatives.
I’m not sure, it is possible to do it in a similar way for SVD, but probably, it is.