Argument order for rand!

[simonster]

Is there a reason that the output array comes last instead of first as is the case for most other mutating functions? Right now we could change the argument order and add a deprecation, but we should decide whether we want to change this before #6003.

[andreasnoack]

This has been discussed somewhere before and I think @johnmyleswhite gave the answer. Personally, I don't like the rule that the mutated argument should be the first one. It could be awkward for some functions and for all the matrix multiplication and division functions I'd prefer that we followed the BLAS conventions strictly and that would mean that the mutating argument was last (as well as adding scalar α and β arguments). For some linear algebra functions you'd also want to give a workspace array as argument to avoid reallocation and such an argument would make more sense as the last argument.

[johnmyleswhite]

I'm happy with any rule personally. Whether the mutated arguments are first or last doesn't matter much to me, although the argument against putting mutated arguments last is that it breaks functions with varargs.

[simonster]

I don't have a preference on the rule, just that there should be one, or that we should get #249. But I also think at this point we don't want to break all code that uses mutating functions by moving the output last.

As far as matrix functions, maybe α, β, and workspace should be keyword arguments? Although right now that incurs kwsorter overhead.

[iamed2]

Perhaps there could be a convention where the names of mutating arguments end in an exclamation point?

rand!{T}(r::Range{T},A::AbstractArray{T,N})

would become

rand!{T}(r::Range{T},A!::AbstractArray{T,N})

in docs and when calling methods().

[timholy]

I don't think there are any established conventions, here or elsewhere. While BLAS puts it last, more often than not libc puts them first (memcpy, fgets, strcpy). But even libc is inconsistent (sscanf, although there's a reason for that violation). I have a slight preference for first, but above all would prefer reasonable consistency (whatever that may be)

I would also conceptually distinguish output from mutating. mutating to me means it's OK to destroy it. output means just that, it's acting as a return value. As discussed elsewhere, there are sometimes performance advantages in passing in a dedicated output, even if you could just use the input to return the output. This is one of the things I definitely don't like about the BLAS/LAPACK API.

[timholy]

@iamed2, it's a good suggestion, but keyword arguments have an amusing gotcha on this: #8020.

[timholy]

Better link for the reason for separating output from input: #7513 (comment)

[andreasnoack]

@timholy I can't see the performance advantage for triangular solve on my machine

julia> minimum([@elapsed (for i = 1:100;MyTest.mysolve1(A, b, x);end) for j = 1:10])
0.306543869
julia> minimum([@elapsed (for i = 1:100;copy!(x,b);MyTest.mysolve2(A, x);end) for j = 1:10])
0.305641564

I'm not that surprised about the result. Efficiency seems to have been a pretty high priority in the development of BLAS.

The code is

module MyTest

    function mysolve1{T,S}(A::Triangular{T,S,:L,false}, b, x = b)
        n = size(A, 1)
        Ad = A.data
        @inbounds begin
            for j = 1:n
                xj = b[j]
                for i = 1:j-1
                    xj -= Ad[j,i]*x[i]
                end
                x[j] = xj/A[j,j]
            end
        end
        x
    end

    function mysolve2{T,S}(A::Triangular{T,S,:L,false}, b)
        n = size(A, 1)
        Ad = A.data
        @inbounds begin
            for j = 1:n
                bj = b[j]
                for i = 1:j-1
                    bj -= Ad[j,i]*b[i]
                end
                b[j] = bj/A[j,j]
            end
        end
        b
    end
end

[andreasnoack]

@simonster I really think it is a mistake that we haven't followed the BLAS convention completely and I also hope that we can still break code to make things right.

Actually, I began the work to change the convention some time ago, but the development stalled in discussion.

[timholy]

Sure, for this kind of algorithm it won't matter. Compare the running time to that of copy!ing x---for this algorithm it's access to the elements of A, not x, that matters.

But that doesn't mean that it's a good model to apply to other algorithms that may be cache-limited on the output.

[rfourquet]

All rand methods are of the form rand([rng], [out]), with:

• rng is the the generator, and defaults to the global one
• out is the type information of the output:
  • first the type of randoms, which defaults to Float64
  • and the dimensions of the generated array, or nothing for a scalar

Note that if the rng is given, the type of randoms can not be specified, it is assumed that it's the responsibility of rng.

So it makes sense that an output array A comes last in rand!. Using @timholy's terminology (#8246 (comment)), A can be seen as a "mutating" argument which provides both the type and dimensions (out above) to the rand! function. It also happens to be an "output" argument. Hence the signature: rand!([rng], out).

[rfourquet]

To insist on the "output" characteristic of A, why not A[:] = randgen([rng], T) where randgen([rng]]) an infinite iterator producing random T's (or maybe with a bit of magic in the relevant getindex to deduce automatically T=eltype(A)) ?

[simonster]

I suppose there is also some analogy to fill!, which takes the array last.