Specifying the numerical precision and ranks of the array arguments of a function

New to Julia
1

Following is a Fortran function for matrix vector multiplication. It specifies that a(:,:) is a rank-2 array of 64-bit floats and that x(:) and b(:) are rank-1 arrays of 64-bit floats.

function matrix_vec_mult(a,x) result(b)
! return the product of matrix a and vector x
use iso_fortran_env, only: real64
real(kind=real64), intent(in) :: a(:,:)
real(kind=real64), intent(in) :: x(:)
real(kind=real64)             :: b(size(a,1))
integer                       :: i
do i=1,size(a,1)
   b(i) = sum(a(i,:)*x)
end do
end function matrix_vec_mult

I know that arguments can be declared Array{Number} in Julia, but is it possible to be as specific as the code above regarding argument rank and numerical precision?

The intent(in) designation means that the argument cannot be changed within the procedure. Is there a Julia analog for that?

2

Yes, that would be Array{Float64, 2} and Array{Float64, 1} respectively (or the aliases Matrix{Float64} and Vector{Float64}).

No, this isn’t really possible in Julia. Immutable structs can help organize data and control mutability, but there’s nothing you can do in a function definition to enforce that the function won’t mutate some mutable input (other than documenting that fact).

5 Likes
3

Note that while you can specify what types your functions take, there is not a performance advantage in doing so.

function matrix_vec_mult(a,x)
    b = Vector{eltype(a)}(undef, size(a,1))
    for i=1:size(a,1)
        b[i] = sum(a[i,:]*x)
    end
    return b
end

will have the exact same performance as

function matrix_vec_mult(a::Matrix{Float64},x::Vector{Float64})
    b = Vector{Float64}(undef, size(a,1))
    for i=1:size(a,1)
        b[i] = sum(a[i,:]*x)
    end
    return b
end

because Julia always compiles a specialized version of your code for the argument types provided.

5 Likes
4

Note that this will give an error in Julia, because a 1d slice like a[i,:] creates a 1d array, not a row vector.

You can use sum(a[i,:] .* x), though this is suboptimal because it creates a temporary array before summing. To avoid that temporary array, you can use dot(a[i,:], x) from the LinearAlgebra package (noting that dot conjugates the first argument if it is complex), or a[i,:]' * x (same as dot) or transpose(a[i,:]) * x (which doesn’t conjugate). a[i,:] also creates a copy for the slice, but you can avoid that by putting @views (e.g. before function).

(Of course, in reality you would just use b = a * x for a matrix–vector multiplication.)

3 Likes