Add rot to BLAS in stdlib/LinearAlgebra

[amontoison]

Applies the plane rotation :

$$[ c s] [x] [-s̅ c] [y]$$

[andreasnoack]

I don't really mind adding this but notice that the functionality is already covered by LinearAlgebra.Givens, i.e. you can do

julia> LinearAlgebra.Givens(1, 3, 0.3, 0.5)*[1, 1, 1]
3-element Array{Float64,1}:
  0.8
  1.0
 -0.2

julia> lmul!(LinearAlgebra.Givens(1, 3, 0.3, 0.5), [1.0, 1, 1])
3-element Array{Float64,1}:
  0.8
  1.0
 -0.2

and it's probably more efficient than calling out to BLAS.

[StefanKarpinski]

Maybe worth a speed comparison then?

[amontoison]

@andreasnoack x and y are vectors here, can we do the same thing with LinearAlgebra.Givens?
An example of use is inside an implementation of MINRES-QLP (Krylov solver for linear systems):

# [ẘₖ₋₂ w̄ₖ₋₁ vₖ] [cpₖ  0   spₖ] [1   0    0 ] = [wₖ₋₂ ẘₖ₋₁ w̄ₖ] ⟷ wₖ₋₂ = cpₖ * ẘₖ₋₂ + spₖ * vₖ
#                [ 0   1    0 ] [0  cdₖ  sdₖ]                  ⟷ ẘₖ₋₁ = cdₖ * w̄ₖ₋₁ + sdₖ * (spₖ * ẘₖ₋₂ - cpₖ * vₖ)
#                [spₖ  0  -cpₖ] [0  sdₖ -cdₖ]                  ⟷ w̄ₖ   = sdₖ * w̄ₖ₋₁ - cdₖ * (spₖ * ẘₖ₋₂ - cpₖ * vₖ)

# Compute wₐᵤₓ = spₖ * ẘₖ₋₂ - cpₖ * vₖ
axpby!(n, -cpₖ, vₖ, spₖ, ẘₖ₋₂)
wₐᵤₓ = ẘₖ₋₂
# Compute ẘₖ₋₁ and w̄ₖ
ref!(n, w̄ₖ₋₁, wₐᵤₓ, cdₖ, sdₖ)

with

ref!(n, x, y, c, s) = rot!(n, y, x, s, c)

I update ẘₖ₋₁ and w̄ₖ without allocating extra memory. Does LinearAlgebra.Givens can do that without allocating?

And in general we like to use Givens reflections (symmetric) instead of Givens rotations.

[dkarrasch]

I'd suggest to include testing that x and y are actually overwritten, something along the lines

Suggested change

	x2, y2 = BLAS.rot!(n,copy(x),1,copy(y),1,c,s)
	@test x2 ≈ c*x + s*y
	@test y2 ≈ -s*x + c*y
	x2 = copy(x)
	y2 = copy(y)
	BLAS.rot!(n, x, 1, y, 1, c, s)
	@test x ≈ c*x2 + s*y2
	@test y ≈ -s*x2 + c*y2

and similarly below.

[andreasnoack]

The current methods for Givens might not do exactly what you want since it's written for the situation where x and y are stored in the same array. Your example doesn't render correctly so I'm not sure exactly what you are doing. In any case, a version that works for two vectors instead of two rows in a matrix should be fairly simple to write in Julia. I still think it could easily be as fast or faster than the BLAS.

[dkarrasch]

Indeed, I found and adopted this one

julia/stdlib/LinearAlgebra/src/givens.jl

Lines 339 to 351 in 2a5bb59

	@inline function lmul!(G::Givens, A::AbstractVecOrMat)
	require_one_based_indexing(A)
	m, n = size(A, 1), size(A, 2)
	if G.i2 > m
	throw(DimensionMismatch("column indices for rotation are outside the matrix"))
	end
	@inbounds for i = 1:n
	a1, a2 = A[G.i1,i], A[G.i2,i]
	A[G.i1,i] = G.c *a1 + G.s*a2
	A[G.i2,i] = -conj(G.s)*a1 + G.c*a2
	end
	return A
	end

to this:

function rot!(x::AbstractVector, y::AbstractVector, c, s)
    LinearAlgebra.require_one_based_indexing(x, y)
    n = length(x)
    n == length(y) || throw(DimensionMismatch("different lenghts"))
    @inbounds for i = 1:n
        xi, yi = x[i], y[i]
        x[i] =       c *xi + s*yi
        y[i] = -conj(s)*xi + c*yi
    end
    return x, y
end

which passes the tests. This is quick-and-dirty, one could be more permissive in terms of lengths (like pass an extra n and apply the rotation to the first n elements only, or determine n as the minimum of the lengths of the vectors) etc. @amontoison Want to give it a shot? It has one allocation that I don't see. Here is a little comparison with how one would obtain the same result with Givens:

using BenchmarkTools, Test
c, s = rand(), rand()
x0, y0 = rand(1000), rand(1000);
@btime rot!(x, y, $c, $s) setup=(x = copy(x0); y = copy(y0)); # 245 ns
@btime lmul!(LinearAlgebra.Givens(1, 2, $c, $s), A) setup=(A = [reshape(x0, 1, :); reshape(y0, 1, :)]); # 721 ns
x, y = rand(1000), rand(1000);
A = [reshape(x, 1, :); reshape(y, 1, :)];
rot!(x, y, c, s)
lmul!(LinearAlgebra.Givens(1, 2, c, s), A)
@test A ≈ [reshape(x, 1, :); reshape(y, 1, :)]

[amontoison]

Sorry for the delay, I did some benchmarks @dkarrasch and your rot! performs as well as the BLAS version for me. I decided to add it in generic.jl with a variant ref! that applies the reflection instead of the rotation.

The BLAS version could be still relevant if you use mkl instead of openblas for example.

About the allocation, it's the pointers for x and y that are returned at the end 😉

[dkarrasch]

rot! -> ref!

[dkarrasch]

Suggested change

	rot!(n, x, y, c, s)
	rot!(x, y, c, s)

[dkarrasch]

Suggested change

	ref!(n, x, y, c, s)
	ref!(x, y, c, s)

[andreasnoack]

Great with generic implementations. Since they are no longer wrappers of BLAS functions, I'll suggest that the names are made more informative since the F77 variable name restrictions don't apply. So maybe rotate! and reflect!?

[dkarrasch]

LGTM

[StefanKarpinski]

This seems good. @dkarrasch, anyone else you'd want to have look it over?