Function JacobiSVD

src/backend/cpu/homography.cpp:66–136

template<typename T, int M, int N>
void JacobiSVD(T* S, T* V) {
    const int iterations = 30;
    array<T, N> d{};

    for (int i = 0; i < N; i++) {
        T sd = 0;
        for (int j = 0; j < M; j++) {
            T t = S[i * M + j];
            sd += t * t;
        }
        d[i] = sd;

        V[i * N + i] = 1;
    }

    for (int it = 0; it < iterations; it++) {
        bool converged = false;

        for (int i = 0; i < N - 1; i++) {
            for (int j = i + 1; j < N; j++) {
                T* Si = S + i * M;
                T* Sj = S + j * M;
                T* Vi = V + i * N;
                T* Vj = V + j * N;

                T p = static_cast<T>(0);
                for (int k = 0; k < M; k++) { p += Si[k] * Sj[k]; }

                if (abs(p) <= M * EPS<T>::eps() * sqrt(d[i] * d[j])) {
                    continue;
                }

                T y  = d[i] - d[j];
                T r  = hypot(p * 2, y);
                T r2 = r * 2;
                T c, s;
                if (y >= 0) {
                    c = sqrt((r + y) / r2);
                    s = p / (r2 * c);
                } else {
                    s = sqrt((r - y) / r2);
                    c = p / (r2 * s);
                }

                T a = 0, b = 0;
                for (int k = 0; k < M; k++) {
                    T t0  = c * Si[k] + s * Sj[k];
                    T t1  = c * Sj[k] - s * Si[k];
                    Si[k] = t0;
                    Sj[k] = t1;

                    a += t0 * t0;
                    b += t1 * t1;
                }
                d[i] = a;
                d[j] = b;

                for (int l = 0; l < N; l++) {