cholesky_update.h

Go to the documentation of this file.

#ifndef FUSED_ALGORITHMS_CHOLESKY_UPDATE_H
#define FUSED_ALGORITHMS_CHOLESKY_UPDATE_H
 
#include "config.h"
#include "fused/fused_matrix.h"
#include "fused/fused_vector.h"
#include "math/math_utils.h" // math::sqrt
 
namespace matrix_algorithms
{
 
    template <typename T, my_size_t N>
    FusedMatrix<T, N, N> cholesky_rank1_update(
        const FusedMatrix<T, N, N> &L,
        const FusedVector<T, N> &v)
    {
        static_assert(is_floating_point_v<T>,
                      "cholesky_rank1_update requires a floating-point scalar type");
 
        FusedMatrix<T, N, N> Lp = L; // work on a copy
        FusedVector<T, N> p = v;     // work vector
 
        for (my_size_t k = 0; k < N; ++k)
        {
            T r = math::sqrt(Lp(k, k) * Lp(k, k) + p(k) * p(k));
            T c = Lp(k, k) / r;
            T s = p(k) / r;
 
            Lp(k, k) = r;
 
            for (my_size_t i = k + 1; i < N; ++i)
            {
                T tmp = Lp(i, k);
                Lp(i, k) = c * tmp + s * p(i);
                p(i) = c * p(i) - s * tmp;
            }
        }
 
        return Lp;
    }
 
} // namespace matrix_algorithms
 
#endif // FUSED_ALGORITHMS_CHOLESKY_UPDATE_H