Porting sve.jl to C++

[terasakisatoshi]

https://github.com/SpM-lab/SparseIR.jl/blob/main/src/sve.jl

[terasakisatoshi]

I'm working on it.

[terasakisatoshi]

#36 partially resolves

[terasakisatoshi]

We could write postprocess(sve::AbstractSVE, u, s, v) without using complicated broadcasting techniques.

"""
    postprocess(sve::AbstractSVE, u, s, v)

Construct the SVE result from the SVD.
"""
function postprocess(sve::SamplingSVE, (u,), (s,), (v,))
    s = Float64.(s)
    u_x = u ./ sqrt.(sve.gauss_x.w)
    v_y = v ./ sqrt.(sve.gauss_y.w)

    u_x3d = reshape(u_x, (sve.n_gauss, length(sve.segs_x) - 1, length(s)))
    v_y3d = reshape(v_y, (sve.n_gauss, length(sve.segs_y) - 1, length(s)))

    cmat = legendre_collocation(sve.rule)
    u_data = reshape(cmat * reshape(u_x3d, (size(u_x3d, 1), :)),
        (:, size(u_x3d, 2), size(u_x3d, 3)))
    v_data = reshape(cmat * reshape(v_y3d, (size(v_y3d, 1), :)),
        (:, size(v_y3d, 2), size(v_y3d, 3)))

    dsegs_x = diff(sve.segs_x)
    dsegs_y = diff(sve.segs_y)
    # Using nested for loops to multiply u_data
    for j in 1:size(u_data, 2)
        for k in 1:size(u_data, 3)
            for i in 1:size(u_data, 1)
                u_data[i, j, k] *= sqrt(0.5 * dsegs_x[j])
            end
        end
    end

    # Using nested for loops to multiply v_data
    for j in 1:size(v_data, 2)
        for k in 1:size(v_data, 3)
            for i in 1:size(v_data, 1)
                v_data[i, j, k] *= sqrt(0.5 * dsegs_y[j])
            end
        end
    end
    # Construct polynomials
    ulx = PiecewiseLegendrePolyVector(Float64.(u_data), Float64.(sve.segs_x))
    vly = PiecewiseLegendrePolyVector(Float64.(v_data), Float64.(sve.segs_y))
    canonicalize!(ulx, vly)
    return SVEResult(ulx, s, vly, sve.kernel, sve.ε)
end

[terasakisatoshi]

To port the postprocess function to C++11, I rewrite it without using complex matrix multiplication and broadcasting:

"""
    postprocess(sve::AbstractSVE, u, s, v)

Construct the SVE result from the SVD.
"""
function postprocess(sve::SamplingSVE, (u,), (s,), (v,))
    s = Float64.(s)
    u_x = u ./ sqrt.(sve.gauss_x.w)
    v_y = v ./ sqrt.(sve.gauss_y.w)

    u_x3d = reshape(u_x, (sve.n_gauss, length(sve.segs_x) - 1, length(s)))
    v_y3d = reshape(v_y, (sve.n_gauss, length(sve.segs_y) - 1, length(s)))

    cmat = legendre_collocation(sve.rule)
    T = eltype(cmat)
    # Using for loop instead of matrix multiplication for u_data
    u_data = Array{T}(undef, size(cmat, 1), size(u_x3d, 2), size(u_x3d, 3))
    for j in 1:size(u_x3d, 2)
        for k in 1:size(u_x3d, 3)
            for i in 1:size(cmat, 1)
                u_data[i, j, k] = zero(T)
                for l in 1:size(cmat, 2)
                    u_data[i, j, k] += cmat[i, l] * u_x3d[l, j, k]
                end
            end
        end
    end

    # Using for loop instead of matrix multiplication for v_data
    v_data = Array{T}(undef, size(cmat, 1), size(v_y3d, 2), size(v_y3d, 3))
    for j in 1:size(v_y3d, 2)
        for k in 1:size(v_y3d, 3)
            for i in 1:size(cmat, 1)
                v_data[i, j, k] = zero(T)
                for l in 1:size(cmat, 2)
                    v_data[i, j, k] += cmat[i, l] * v_y3d[l, j, k]
                end
            end
        end
    end

    # Replace diff with a for loop for segs_x
    dsegs_x = Vector{T}(undef, length(sve.segs_x) - 1)
    for i in 1:length(dsegs_x)
        dsegs_x[i] = sve.segs_x[i + 1] - sve.segs_x[i]
    end

    # Replace diff with a for loop for segs_y
    dsegs_y = Vector{T}(undef, length(sve.segs_y) - 1)
    for i in 1:length(dsegs_y)
        dsegs_y[i] = sve.segs_y[i + 1] - sve.segs_y[i]
    end

    # Using nested for loops to multiply u_data
    for j in 1:size(u_data, 2)
        for k in 1:size(u_data, 3)
            for i in 1:size(u_data, 1)
                u_data[i, j, k] *= sqrt(0.5 * dsegs_x[j])
            end
        end
    end

    # Using nested for loops to multiply v_data
    for j in 1:size(v_data, 2)
        for k in 1:size(v_data, 3)
            for i in 1:size(v_data, 1)
                v_data[i, j, k] *= sqrt(0.5 * dsegs_y[j])
            end
        end
    end
    # Construct polynomials
    ulx = PiecewiseLegendrePolyVector(Float64.(u_data), Float64.(sve.segs_x))
    vly = PiecewiseLegendrePolyVector(Float64.(v_data), Float64.(sve.segs_y))
    canonicalize!(ulx, vly)
    return SVEResult(ulx, s, vly, sve.kernel, sve.ε)
end