Avoiding transpose

[MichelJuillard]

1. Transposing a matrix is an expensive operation that should be avoided if at all possible
2. The Tasmanian C library uses only vectors as input or output arguments
3. In Julia, one can pass a matrix to ccall vector argument. The vector is just the continuous memory space occupied by the Julia matrix.
4. The Tasmanian vector is formed by a succession of strips of data piled on top of each others. As a Julia matrix is stored column wise, the most efficient is to have the Julia matrix columns to correspond to the strips of data

Example:

# load needed points
"""
    loadNeededPoints!(tsg::TasmanianSG, vals::Matrix{Float64})

loads the values of the target function at the needed points
if there are no needed points, this reset the currently loaded
values

vals: a matrix
      with dimensions iOutputs X getNumNeeded() 
      each column corresponds to the values of the outputs at
      the corresponding needed point. The order and leading
      dimension must match the points obtained from
      getNeededPoints()
"""
function loadNeededPoints!(tsg::TasmanianSG,vals::Array{Float64})
    if ndims(vals) != 2
        throw(ArgumentError("vals should be a 2-D array, instead it has $(ndims(vals)) dimensions"))
    end
    n = size(vals)
    if n[1] != getNumOutputs(tsg)
        throw(ArgumentError("leading dimension of vals is $(n[1]) but the number of outputs is set to $(getNumOutputs(tsg))"))
    end
    if getNumNeeded(tsg) == 0
        if n[2] != getNumLoaded(tsg)
            throw(ArgumentError("the second dimension of vals is $(n[2]) but the number of current points is $(getNumLoaded(tsg))"))
        end
    elseif n[2] != getNumNeeded(tsg)
        throw(ArgumentError("the second dimension of vals is $(n[2]) but the number of needed points is $(getNumNeeded(tsg))"))
    end
    ccall((:tsgLoadNeededPoints,TASlib),Nothing,(Ptr{Nothing},Ptr{Cdouble}),tsg.pGrid,vals)
end

[floswald]

Ok I see, that starts to ring a bell, yes. One thing that would need to be clarified maybe is how to deal with more than 2D. Like what if you have a 3D object. then one needs to reshape ones array into 2D first. correct? Also, in which way is the current implementation "wrong" here - I can't see any transpose operation:

https://github.com/floswald/Tasmanian.jl/blob/81ce15844119e7cd3437078163c77de1b3e64308/src/TSG.jl#L124C5-L124C5

[MichelJuillard]

I haven't looked at the use of more than 2D objects in Tasmanian, but I see they have multi-index support and hierarchical grids. We probably shouldn't make assumption about the dimensions of the arrays that Tasmanian.jl accept and leave it to the user to prepare that memory block correctly. In the current implementation, you have transpose on lines 106, 118 and 164 Best Michel

…

________________________________ From: Florian Oswald ***@***.***> Sent: Tuesday, November 14, 2023 8:49 AM To: floswald/Tasmanian.jl Cc: Michel Juillard; Author Subject: Re: [floswald/Tasmanian.jl] Avoiding transpose (Issue #7) Ok I see, that starts to ring a bell, yes. One thing that would need to be clarified maybe is how to deal with more than 2D. Like what if you have a 3D object. then one needs to reshape ones array into 2D first. correct? Also, in which way is the current implementation "wrong" here - I can't see any transpose operation: https://github.com/floswald/Tasmanian.jl/blob/81ce15844119e7cd3437078163c77de1b3e64308/src/TSG.jl#L124C5-L124C5 — Reply to this email directly, view it on GitHub<#7 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAEOULGG2TBS34UZSXNT67TYEMO7JAVCNFSM6AAAAAA7JYSAMSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMBZGY4TGMBZGA>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[floswald]

ah yes, I see. I was looking at loadNeededPoints! only. indeed i'm seeing this one

Tasmanian.jl/src/TSG.jl

Line 164 in 81ce158

v = Array(reshape(vals',nx*nd,1))

which seems hard to avoid. in the end you always need to strip away all dimension info from the array and supply a 1D vector to tasmanian.

[MichelJuillard]

Not really, in Julia you pass to ccall()l a pointer to a Julia array with any number of dimensions. Tasmanian C library will reinterpret it as a vector. So what matters is only the order of the elements in the continuous memory block containing the Julia array. For evaluateBatch() we can avoid the transpose by requesting that the columns and rows of the vals matrix be exchanged Best Michel

…

________________________________ From: Florian Oswald ***@***.***> Sent: Tuesday, November 14, 2023 9:25 AM To: floswald/Tasmanian.jl Cc: Michel Juillard; Author Subject: Re: [floswald/Tasmanian.jl] Avoiding transpose (Issue #7) ah yes, I see. I was looking at loadNeededPoints! only. indeed i'm seeing this one https://github.com/floswald/Tasmanian.jl/blob/81ce15844119e7cd3437078163c77de1b3e64308/src/TSG.jl#L164 which seems hard to avoid. in the end you always need to strip away all dimension info from the array and supply a 1D vector to tasmanian. — Reply to this email directly, view it on GitHub<#7 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAEOULCYWOBZLV3HEYCSERLYEMTHNAVCNFSM6AAAAAA7JYSAMSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMBZG4ZTQMBUHE>. You are receiving this because you authored the thread.Message ID: ***@***.***>