Implement a fast SVD routine

[romeric]

At the moment Fastor provides the fastest LU decomposition, inversion and solve routines for small tensors but lacks a fast SVD routine

[dynamic-queries]

I am interested in contributing to this issue.
Did you have anything specific in mind?

[romeric]

Absolutely, go ahead!
The only requirement we have is that you need to add unittests for it and it has to be sufficiently fast. Look at the signature of lu and qr function for reference here.

A very naive, inefficient and perhaps buggy solution would be (just for your reference)

template<typename T, size_t M, size_t N>
inline void svd(const Tensor<T,M,N> &a, Tensor<T,M,M> &u, Tensor<T,N,M> &s, Tensor<T,N,N> &v, const T tol=1e-10) {

    // reserve space in advance
    size_t maxiter = 100*std::max(M,N);
    size_t counter = 0;
    u.eye2();
    s = transpose(a);
    v.eye2();
    Tensor<T,M,M> q;
    T err = std::numeric_limits<T>::max();
    while (err > tol && counter < maxiter) {
        qr(transpose(s),q,s); u = u % q;
        qr(transpose(s),q,s); v = v % q;
        T F = sum(diag(s));
        if (std::abs(F) < tol) F = 1;
        err = E/F;
        counter++;
    }
    return;
}