Fast qr for `BlockSkylineMatrix` and Cholesky for `Symmetric{<:Any,<:BlockSkylineMatrix}`

[dlfivefifty]

@jagot I'm going to tackle this soon which should speed up \ a lot. Essentially its not too hard to do this in-place, calling LAPACK to do the heavy work.

[jagot]

Since I'm at the point where I've transitioned FEDVRQuasi to Applied fast factorizations is really high on my wish-list, for solving Poisson's problem and eigenproblems.

[jagot]

One thing I forgot to mention, I need to update my factorizations every iteration, so it would be nice to be able to precompute the structure of the factorization (since that will not change), and every iteration do something along the lines of

copyto!(factorization, Factorization(H)

where H is a BlockSkylineMatrix.

An example that works with latest FEDVRQuasi, which I'd like to be fast:

using FEDVRQuasi
using LinearAlgebra
using ArnoldiMethod

L = 40.0
ℓ = 1 # Angular momentum
Z = 1 # Nuclear charge

t = range(0.0,stop=L,length=10)
R = FEDVR(t, 10)[:,2:end-1]

D = Derivative(axes(R,1))

∇² = R' * D' * D * R

T = -0.5∇²
V = Matrix(r -> -Z/r + ℓ*(ℓ+1)/2r^2, R)

H = T + V
# Aim slightly below the ground state eigenvalue
σ = 1.1*(-Z^2/2)
H⁻¹ = factorize(Symmetric(Matrix(H - σ*I)))

Diagonalizing:

struct ShiftInvert{Factorization, T}
    A⁻¹::Factorization
    σ::T
end

Base.size(S::ShiftInvert, args...) = size(S.A⁻¹, args...)
Base.eltype(S::ShiftInvert) = eltype(S.A⁻¹)

LinearAlgebra.mul!(y, si::ShiftInvert, x) = ldiv!(y, si.A⁻¹, x)

println("Direct:")
schurQR,history = @time partialschur(H,nev=1,which=SR())
println(history)

println("Shift-and-invert:")
schurQR,history = @time partialschur(ShiftInvert(H⁻¹,σ),nev=1,which=LR())
println(history)

gives

Direct:
  0.004918 seconds (61.18 k allocations: 3.942 MiB)
Converged: 1 of 1 eigenvalues in 230 matrix-vector products
Shift-and-invert:
  0.157293 seconds (262.31 k allocations: 12.258 MiB)
Converged: 1 of 1 eigenvalues in 20 matrix-vector products

[dlfivefifty]

Perhaps I'm not understanding but factorizations are typically in-place. So I think the code would actually be something like:

F_data .= H - σ*I # or applied(-, H, σ*I) to avoid allocation 
F = cholesky!(Symmetric(F_data))

PS Is Cholesky more useful for you? That is, do you expect to mostly work with positive semi-definite? It should actually be slightly easier to implement than QR and much faster.

[jagot]

For hydrogen (or any one-electron problem), you can always choose a shift such that the Hermitian matrix becomes positive-definite, it is more complicated for multi-electron problems, where you need to optimize multiple (coupled) eigenvectors, belonging to different eigenvalues at the same time. Then, for some of the interior eigenvalues, the shifted matrix is not positive definite any more. Furthermore, the integral part of the integro-differential operator is not easily factorizable, so for problems where integral operators arise, my strategy is to use IterativeSolvers.jl with a preconditioner built from the part that I can factorize, i.e. the one-body operator above and some extra potentials.

At the moment I'm doing straight Arnoldi, but that has convergence issues and high condition numbers. Simple shift-and-invert does not work very well either, since the spectrum is very dense. I'm thinking maybe I need inverse iterations, which also needs factorizations as above.

EDIT: TL-DR, Cholesky is nice, but I need the more general case as well.

[MikaelSlevinsky]

I think you want LDL^T, which should be even faster than Cholesky for symmetric block-skyline matrices.

[jagot]

Yes, that's what I get back for finite-differences, where the Hamiltonian is symmetric tridiagonal.