WIP: cleanup blas.jl scal!

[nwh]

This pull request addresses the issue discussed in #8415.

1. scal! is no longer a method for StridedArrays
2. helper methods are added for StridedVectors
3. an issue with scal! when scaling a complex array with a real number is uncovered
   the technique to save flops assumes incx == 1
   a check for this is added, however this may not be a desirable solution
4. another solution to issue 3 above is to use 2 calls to the appropriate blas function. The method could then be extended to StridedVectors

[nwh]

Interesting. The build for osx passed, while the build for linux failed.

[nwh]

The issue with scal! with real scalar and complex array is described
with the following example:

using Base.LinAlg.BLAS
n = 5
A = ones(n) + im*ones(n)
scal!(3,0.0,A,2)

Result:

5-element Array{Complex{Float64},1}:
 0.0+1.0im
 0.0+1.0im
 0.0+1.0im
 0.0+1.0im
 0.0+1.0im

This is also an out of bounds memory access.

[andreasnoack]

Isn't this one ready for merging? The Travis test error is not related to your changes.

[nwh]

Yes, it can be merged. I thought the solution for the special case of real scalar, complex vector might not be satisfactory. Anyway, it is a step in the right direction.

[andreasnoack]

Okay. I can see what you mean now. Good catch.

Please have a look at this method which does something similar. It might be better to do something similar for the more high levelscale instead of BLAS.scal! and avoid such tricks in blas.jl. What do you think?

[nwh]

I think it would be a good idea to move the real-complex scal! method out of blas.jl.

The A_mul_B! method you mention also has a problem. It does not seem possible to reinterpret a SubArray:

A = rand(Complex128,10,10)
B = sub(A,1:2:9,1:2:9)
C = reinterpret(Float64,B,(2size(B,1),size(B,2)))

Result:

ERROR: `reinterpret` has no method matching reinterpret(::Type{Float64}, ::SubArray{Complex{Float64},2,Array{Complex{Float64},2},(StepRange{Int64,Int64},StepRange{Int64,Int64})}, ::(Int64,Int64))

[andreasnoack]

Good point. You could consider to define reinterpret for this operation, but that operation would have to change the range parameters of the new object an maybe that is too much side effect for reinterpret. A more transparent solution could be to define something like convert(SubMatrix{Float64}, A::StridedMatrix{Complex128}) methods for this operation.

[timholy]

Can't even do it at all with current SubArrays, because in general it would be "non-strided": if your step along the first dimension were 4, then when reinterpreting as Float64 you'd have to take 2, then skip 6, take another 2, then skip 6, etc.

However, this is doable with the implementation in #8501.

[timholy]

You could do it with current SubArrays if you increased the dimensionality, though, inserting a new first dimension of size 2.

[andreasnoack]

The requirement that stride(A,1)=1 is already in @nwh's proposal. BLAS can only handle contiguous block for matrices and regular step sizes for vectors so I'd be fine with convert throwing an error when the first stride is not one.

[timholy]

Sorry, I didn't read the whole thing; I was just looking at his B = sub(A,1:2:9,1:2:9) call in the post two above mine, and replying to whether reinterpret could be defined on SubArrays.

[andreasnoack]

I see. Then you can get a sorry back from here. I didn't see that @nwh had step sizes of two in his example. That wouldn't be possible to handle for BLAS.

[timholy]

Indeed. Our generic matrix multiplication could handle it, but you're right it's problematic for BLAS.

I suspect this is going to become a recurring theme with our explorations on returning views. There are a lot of operations (like reinterpret and reshape) that we take for granted with Arrays. #8501 is flexible enough to implement all of those even for views, but sometimes the cost (like in this example) will be to transition to a non-strided representation. In such cases we have to decide whether to call a pure-julia fallback or make a copy; either choice has its downsides. The performance gap between BLAS and julia is narrowing rapidly (and will narrow further once we have threading), but LAPACK is huge.

[nwh]

For the case of scal! in blas.jl: it might be easiest to just remove (or deprecate) the specialized method for real-scalar, complex-array. I agree with @andreasnoack, this behavior would be a better fit for a higher level function.

[ViralBShah]

Yes, simpler to remove it, since it is not natively supported by BLAS.

[ViralBShah]

Am I right in saying that it is good to merge this after removing support for the real-scalar complex-array case?

[nwh]

Yes, I agree that real-scalar complex-array case should be removed. I was going to make a new PR on the updated blas.jl (with the 64_ suffix to openblas calls).

[ViralBShah]

Thanks! Closing this one in that case, and awaiting the new one.