import Pkg; Pkg.add(["GenericLinearAlgebra", "SparseArrays", "Profile", "ProfileView", "BenchmarkTools", "PProf"])
using LinearAlgebra
using SparseArrays

# Random real matrix -> will do an LU
A = randn(10, 10)
@show typeof(factorize(A))

# Real-symmetric matrix ->  will do a Bunch-Kaufman
Am = Symmetric(A + A')
@show typeof(factorize(Am))

# Symmetric tridiagonal -> will do a LDLt
Am = SymTridiagonal(A + A')
@show typeof(factorize(Am))

# Random sparse matrix -> will do sparse LU
S = sprandn(50, 50, 0.3)
@show typeof(factorize(S))

# ... and so on ...

function solve_many(A, xs)
    F = factorize(A)
    [F \ rhs for rhs in xs]
end

using LinearAlgebra

function davidson(A, SS::AbstractArray; maxiter=100, prec=I,
                  tol=20size(A,2)*eps(eltype(A)),
                  maxsubspace=8size(SS, 2))
    m = size(SS, 2)
    for i in 1:100
        Ass = A * SS

        # Use eigen specialised for Hermitian matrices
        rvals, rvecs = eigen(Hermitian(SS' * Ass))
        rvals = rvals[1:m]
        rvecs = rvecs[:, 1:m]
        Ax = Ass * rvecs

        R = Ax - SS * rvecs * Diagonal(rvals)
        if norm(R) < tol
            return rvals, SS * rvecs
        end

        println(i, "  ", size(SS, 2), "  ", norm(R))

        # Use QR to orthogonalise the subspace.
        if size(SS, 2) + m > maxsubspace
            SS = typeof(R)(qr(hcat(SS * rvecs, prec * R)).Q)
        else
            SS = typeof(R)(qr(hcat(SS, prec * R)).Q)
        end
    end
    error("not converged.")
end

nev = 2
A = randn(20, 20); A = A + A' + I;

# Generate two random orthogonal guess vectors
x0 = randn(size(A, 2), nev)
x0 = Array(qr(x0).Q)

# Run the problem
davidson(A, x0)

using GenericLinearAlgebra  # (Slow) generic fallbacks for eigen and qr

λ, v = davidson(Matrix{Float16}(A), Float16.(x0), tol=1)
println()
λ, v = davidson(Matrix{Float32}(A), v, tol=1e-3)
println()
λ, v = davidson(Matrix{Float64}(A), v, tol=1e-13)
println()
λ, v = davidson(Matrix{BigFloat}(A), v, tol=1e-25)
λ

nev = 2
spA = sprandn(100, 100, 0.3); spA = spA + spA' + 2I
spx0 = randn(size(spA, 2), nev)
spx0 = Array(qr(spx0).Q)

davidson(spA, spx0, tol=1e-6)

# Use the GPU
using CuArrays
using CuArrays.CUSOLVER

function LinearAlgebra.eigen(A::Hermitian{T, CuArray{T, 2}}) where {T <: Real}
    CUSOLVER.syevd!('V', 'U', A.data)
end

davidson(cu(A), 2)

using Profile
using ProfileView
using PProf

# Setup of the problem:
nev = 2
A = randn(20, 20); A = A + A' + I;
x0 = randn(size(A, 2), nev)
x0 = Array(qr(x0).Q)

# Run once to compile everything ... this should be ignored
@profview davidson(A, x0)
@profview davidson(A, x0)

# Run once to compile everything ... this should be ignored
@pprof davidson(A, x0)
@pprof davidson(A, x0)

using BenchmarkTools

SS = randn(20, 16)
precR = randn(20, 2)

catted = @btime hcat(SS, precR)
fac = @btime qr(catted)
@btime Array(fac.Q);