linalg: Singular Value Decomposition

[perazz]

Singular Value Decomposition

of $A = U \cdot S \cdot V^T$, for real and complex matrices, based on LAPACK *GESDD functions.

• [x] base implementation
• [x] tests
• [x] exclude unsupported xdp
• [x] documentation
• [x] submodule
• [x] examples
• [x] subroutine interface

Prior art

• NumPy: svd(a, full_matrices=True, compute_uv=True, hermitian=False)
• Scipy: svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, check_finite=True, lapack_driver='gesdd')
• Scipy: svdvals(a, overwrite_a=False, check_finite=True)

Proposed implementation

• Subroutine interface call svd(a,s [,u,vt,overwrite_a,full_matrices,err])

• Function interface s = svdvals(a [, err])

GESDD allows to overwrite a with either u or v vectors, saving storage and matrix copies. The internal GESDD task is chosen based on user inputs in order to minimize the number of allocations in particular:

• Whether u/v are provided, and if a can be overwritten
• Whether storage is for full matrices or partial matrices (min(m,n))

cc: @jvdp1 @jalvesz @zoziha: I believe this is ready for consideration. I'm investigating what is causing ifx to segfault, although I don't have a local installation (on Apple silicon) EDIT: segfault fixed (ifx compiler bug)

[perazz]

ifx has an internal compiler error on the test: will need to investigate

[perazz]

Fixed:

• ifx has a compiler bug setting complex(real128) values from the result of an arithmetic operation, as shown in Compiler Explorer
• discussion at Fortran Discourse
• Issue solved in 9c2f133; the outstanding CI failure is fixed by #809

[jalvesz]

@perazz I just started to read the PR, brilliant! I would like to bring an idea for discussion regarding the naming for optional variables and their internal counterparts. I've seen elsewhere and started to stick to the following naming style:

subroutine my_sub( a, some_optional )
 <type>, intent(<>) :: a
 <type>, intent(in), optional :: some_optional
 <type>                 :: some_optional_ !> same name as the user-exposed name with the "_" suffix
end subroutine

I found this useful for two reasons (1) less head-scratching to find a name for the corresponding internal counterpart (2) easy to find the variable or change their names if needed.

Following the use of optionals, I also found that is easier to read the implementations using the following assignment pattern (say for a logical):

some_optional_ = .false.
if(present(some_optional)) some_optional_ = some_optional

Which I feel, makes it clear what is the default at a glance when reading the code.

I've being thinking about the overwrite_a boolean flag. This design currently implies:

• If .false., the matrix is not touched, but an internal allocation is incurred.
• If .true., no internal allocation is incurred, but the matrix is lost as it is recycled as a temp buffer.

What if instead of managing this choice in this manner, an optional temp array is allowed? such that:

• if(present(temp)) the matrix a is untouched, temp is used (and recycled) for internal processing.
• else internal allocation is done, a remains untouched.

[perazz]

Appreciate your reviewing @jalvesz. Though I'had stuck to the similar logic overwrite_a -> copy_a of the other linear algebra routines, I agree it's not super clear. So I am fixing it here and will update the other instances to have the same new approach.

Regarding optional present arguments, things were a bit tricky in the SVD because the LAPACK backend allows to return this data on top of A. So if A can be destroyed (i.e. it is an internal copy, or, the user said it can be destroyed), and U or V is not requested on output (not present) but must be computed anyways (the other matrix was requested), we don't want to allocate another temporary for it. We can just spare this additional allocation by pointing the routine to replace it on top of A.

So yes, it is a bit hard to read in the code, because the backend logic is also pretty cumbersome. But ultimately, the goal here was to strive for performance i.e. to avoid allocations wherever possible

[perazz]

Thank you @jvdp1 @jalvesz for your reviews, it looks much better now!

[perazz]

Thank you @jvdp1. I'll wait for a few more days, and if there are no comments, will merge this one. Thank you all for the reviews.