Extending inner products

[cortner]

Two related questions:

1. I have the need for a generalised euclidean inner product

   dot(x, M, y)  =  dot(x, M * y)

which I don't want to hard-code prefer an ability to overload it. Is this functionality sufficiently general (I think it is) that it should be in Base by default for the standard arrays?

EDIT: second question is deleted; it needs more thinking.

[andreasnoack]

I that this would be useful to have in base and I've thought about adding it several times. One of the benefits is that you can completely avoid a temporary M*y.

[cortner]

yes - that was exactly the origin of this - I should have said, sorry.

[stevengj]

Would be nice to get some benefit from BLAS here, though it's only a BLAS2 operation. @xianyi, do you have a suggestion on the most efficient way to use BLAS for x'*M*y?

One option would be to break x into chunks of say, 16 elements (big enough to exceed a cache line but small enough that the memory allocation of a length-16 temporary vector is an issue, and compute x'*M on that chunk via BLAS2 using a subset of the columns of M (since Julia is column-major, this doesn't require any copying), then accumulate that chunk of the dot product with y. However, this is suboptimal because it doesn't fully exploit multi-threading: each of the chunk dot products could ideally be computed in a parallel thread.

[stevengj]

Ideally you'd also want a specialized method for the common case where M is a diagonal matrix, especially a diagonal real matrix.

And tridiagonal "mass" matrices are also common, as well as general sparse matrices M, in problems arising from discretized PDEs.

[cortner]

I was mostly thinking of sparse cases for my own use. These should be relatively easy to implement in Julia?

[andreasnoack]

For the matvec, it seems that multithreading is almost the only speedup you'll get from BLAS. For a 500x500 problem and a Julia double loop with @inbounds and @simd I get

julia> blas_set_num_threads(1);

julia> @time for i = 1:10000; Ac_mul_B!(y, A, x); end # OpenBLAS 1 thread
  0.484342 seconds

julia> blas_set_num_threads(4);

julia> @time for i = 1:10000; Ac_mul_B!(y, A, x); end # OpenBLAS 4 threads
  0.184681 seconds

julia> @time for i = 1:10000; mymul(y, A, x); end # Julia
  0.509110 seconds

so if/when the @threadsoverhead is gone, a pure Julia version should probably be okay.

[stevengj]

I'm surprised there isn't some SIMD speedup.

[mschauer]

Probably also worthwhile to special case x === y and issym(A).

[andreasnoack]

@stevengj I don't understand. Both OpenBLAS and Julia use SIMD so where would you have expected to have seen the speedup?

[stevengj]

In OpenBLAS the SIMD is hard-coded (I thought), whereas in Julia we rely on the compiler to do its magic, which doesn't always work as well. Maybe this is one of the cases where the compiler does a good job in exploiting SIMD.

[JaredCrean2]

In my experiments with mat-vec and mat-mat, Julia has done a good job exploiting SIMD. For mat-vec my implemention was asymptotically as fast as OpenBlas, although they use a better algorithm for mat-mat, so theirs was faster in that case (about about 3x IIRC).

[simonbyrne]

To follow up on this, for when everything is sparse we can also do much better by interleaving the computations. e.g. in a recent project I did:

function mul_ut_X_v(u::SparseVector, X::SparseMatrixCSC, v::SparseVector)
    u.n == X.m && X.n == u.n || throw(DimensionMismatch())
    
    z = zero(promote_type(eltype(u),eltype(X),eltype(v)))
    for (vi,vv) in zip(v.nzind, v.nzval)        
        X_ptr_lo = X.colptr[vi] 
        X_ptr_hi = X.colptr[vi+1] - 1
        if X_ptr_lo <= X_ptr_hi        
            z += Base.SparseArrays._spdot(dot, 1, length(u.nzind), u.nzind, u.nzval,
                                            X_ptr_lo, X_ptr_hi, X.rowval, X.nzval) * vv
        end
    end
    z
end

Once we have lazy transposes, we could overload the ternary multiplication operator u' * X * v for this. Until then, I would be fine with overloading dot, or adding a bilinear function.