multishift qr and blas3

01_example.txt

Chasing two double-shift-bulges one step forward using two 
reflections G1 and G2 of size 3 (they are each composed of 
two Givens rotations).

x x x x x x x x                         x x x x x x x x
x x x x x x x x ┐                       x x x x x x x x
x x x x x x x x │ double shift G2       . x x x x x x x
x x x x x x x x ┘                       . x x x x x x x
. . . x x x x x ┐                       . x x x x x x x
. . . x x x x x │ double shift G1       . . . . x x x x
. . . x x x x x ┘                       . . . . x x x x
. . . . . . x x                         . . . . x x x x
  └───┘ └───┘                         
   G2    G1                           
        
In general one can move `b` bulges `k` steps further (in the 
example b = 2 and k = 1) and then accumulate all the 
reflections in a big matrix U.

02_example.jl

using LinearAlgebra
using LinearAlgebra: givensAlgorithm
using Base: OneTo

import LinearAlgebra: lmul!, rmul!

abstract type SmallRotation end

"""
Two Given's rotations acting on rows i:i+2. This could also 
be implemented as one Householder reflector, but let's assume 
we want good stability and we now how to get that just with
Given's rotations.
"""
struct Rotation3{Tc,Ts} <: SmallRotation
    c₁::Tc
    s₁::Ts
    c₂::Tc
    s₂::Ts
    i::Int
end

"""
    get_rotation(p₁, p₂, p₃, i) -> Rotation3, τ

Returns a double Givens rotation G such that G * [p₁; p₂; p₃] = [τ; 0; 0].
"""
function get_rotation(p₁, p₂, p₃, i::Int)
    c₁, s₁, nrm₁ = givensAlgorithm(p₂, p₃)
    c₂, s₂, nrm₂ = givensAlgorithm(p₁, nrm₁)
    Rotation3(c₁, s₁, c₂, s₂, i), nrm₂
end

@inline lmul!(G::SmallRotation, A::AbstractMatrix) = lmul!(G, A, axes(A, 2))
@inline rmul!(A::AbstractMatrix, G::SmallRotation) = rmul!(A, G, axes(A, 1))

@inline function lmul!(G::Rotation3, A::AbstractMatrix, range)
    @inbounds for j = range
        a₁ = A[G.i+0,j]
        a₂ = A[G.i+1,j]
        a₃ = A[G.i+2,j]

        a₂′ = G.c₁ * a₂ + G.s₁ * a₃
        a₃′ = -G.s₁' * a₂ + G.c₁ * a₃

        a₁′′ = G.c₂ * a₁ + G.s₂ * a₂′
        a₂′′ = -G.s₂' * a₁ + G.c₂ * a₂′
        
        A[G.i+0,j] = a₁′′
        A[G.i+1,j] = a₂′′
        A[G.i+2,j] = a₃′
    end

    A
end

@inline function rmul!(A::AbstractMatrix, G::Rotation3, range)
    @inbounds for j = range
        a₁ = A[j,G.i+0]
        a₂ = A[j,G.i+1]
        a₃ = A[j,G.i+2]

        a₂′ = a₂ * G.c₁ + a₃ * G.s₁'
        a₃′ = a₂ * -G.s₁ + a₃ * G.c₁

        a₁′′ = a₁ * G.c₂ + a₂′ * G.s₂'
        a₂′′ = a₁ * -G.s₂ + a₂′ * G.c₂

        A[j,G.i+0] = a₁′′
        A[j,G.i+1] = a₂′′
        A[j,G.i+2] = a₃′
    end
    A
end

function multiple_shifts(b::Int, k::Int)
    # k is the number of steps the shifts will move
    # b is the number of bulges
    # each bulge spans 3 columns
    # so we need pre-allocate a 3b + k matrix.
    n = 3b + k
    Q = Matrix(1.0I, n - 1, n - 1)

    # loop over the columns where each bulge start from front to back
    for start = 3b-2:-3:1
        # move the current bulge k steps forward 
        for j = start:start+k-1
            # generate a random rotation.
            G, = get_rotation(rand(), rand(), rand(), j)
            rmul!(Q, G)
        end
    end

    Q
end

function matrix_structure(A)
    nz = 0
    @inbounds for i = 1:size(A, 1)
        for j = 1:size(A, 2)
            if iszero(A[i, j])
                print(". ")
            else
                print("x ")
                nz += 1
            end
        end
        println()
    end
    println("Nonzero ratio: ", round.(nz / length(A), digits = 2))
end

03_example.jl

julia> matrix_structure(multiple_shifts(2, 6))
x x x x x x x . . . .
x x x x x x x x . . .
x x x x x x x x . . .
. x x x x x x x x x .
. x x x x x x x x x x
. x x x x x x x x x x
. . x x x x x x x x x
. . . x x x x x x x x
. . . . x x x x x x x
. . . . . x x x x x x
. . . . . . . . x x x
Nonzero ratio: 0.7

julia> matrix_structure(multiple_shifts(6, 10))
x x x x x x x x x x x . . . . . . . . . . . . . . . .
x x x x x x x x x x x x . . . . . . . . . . . . . . .
x x x x x x x x x x x x . . . . . . . . . . . . . . .
. x x x x x x x x x x x x x . . . . . . . . . . . . .
. x x x x x x x x x x x x x x . . . . . . . . . . . .
. x x x x x x x x x x x x x x . . . . . . . . . . . .
. . x x x x x x x x x x x x x x x . . . . . . . . . .
. . x x x x x x x x x x x x x x x x . . . . . . . . .
. . x x x x x x x x x x x x x x x x . . . . . . . . .
. . . x x x x x x x x x x x x x x x x x . . . . . . .
. . . x x x x x x x x x x x x x x x x x x . . . . . .
. . . x x x x x x x x x x x x x x x x x x . . . . . .
. . . . x x x x x x x x x x x x x x x x x x x . . . .
. . . . x x x x x x x x x x x x x x x x x x x x . . .
. . . . x x x x x x x x x x x x x x x x x x x x . . .
. . . . . x x x x x x x x x x x x x x x x x x x x x .
. . . . . x x x x x x x x x x x x x x x x x x x x x x
. . . . . x x x x x x x x x x x x x x x x x x x x x x
. . . . . . x x x x x x x x x x x x x x x x x x x x x
. . . . . . . x x x x x x x x x x x x x x x x x x x x
. . . . . . . . x x x x x x x x x x x x x x x x x x x
. . . . . . . . . x x x x x x x x x x x x x x x x x x
. . . . . . . . . . . . x x x x x x x x x x x x x x x
. . . . . . . . . . . . . . . x x x x x x x x x x x x
. . . . . . . . . . . . . . . . . . x x x x x x x x x
. . . . . . . . . . . . . . . . . . . . . x x x x x x
. . . . . . . . . . . . . . . . . . . . . . . . x x x
Nonzero ratio: 0.58

julia> matrix_structure(multiple_shifts(4, 20))
x x x x x x x x x x x x x x x x x x x x x . . . . . . . . . .
x x x x x x x x x x x x x x x x x x x x x x . . . . . . . . .
x x x x x x x x x x x x x x x x x x x x x x . . . . . . . . .
. x x x x x x x x x x x x x x x x x x x x x x x . . . . . . .
. x x x x x x x x x x x x x x x x x x x x x x x x . . . . . .
. x x x x x x x x x x x x x x x x x x x x x x x x . . . . . .
. . x x x x x x x x x x x x x x x x x x x x x x x x x . . . .
. . x x x x x x x x x x x x x x x x x x x x x x x x x x . . .
. . x x x x x x x x x x x x x x x x x x x x x x x x x x . . .
. . . x x x x x x x x x x x x x x x x x x x x x x x x x x x .
. . . x x x x x x x x x x x x x x x x x x x x x x x x x x x x
. . . x x x x x x x x x x x x x x x x x x x x x x x x x x x x
. . . . x x x x x x x x x x x x x x x x x x x x x x x x x x x
. . . . . x x x x x x x x x x x x x x x x x x x x x x x x x x
. . . . . . x x x x x x x x x x x x x x x x x x x x x x x x x
. . . . . . . x x x x x x x x x x x x x x x x x x x x x x x x
. . . . . . . . x x x x x x x x x x x x x x x x x x x x x x x
. . . . . . . . . x x x x x x x x x x x x x x x x x x x x x x
. . . . . . . . . . x x x x x x x x x x x x x x x x x x x x x
. . . . . . . . . . . x x x x x x x x x x x x x x x x x x x x
. . . . . . . . . . . . x x x x x x x x x x x x x x x x x x x
. . . . . . . . . . . . . x x x x x x x x x x x x x x x x x x
. . . . . . . . . . . . . . x x x x x x x x x x x x x x x x x
. . . . . . . . . . . . . . . x x x x x x x x x x x x x x x x
. . . . . . . . . . . . . . . . x x x x x x x x x x x x x x x
. . . . . . . . . . . . . . . . . x x x x x x x x x x x x x x
. . . . . . . . . . . . . . . . . . x x x x x x x x x x x x x
. . . . . . . . . . . . . . . . . . . x x x x x x x x x x x x
. . . . . . . . . . . . . . . . . . . . . . x x x x x x x x x
. . . . . . . . . . . . . . . . . . . . . . . . . x x x x x x
. . . . . . . . . . . . . . . . . . . . . . . . . . . . x x x
Nonzero ratio: 0.65