14544069614

Committed 19 Apr 2025 12:57AM UTC coverage: 94.777% (-0.009%) from 94.786%

Build # 14544069614

Build Type

push

github

Committed by

web-flow

Commit Message

Better exact diagonalisation (#306)

# Changes
* Implement a `LinearMap` from
[LinearMaps.jl](https://github.com/JuliaLinearAlgebra/LinearMaps.jl),
which allows us to use operators as matrices without storing the
matrices in memory.
* Use the `LinearMap` to implement matrix-free variants of
`KrylovKitSolver`, `ArpackSolver`, and `LOBPCGSolver`.
* Refactor of the `ExactDiagonalizationProblem` and its surrounding
infrastructure.

# Changed behaviour
* `ArpackSolver` and `LOBPCGSolver` now use the matrix-free algorithm by
default.

Run Details

208 of 209 new or added lines in 9 files covered. (99.52%)

19 existing lines in 3 files now uncovered.

6932 of 7314 relevant lines covered (94.78%)

11143439.58 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

83.72

/src/ExactDiagonalization/solve.jl

# a lazy type for iterating over coefficient vectors starting from a vector of DVecs
struct LazyCoefficientVectorsDVecs{T,VDV<:Vector{<:AbstractDVec},B} <: AbstractVector{T}
    vecs::VDV
    basis::B
end
function LazyCoefficientVectorsDVecs(vecs, basis)
    T = valtype(vecs[1])
    return LazyCoefficientVectorsDVecs{T,typeof(vecs),typeof(basis)}(vecs, basis)
end
Base.size(v::LazyCoefficientVectorsDVecs) = (length(v.vecs),)
function Base.getindex(lcv::LazyCoefficientVectorsDVecs{T}, i::Int) where {T}
    return T[lcv.vecs[i][address] for address in lcv.basis]
end

# lazy type for iterating over constructed DVecs
struct LazyDVecs{DV,C<:AbstractVector{<:AbstractVector},B} <: AbstractVector{DV}
    coefficient_vecs::C
    basis::B
end
function LazyDVecs(vs::LCV, basis::B) where {LCV,B}
    K = eltype(basis)
    V = eltype(eltype(vs))
    return LazyDVecs{typeof(PDVec{K,V}()),LCV,B}(vs, basis)
end
Base.size(lv::LazyDVecs) = size(lv.coefficient_vecs)
Base.getindex(lv::LazyDVecs, i::Int) = PDVec(zip(lv.basis, lv.coefficient_vecs[i]))

# a generic result type for ExactDiagonalizationProblem
struct EDResult{A,P,VA<:AbstractVector,VE<:AbstractVector,CV,B,I,R}
    algorithm::A
    problem::P
    values::VA
    vectors::VE
    coefficient_vectors::CV
    basis::B
    info::I
    howmany::Int
    raw::R # algorithm-specific raw result, e.g. the matrix of eigenvectors
    success::Bool
end
# iteration for destructuring into components
Base.iterate(S::EDResult) = (S.values, Val(:vectors))
Base.iterate(S::EDResult, ::Val{:vectors}) = (S.vectors, Val(:success))
Base.iterate(S::EDResult, ::Val{:success}) = (S.info, Val(:done))
Base.iterate(::EDResult, ::Val{:done}) = nothing

function Base.show(io::IO, r::EDResult)
    io = IOContext(io, :compact => true)
    n = length(r.values)
    println(io, "EDResult for algorithm $(r.algorithm) with $n eigenvalue(s),")
    print(io, "  values = ")
    show(io, r.values)
    print(io, "\n  Convergence info: ")
    show(io, r.info)
    print(io, ", with howmany = $(r.howmany) eigenvalues requested.")
    print(io, "\n  success = $(r.success).")
end

# solve directly on the ExactDiagonalizationProblem
"""
    solve(p::ExactDiagonalizationProblem, [algorithm]; kwargs...)

Solve an [`ExactDiagonalizationProblem`](@ref) `p` directly. Optionally specify an
`algorithm.` Returns a result type with the eigenvalues, eigenvectors, and convergence
information.

For a description of the keyword arguments, see the documentation for
[`ExactDiagonalizationProblem`](@ref).

See also
[`solve(::ProjectorMonteCarloProblem)`](@ref CommonSolve.solve(::ProjectorMonteCarloProblem)).
"""
function CommonSolve.solve(p::ExactDiagonalizationProblem; kwargs...)
    s = init(p; kwargs...)
    return solve(s)
end
function CommonSolve.solve(p::ExactDiagonalizationProblem, algorithm; kwargs...)
    s = init(p, algorithm; kwargs...)
    return solve(s)
end

# The code for `CommonSolve.solve(::IterativeEDSolver{<:ALG}; ...)` for
# - ALG<:KrylovKitSolver is part of the `KrylovKitExt.jl` extension
# - ALG<:ArpackSolver is part of the `ArpackExt.jl` extension
# - ALG<:LOBPCGSolver is part of the `IterativeSolversExt.jl` extension


function CommonSolve.solve(s::DenseEDSolver; kwargs...)
    # Combine and clean keyword arguments
    kwargs = (; verbose=false, s.solver_kwargs..., kwargs...)
    eigen_kwargs, rest = split_keys(kwargs, :permute, :scale, :sortby)
    (; verbose) = clean_and_warn_if_others_present(rest, (:verbose,))

    eigen_factorization = eigen!(Matrix(s.basis_set_rep.sparse_matrix); eigen_kwargs...)

    verbose && @info "LinearAlgebra.eigen: success"

    coefficient_vectors = eachcol(eigen_factorization.vectors)
    vectors = LazyDVecs(coefficient_vectors, s.basis_set_rep.basis)
    return EDResult(
        s.algorithm,
        s.problem,
        eigen_factorization.values,
        vectors,
        coefficient_vectors,
        s.basis_set_rep.basis,
        "Dense matrix eigensolver solution from `LinearAlgebra.eigen`",
        dimension(s.basis_set_rep),
        eigen_factorization.vectors,
        true # successful if no exception was thrown
    )
end

1	# a lazy type for iterating over coefficient vectors starting from a vector of DVecs
2	struct LazyCoefficientVectorsDVecs{T,VDV<:Vector{<:AbstractDVec},B} <: AbstractVector{T}
3	vecs::VDV
4	basis::B
5	end
UNCOV 6	function LazyCoefficientVectorsDVecs(vecs, basis)	×
UNCOV 7	T = valtype(vecs[1])	×
UNCOV 8	return LazyCoefficientVectorsDVecs{T,typeof(vecs),typeof(basis)}(vecs, basis)	×
9	end
UNCOV 10	Base.size(v::LazyCoefficientVectorsDVecs) = (length(v.vecs),)	×
UNCOV 11	function Base.getindex(lcv::LazyCoefficientVectorsDVecs{T}, i::Int) where {T}	×
UNCOV 12	return T[lcv.vecs[i][address] for address in lcv.basis]	×
13	end
14
15	# lazy type for iterating over constructed DVecs
16	struct LazyDVecs{DV,C<:AbstractVector{<:AbstractVector},B} <: AbstractVector{DV}
17	coefficient_vecs::C	253✔
18	basis::B
19	end
20	function LazyDVecs(vs::LCV, basis::B) where {LCV,B}	253✔
21	K = eltype(basis)	253✔
22	V = eltype(eltype(vs))	253✔
23	return LazyDVecs{typeof(PDVec{K,V}()),LCV,B}(vs, basis)	253✔
24	end
25	Base.size(lv::LazyDVecs) = size(lv.coefficient_vecs)	72✔
26	Base.getindex(lv::LazyDVecs, i::Int) = PDVec(zip(lv.basis, lv.coefficient_vecs[i]))	544✔
27
28	# a generic result type for ExactDiagonalizationProblem
29	struct EDResult{A,P,VA<:AbstractVector,VE<:AbstractVector,CV,B,I,R}
30	algorithm::A	253✔
31	problem::P
32	values::VA
33	vectors::VE
34	coefficient_vectors::CV
35	basis::B
36	info::I
37	howmany::Int
38	raw::R # algorithm-specific raw result, e.g. the matrix of eigenvectors
39	success::Bool
40	end
41	# iteration for destructuring into components
42	Base.iterate(S::EDResult) = (S.values, Val(:vectors))	1✔
43	Base.iterate(S::EDResult, ::Val{:vectors}) = (S.vectors, Val(:success))	1✔
44	Base.iterate(S::EDResult, ::Val{:success}) = (S.info, Val(:done))	1✔
45	Base.iterate(::EDResult, ::Val{:done}) = nothing	×
46
47	function Base.show(io::IO, r::EDResult)	25✔
48	io = IOContext(io, :compact => true)	25✔
49	n = length(r.values)	25✔
50	println(io, "EDResult for algorithm $(r.algorithm) with $n eigenvalue(s),")	25✔
51	print(io, " values = ")	25✔
52	show(io, r.values)	25✔
53	print(io, "\n Convergence info: ")	25✔
54	show(io, r.info)	25✔
55	print(io, ", with howmany = $(r.howmany) eigenvalues requested.")	25✔
56	print(io, "\n success = $(r.success).")	25✔
57	end
58
59	# solve directly on the ExactDiagonalizationProblem
60	"""
61	solve(p::ExactDiagonalizationProblem, [algorithm]; kwargs...)
62
63	Solve an [`ExactDiagonalizationProblem`](@ref) `p` directly. Optionally specify an
64	`algorithm.` Returns a result type with the eigenvalues, eigenvectors, and convergence
65	information.
66
67	For a description of the keyword arguments, see the documentation for
68	[`ExactDiagonalizationProblem`](@ref).
69
70	See also
71	[`solve(::ProjectorMonteCarloProblem)`](@ref CommonSolve.solve(::ProjectorMonteCarloProblem)).
72	"""
73	function CommonSolve.solve(p::ExactDiagonalizationProblem; kwargs...)	58✔
74	s = init(p; kwargs...)	29✔
75	return solve(s)	29✔
76	end
77	function CommonSolve.solve(p::ExactDiagonalizationProblem, algorithm; kwargs...)	422✔
78	s = init(p, algorithm; kwargs...)	211✔
79	return solve(s)	211✔
80	end
81
82	# The code for `CommonSolve.solve(::IterativeEDSolver{<:ALG}; ...)` for
83	# - ALG<:KrylovKitSolver is part of the `KrylovKitExt.jl` extension
84	# - ALG<:ArpackSolver is part of the `ArpackExt.jl` extension
85	# - ALG<:LOBPCGSolver is part of the `IterativeSolversExt.jl` extension
86
87
88	function CommonSolve.solve(s::DenseEDSolver; kwargs...)	124✔
89	# Combine and clean keyword arguments
90	kwargs = (; verbose=false, s.solver_kwargs..., kwargs...)	124✔
91	eigen_kwargs, rest = split_keys(kwargs, :permute, :scale, :sortby)	62✔
92	(; verbose) = clean_and_warn_if_others_present(rest, (:verbose,))	62✔
93
94	eigen_factorization = eigen!(Matrix(s.basis_set_rep.sparse_matrix); eigen_kwargs...)	94✔
95
96	verbose && @info "LinearAlgebra.eigen: success"	62✔
97
98	coefficient_vectors = eachcol(eigen_factorization.vectors)	98✔
99	vectors = LazyDVecs(coefficient_vectors, s.basis_set_rep.basis)	83✔
100	return EDResult(	62✔
101	s.algorithm,
102	s.problem,
103	eigen_factorization.values,
104	vectors,
105	coefficient_vectors,
106	s.basis_set_rep.basis,
107	"Dense matrix eigensolver solution from `LinearAlgebra.eigen`",
108	dimension(s.basis_set_rep),
109	eigen_factorization.vectors,
110	true # successful if no exception was thrown
111	)
112	end

RimuQMC / Rimu.jl / 14544069614

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous