- plan the following refactor of exa_linalg:
- rename Vec1DReal to VecReal (a Vec is a 1D object bloody hell!), Vec1DNode to VecNode, Mat2D to Mat
const VecReal{T<:Real} = Union{Vector{T}, DenseSubArray{T,1}, DenseReshapedArray{T,1}, Base.ReinterpretArray{T,1}}
const VecNode{T<:ExaModels.AbstractNode} = Union{Vector{T}, DenseSubArray{T,1}, DenseReshapedArray{T,1}}
const Mat{T} = Union{Matrix{T}, DenseSubArray{T,2}, DenseReshapedArray{T,2}}
- for the recall, these types are used to extend standard linear algebra functions to combinations of Real and AbstractNode types without creating any ambiguity
- rule no. 1: as for * (which is extended and optimised to implement 0 * e = 0 and used everywhere possible as a building block), use an optimised __dot function to define ALL operations; here is the code of the optimised __dot function:
function __dot(v::AbstractVector{T}, w::AbstractVector{S}) where {T<:Real, S<:ExaModels.AbstractNode}
@assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
return sum(ExaModels.Null(v[i]) * w[i] for i in eachindex(v))
end
function __dot(v::AbstractVector{T}, w::AbstractVector{S}) where {T<:ExaModels.AbstractNode, S<:Real}
@assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
return sum(v[i] * ExaModels.Null(w[i]) for i in eachindex(v))
end
function __dot(v::AbstractVector{T}, w::AbstractVector{S}) where {T<:ExaModels.AbstractNode, S<:ExaModels.AbstractNode}
@assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
return sum(v[i] * w[i] for i in eachindex(v))
end
- this function is in particular used to redefine dot for our specific types (not on AbstractVector to prevent type ambiguities!):
dot(v::VecReal{T}, w::VecNode{S}) where {T<:Real, S<:ExaModels.AbstractNode} = __dot(v, w)
dot(v::VecNode{T}, w::VecReal{S}) where {T<:ExaModels.AbstractNode, S<:Real} = __dot(v, w)
dot(v::VecNode{T}, w::VecNode{S}) where {T<:ExaModels.AbstractNode, S<:ExaModels.AbstractNode} = __dot(v, w)
- rule no. 2: do not use slices like A[:, j] that make copies (allocations!) but ALWAYS views instead, like view(A, :, j) (never the macro @view)
- if necessary, consider adjoint x vector (currently there is only adjoint x matrix)