Static version of lowrankdowndate

staticlowrankdowndate.jl

# Static version of lowrankdowndate.
# The source code is derived from stdlib/LinearAlgebra/src/cholesky.jl.
# License is MIT (https://julialang.org/license)

using StaticArrays: SMatrix
using LinearAlgebra: LowerTriangular, PosDefException

"""
    lowrankdowndate(L, v)

Compute the cholesky factor of `L*L' - v*v'`, where `L` is a static matrix
wrapped by `LowerTriangular`.

This function does not allocate heap memory and therefore is faster than
`LinearAlgebra.lowrankdowndate` for small matrices.
"""
@generated function lowrankdowndate(L::LowerTriangular{T,<:SMatrix{M,M,T}}, v::AbstractVector) where {T,M}
    q = quote
        $M == length(v) ||
            throw(DimensionMismatch("updating vector must fit size of factorization"))
        A = L.data
    end
    emit(ex) = push!(q.args, ex)
    for i in 1:M
        vi = Symbol(:v, i)
        emit(:($vi = v[$i]))
        for j in i:M
            Aji = Symbol(:A, j, i)
            emit(:($Aji = A[$j,$i]))
        end
    end
    for i in 1:M
        vi = Symbol(:v, i)
        Aii = Symbol(:A, i, i)
        emit(quote
            s = conj($vi / $Aii)
            s2 = abs2(s)
            s2 > 1 && throw(PosDefException($i))
            c = sqrt(1 - s2)
            $Aii = c*$Aii
        end)
        for j in i+1:M
            vj = Symbol(:v, j)
            Aji = Symbol(:A, j, i)
            emit(quote
                tmp = ($Aji - s*$vj)/c
                $Aji = tmp
                $vj = -s'tmp + c*$vj
            end)
        end
    end
    elms = [i ≥ j ? Symbol(:A, i, j) : zero(T) for j in 1:M for i in 1:M]
    emit(:(return SMatrix{$M,$M}($(elms...))))
    return Expr(:block, Expr(:meta, :inline), q)
end

# quick test
#using StaticArrays: @SMatrix
#using LinearAlgebra: cholesky
#L = LowerTriangular(@SMatrix [1.0 0.0; 2.0 3.0])
#v = [0.5, 1.0]
#@assert lowrankdowndate(L, v) ≈ cholesky(Matrix(L*L') - v*v').L