From ded786d1d61202fcd7f4cab85169b349e17b33f6 Mon Sep 17 00:00:00 2001 From: Nick Mayhall Date: Sat, 10 Jun 2023 09:26:35 -0400 Subject: [PATCH] preconditioning working ok it seems, still tough to get very high accuracy --- src/BlockDavidson.jl | 105 ++++------------------------------------ test/runtests.jl | 1 + test/test_precond_01.jl | 33 +++++++++++-- 3 files changed, 38 insertions(+), 101 deletions(-) diff --git a/src/BlockDavidson.jl b/src/BlockDavidson.jl index 757cfe7..b535976 100644 --- a/src/BlockDavidson.jl +++ b/src/BlockDavidson.jl @@ -135,89 +135,6 @@ function apply_diagonal_precond!(res_s::Vector{T}, Adiag::Vector{T}, denom::T) w end function iteration(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thresh=1.0) - - # - # perform sig_curr = A*vec_curr - solver.sig_curr = Matrix(solver.op * solver.vec_curr) - - # - # add these new vectors to previous quantities - solver.sig_prev = hcat(solver.sig_prev, solver.sig_curr) - solver.vec_prev = hcat(solver.vec_prev, solver.vec_curr) - - # - # Check orthogonality - ss_metric = solver.vec_prev'*solver.vec_prev - solver.lindep = abs(1.0 - det(ss_metric)) - - # - # form H in current subspace - Hss = solver.vec_prev' * solver.sig_prev - F = eigen(Hss) - #ritz_e = F.values[1:solver.nroots] - #ritz_v = F.vectors[:,1:solver.nroots] - ritz_e = F.values - ritz_v = F.vectors - ss_size = length(F.values) - if ss_size >= solver.max_ss_vecs - ritz_v = ritz_v[:,sortperm(ritz_e)][:,1:solver.nroots] - ritz_e = ritz_e[sortperm(ritz_e)][1:1:solver.nroots] - #ritz_v = ritz_v[:,sortperm(ritz_e)][:,1:solver.max_ss_vecs] - #ritz_e = ritz_e[sortperm(ritz_e)][1:1:solver.max_ss_vecs] - else - ritz_v = ritz_v[:,sortperm(ritz_e)] - ritz_e = ritz_e[sortperm(ritz_e)] - end - - solver.ritz_e = ritz_e - solver.ritz_v = ritz_v - - solver.sig_prev = solver.sig_prev * ritz_v - solver.vec_prev = solver.vec_prev * ritz_v - Hss = ritz_v' * Hss * ritz_v - - res = deepcopy(solver.sig_prev[:,1:solver.nroots]) - for s in 1:solver.nroots - res[:,s] .-= solver.vec_prev[:,s] * Hss[s,s] - end - - - ss = deepcopy(solver.vec_prev) - - new = zeros(solver.dim, 0) - #solver.statusconv = [false for i in 1:solver.nroots] - for s in 1:solver.nroots - # |v'> = (I-|SS> - # = |v> - |SS> - solver.resid[s] = norm(res[:,s]) - if norm(res[:,s]) <= solver.tol - solver.status[s] = true - continue - else - solver.status[s] = false - end - if Adiag != nothing - #level_shift = 1e-3 - #denom = (ritz_e[s] + level_shift) - tmp = deepcopy(res) - # apply_diagonal_precond!(res[:,s], Adiag, ritz_e[s]) - # for i in 1:length(Adiag) - # res[i,s] = res[i,s] / (ritz_e[s] - Adiag[i]) - # end - end - scr = ss' * res[:,s] - res[:,s] .= res[:,s] .- (ss * scr) - nres = norm(res[:,s]) - if nres > solver.tol - ss = hcat(ss,res[:,s]/nres) - new = hcat(new,res[:,s]/nres) - end - end - return new -end - -#function iteration2(solver::Davidson, v_prev, v_curr, σ_prev; Adiag=nothing, iprint=0) -function iteration2(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thresh=1.0) # # project out v_prev from v_curr @@ -245,9 +162,6 @@ function iteration2(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thr idx = sortperm(F.values) ss_size = min(size(solver.vec_prev,2), solver.max_ss_vecs) - if ss_size >= solver.max_ss_vecs - ss_size = solver.nroots - end idx = idx[1:ss_size] @@ -290,19 +204,25 @@ function iteration2(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thr end end + res = res - solver.vec_prev * (solver.vec_prev' * res) return Matrix(qr(res).Q) end -function eigs(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thresh=1.0) +""" + eigs(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thresh=1e-1) + +TBW +""" +function eigs(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thresh=1e-1) for iter = 1:solver.max_iter #@btime $solver.vec_curr = $iteration($solver) - solver.vec_curr = iteration2(solver, Adiag=Adiag, iprint=iprint, precond_start_thresh=precond_start_thresh) + solver.vec_curr = iteration(solver, Adiag=Adiag, iprint=iprint, precond_start_thresh=precond_start_thresh) solver.iter_curr = iter print_iter(solver) if all(solver.status) - solver.converged == true + solver.converged = true break end if solver.lindep > solver.lindep_thresh && iter < solver.max_iter @@ -322,11 +242,4 @@ function eigs(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thresh=1. #return solver.fritz_e[1:solver.nroots], solver.vec_prev*solver.ritz_v[:,1:solver.nroots] end -#function orthogonalize!(solver) -# -# # |vv_i> = |v_i> - |w -# new_v = zeros(solver.dim, nroots) -# for r in 1:nroots -#end - end diff --git a/test/runtests.jl b/test/runtests.jl index e729f2f..8e6ade1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,4 +3,5 @@ using Test @testset "BlockDavidson.jl" begin include("test01.jl") + include("test_precond_01.jl") end diff --git a/test/test_precond_01.jl b/test/test_precond_01.jl index ffc06ca..d528d3c 100644 --- a/test/test_precond_01.jl +++ b/test/test_precond_01.jl @@ -1,11 +1,34 @@ using BlockDavidson using LinearAlgebra using Random +using Test -Random.seed!(2) +@testset "BlockDavidson preconditioning" begin + Random.seed!(2) -N = 10_000 -A = -10 .* diagm(rand(N)) + .001 * rand(N,N) -A += A' + N = 1_000 + nr = 4 + A = -10 .* diagm(rand(N)) + 0.001 * rand(N, N) + A += A' -v0 = qr(rand(N,2)).Q[:,1:3]; \ No newline at end of file + v0 = qr(rand(N, 2)).Q[:, 1:nr] + + idx = sortperm(diag(A)) + v0 = Matrix(1.0I, N, N)[:, idx[1:nr]] + + + println(" No preconditioning") + dav1 = Davidson(A, max_iter=50, nroots=nr, tol=1e-8, v0=v0) + dav2 = Davidson(A, max_iter=50, nroots=nr, tol=1e-8, v0=v0) + e1, v1 = eigs(dav1) + println() + println(" Preconditioned") + e2, v2 = eigs(dav2, Adiag=diag(A)) + + flush(stdout) + display(dav1.converged) + display(dav2.converged) + @test dav1.converged == false + @test dav2.converged == true + +end \ No newline at end of file