Cleanup

Author: mfherbst

src/response/chi0.jl

@@ -108,7 +108,7 @@ precondprep!(P::FunctionPreconditioner, ::Any) = P
 # /!\ It is assumed (and not checked) that ψk'Hk*ψk = Diagonal(εk) (extra states
 # included).
 function sternheimer_solver(Hk, ψk, ε, rhs;
-                            callback=identity, cg_callback=identity,
+                            callback=identity,
                             ψk_extra=zeros_like(ψk, size(ψk, 1), 0), εk_extra=zeros(0),
                             Hψk_extra=zeros_like(ψk, size(ψk, 1), 0), tol=1e-9,
                             miniter=1, maxiter=100)
@@ -179,10 +179,9 @@ function sternheimer_solver(Hk, ψk, ε, rhs;
     end
     J = LinearMap{eltype(ψk)}(RAR, size(Hk, 1))
     cg_res = cg(J, bb; precon=FunctionPreconditioner(R_ldiv!), tol, proj=R,
-                callback=cg_callback, miniter, maxiter)
+                callback=info -> callback(merge(info, (; basis, kpoint))),
+                miniter, maxiter)
     δψknᴿ = cg_res.x
-    info = (; basis, kpoint, cg_res)
-    callback(info)
 
     # 2) solve for αkn now that we know δψknᴿ
     # Note that αkn is an empty array if there is no extra bands.
@@ -618,10 +617,10 @@ end
 Set the convergence tolerance for each Sternheimer solver to be a fixed factor
 times the desired density tolerance.
 """
-struct SternheimerFactor{T}
+struct BandtolFactor{T}
     factor::T
 end
-SternheimerFactor() = SternheimerFactor(1)
-function determine_band_tolerances(alg::SternheimerFactor, tol_density)
+BandtolFactor() = BandtolFactor(1)
+function determine_band_tolerances(alg::BandtolFactor, tol_density)
     alg.factor * tol_density
 end

src/response/gmres.jl

@@ -77,7 +77,6 @@ function inexact_gmres!(x, A, b::AbstractVector{T};
                  residual_norm=β, maxiter, tol, s, residual=r,
                  restart_history, stage=:restart, krylovdim, k, Axinfo)
 
-
         while (n_iter < maxiter && k < m)  # Arnoldi loop
             n_iter += 1
 
@@ -145,263 +144,3 @@ end
 function inexact_gmres(A, b; kwargs...)
     inexact_gmres!(zeros_like(b, eltype(b), size(A, 2)), A, b; kwargs...)
 end
-
-
-
-
-
-
-
-
-
-
-function gmres(operators::Function, b::AbstractVector{T}; 
-               x₀=nothing,
-               maxiter=min(100, size(operators(0), 1)),
-               restart=min(20, size(operators(0), 1)),
-               tol=1e-6, s_guess=0.8, verbose=1, debug=false) where {T}
-    # pass s/3m 1/\|r_tilde\| tol to operators()
-    normb = norm(b)
-    b = b ./ normb
-    tol = tol / normb
-    if tol < 0.5*eps(T)
-        #throw a warning
-        @warn "The effective tolerance is smaller than half of the machine epsilon, the requested accuracy may not be achieved."
-    end
-
-    n_size = length(b)
-
-    # counters
-    numrestarts::Int64 = 0
-    restart_inds = Int64[]
-    MvTime = Float64[]
-
-    # initialization: Arnoldi basis and Hessenberg matrix
-    V = zeros(n_size, restart + 1)
-    H = zeros(restart + 1, restart)
-
-    # residual history
-    residuals = zeros(maxiter)
-    denominatorvec = zeros(restart)
-    denominator2 = 1.0
-    # record min singular values
-    minsvdvals = zeros(maxiter)
-
-    # solution history
-    if debug
-        X = zeros(n_size, maxiter)
-    end
-
-    # if x0 is zero, set r0 = b
-    x0iszero = isnothing(x₀)
-    if x0iszero
-        β₀ = 1.0
-        V[:, 1] = b #/ β₀
-        tol_new = tol / 2
-        lm = s_guess / restart / 2 * tol
-    else
-        # assert x₀ is a vector of the right size
-        @assert length(x₀) == n_size
-        x₀ ./= normb
-        x0iszero = false
-        tol_new = tol / 3
-        A0_tol = tol / 3
-        lm = s_guess / restart / 3 * tol
-        push!(MvTime, @elapsed r₀ = operators(A0_tol) * x₀)
-        r₀ = b - r₀
-        β₀ = norm(r₀)
-        V[:, 1] = r₀ / β₀
-    end
-
-    Ai_tol = lm / β₀ 
-    i = 0 # iteration counter, set to 0 when restarted
-    for j = 1:maxiter
-        i += 1
-
-        # Arnoldi iteration using MGS
-        push!(MvTime, @elapsed w = operators(Ai_tol) * V[:, i])
-        for j = 1:i
-            H[j, i] = dot(V[:, j], w)
-            w .-= H[j, i] * V[:, j]
-        end
-        H[i+1, i] = norm(w)
-        V[:, i+1] = w / H[i+1, i]
-
-        # compute the min singular value of the upper Hessenberg matrix
-        minsvdvals[j] = svdvals(H[1:i+1, 1:i])[end]
-    
-        # update residual history
-        if i == 1
-            denominatorvec[1] = conj(-H[1, 1] / H[2, 1])
-        else
-            denominatorvec[i] = conj(-(H[1, i] + dot(denominatorvec[1:i-1], H[2:i, i])) / H[i+1, i])
-        end
-        denominator2 += abs2(denominatorvec[i])
-        residuals[j] = β₀ / sqrt(denominator2)
-        Ai_tol = lm / residuals[j]
-
-        (verbose > 0) && println("restart = ", numrestarts + 1, ", iteration = ", i, ", res = ", residuals[j]*normb, "\n")
-
-        # when converged without updating the s_guess, restart once
-        converged = residuals[j] < tol_new
-        happy = converged && ((numrestarts > 0) || (s_guess <= minsvdvals[j]))
-        needrestart = (i == restart) || (converged && !happy) # equal to converged && (numrestarts == 0) && (s_guess > minsvdvals[j])
-        reachmaxiter = (j == maxiter)
-        exitloop = (happy || reachmaxiter)
-
-        solvels = (exitloop || needrestart || debug)
-        if solvels
-            # solve the Hessenberg least squares problem
-            rhs = (UpperTriangular(H[2:i+1, 1:i]) \ denominatorvec[1:i]) * (-β₀ / denominator2)
-            if x0iszero
-                xᵢ = V[:, 1:i] * rhs
-            else
-                xᵢ = x₀ + V[:, 1:i] * rhs
-            end
-
-            debug && (X[:, j] = xᵢ)
-            
-            if exitloop
-                if debug
-                    return (; x=X[:, 1:j].*normb, residuals=residuals[1:j].*normb, MvTime=MvTime, restart_inds=restart_inds, minsvdvals=minsvdvals[1:j])
-                else
-                    return (; x=xᵢ.*normb, residuals=residuals[1:j].*normb, MvTime=MvTime, restart_inds=restart_inds)
-                end
-            end
-
-            if needrestart
-                i = 0
-                x0iszero = false
-                # whenever restart, update the guess of the smallest singular value
-                min_s = minimum(minsvdvals[1:j])
-                if s_guess > min_s
-                    verbose > 0 && println("current s_guess = ", s_guess, "> minsvdval = ", min_s, ", restart with s_guess = minsvdval")
-                else
-                    verbose > 0 && println("current s_guess = ", s_guess, "≤ minsvdval = ", min_s, ", restart with s_guess = minsvdval")
-                end
-                s_guess = min_s
-                lm = s_guess / restart / 3 * tol
-                tol_new = tol / 3
-                A0_tol = tol / 3
-                
-                x₀ = xᵢ
-                push!(MvTime, @elapsed r₀ = operators(A0_tol) * x₀)
-                r₀ = b - r₀
-                β₀ = norm(r₀)
-                V[:, 1] = r₀ / β₀
-                numrestarts += 1
-                denominatorvec = zeros(maxiter)
-                denominator2 = 1.0
-                push!(restart_inds, j)
-                (verbose > -1) && println("restarting: ", "iteration = ", j, ", r₀ = ", β₀, "\n")
-                Ai_tol = lm / β₀ 
-            end
-        end
-    end
-end
-
-
-function gmres(operator::LinearMap{T}, b::AbstractVector{T}; 
-               x₀ = nothing,
-               maxiter=min(100, size(operator, 1)),
-               restart=min(20, size(operator, 1)),
-               tol=1e-6, verbose=0, debug=false) where {T}
-
-    n_size = length(b)
-
-    # counters
-    numrestarts::Int64 = 0
-    restart_inds = Int64[]
-    MvTime = Float64[]
-
-    # initialization: Arnoldi basis and Hessenberg matrix
-    V = zeros(n_size, restart + 1)
-    H = zeros(restart + 1, restart)
-
-    # residual history
-    residuals = zeros(maxiter)
-    denominatorvec = zeros(restart)
-    denominator2 = 1.0
-
-    # solution history
-    if debug
-        X = zeros(n_size, maxiter)
-    end
-
-    # if x0 is zero, set r0 = b
-    if isnothing(x₀)
-        β₀ = norm(b)
-        V[:, 1] = b / β₀
-    else
-        # assert x₀ is a vector of the right size
-        @assert length(x₀) == n_size
-        push!(MvTime, @elapsed r₀ = operator * x₀)
-        r₀ = b - r₀
-        β₀ = norm(r₀)
-        V[:, 1] = r₀ / β₀
-    end
-
-    i = 0 # iteration counter, set to 0 when restarted
-    for j = 1:maxiter
-        i += 1
-        # Arnoldi iteration using MGS
-        push!(MvTime, @elapsed w = operator * V[:, i])
-        for j = 1:i
-            H[j, i] = dot(V[:, j], w)
-            w .-= H[j, i] * V[:, j]
-        end
-        H[i+1, i] = norm(w)
-        V[:, i+1] = w / H[i+1, i]
-
-        # update residual history
-        if i == 1
-            denominatorvec[1] = conj(-H[1, 1] / H[2, 1])
-        else
-            denominatorvec[i] = conj(-(H[1, i] + dot(denominatorvec[1:i-1], H[2:i, i])) / H[i+1, i])
-        end
-        denominator2 += abs2(denominatorvec[i])
-        residuals[j] = β₀ / sqrt(denominator2)
-        
-        (verbose > 0) && println("restart = ", numrestarts + 1, ", iteration = ", i, ", res = ", residuals[j], "\n")
-
-        happy = residuals[j] < tol
-        needrestart = (i == restart)
-        reachmaxiter = (j == maxiter)
-        exitloop = (happy || reachmaxiter)
-
-        solvels = (exitloop || needrestart || debug)
-        if solvels
-            # solve the Hessenberg least squares problem
-            rhs = (UpperTriangular(H[2:i+1, 1:i]) \ denominatorvec[1:i]) * (-β₀ / denominator2)
-            if isnothing(x₀)
-                xᵢ = V[:, 1:i] * rhs
-            else
-                xᵢ = x₀ + V[:, 1:i] * rhs
-            end
-
-            debug && (X[:, j] = xᵢ)
-            
-            if exitloop
-                if debug
-                    return (; x=X[:, 1:j], residuals=residuals[1:j], MvTime=MvTime, restart_inds=restart_inds)
-                else
-                    return (; x=xᵢ, residuals=residuals[1:j], MvTime=MvTime, restart_inds=restart_inds)
-                end
-            end
-
-            if needrestart
-                i = 0
-                x₀ = xᵢ
-                push!(MvTime, @elapsed r₀ = operator * x₀)
-                r₀ = b - r₀
-                β₀ = norm(r₀)
-                V[:, 1] = r₀ / β₀
-                numrestarts += 1
-                denominatorvec = zeros(maxiter)
-                denominator2 = 1.0
-                push!(restart_inds, j)
-                (verbose > 0) && println("restarting: ", "iteration = ", j, ", r₀ = ", β₀, "\n")
-            end
-        end
-    end
-end

src/response/hessian.jl

@@ -247,7 +247,7 @@ Input parameters:
                                                                eigenvalues; occupation_threshold),
                                     q=zero(Vec3{real(T)}),
                                     maxiter_sternheimer=100,
-                                    maxiter=100, krylovdim=10, s=1.0, debug_use_old_solver=true,
+                                    maxiter=100, krylovdim=10, s=10.0,
                                     callback=verbose ? OmegaPlusKDefaultCallback() : identity,
                                     kwargs...) where {T}
     # Using χ04P = -Ω^-1, E extension operator (2P->4P) and R restriction operator:
@@ -280,35 +280,18 @@ Input parameters:
                          occupation_threshold, q)
     end
 
-    if debug_use_old_solver
-        # compute total δρ
-        function dielectric_adjoint(δρ)
-            δV = apply_kernel(basis, δρ; ρ, q)
-            # TODO
-            # Would be nice to play with abstol / reltol etc. to avoid over-solving
-            # for the initial GMRES steps.
-            χ0δV = apply_χ0(ham, ψ, occupation, εF, eigenvalues, δV; miniter=1,
-                            occupation_threshold, tol, bandtolalg, q, kwargs...).δρ
-            δρ - χ0δV
-        end
-        δρ, info_gmres = linsolve(dielectric_adjoint, δρ0;
-                                  ishermitian=false, krylovdim,
-                                  tol, verbosity=(verbose ? 3 : 0))
-        info_gmres.converged == 0 && @warn "Solve_ΩplusK_split solver not converged"
-    else
-        # compute total δρ
-        # TODO Can be smarter here, e.g. use mixing to come up with initial guess.
-        ε = DielectricAdjoint(ham, ρ, ψ, occupation, εF, eigenvalues, occupation_threshold,
-                              bandtolalg, maxiter_sternheimer, q)
-        precon = I  # TODO Use mixing
-        callback_inner(info) = callback(merge(info, (; runtime_ns=time_ns() - start_ns)))
-        info_gmres = inexact_gmres(ε, vec(δρ0);
-                                   tol, precon, krylovdim, maxiter, s,
-                                   callback=callback_inner, kwargs...)
-        δρ = reshape(info_gmres.x, size(ρ))
-        if !info_gmres.converged
-            @warn "Solve_ΩplusK_split solver not converged"
-        end
+    # compute total δρ
+    # TODO Can be smarter here, e.g. use mixing to come up with initial guess.
+    ε = DielectricAdjoint(ham, ρ, ψ, occupation, εF, eigenvalues, occupation_threshold,
+                          bandtolalg, maxiter_sternheimer, q)
+    precon = I  # TODO Use mixing
+    callback_inner(info) = callback(merge(info, (; runtime_ns=time_ns() - start_ns)))
+    info_gmres = inexact_gmres(ε, vec(δρ0);
+                               tol, precon, krylovdim, maxiter, s,
+                               callback=callback_inner, kwargs...)
+    δρ = reshape(info_gmres.x, size(ρ))
+    if !info_gmres.converged
+        @warn "Solve_ΩplusK_split solver not converged"
     end
 
     # Compute total change in Hamiltonian applied to ψ
@@ -335,7 +318,8 @@ Input parameters:
     callback((; stage=:final, runtime_ns=time_ns() - start_ns,
                 Axinfo=(; basis, tol=tol*factor_final, resfinal...)))
 
-    (; resfinal.δψ, δρ, δHψ, δVind, δeigenvalues, resfinal.δoccupation, resfinal.δεF, info_gmres)
+    (; resfinal.δψ, δρ, δHψ, δVind, δρ0, δeigenvalues, resfinal.δoccupation, resfinal.δεF,
+     ε, info_gmres)
 end
 
 function solve_ΩplusK_split(scfres::NamedTuple, rhs; kwargs...)
@@ -360,13 +344,17 @@ struct DielectricAdjoint{Tρ, Tψ, Toccupation, TεF, Teigenvalues, Tq}
     maxiter::Int  # CG maximum number of iterations
     q::Tq
 end
-function inexact_mul(ε::DielectricAdjoint, δρ; tol=0.0)
+function DielectricAdjoint(scfres; bandtolalg=BandtolBalanced(scfres), q=zero(Vec3{Float64}), maxiter=100)
+    DielectricAdjoint(scfres.ham, scfres.ρ, scfres.ψ, scfres.occupation, scfres.εF,
+                      scfres.eigenvalues, scfres.occupation_threshold, bandtolalg, maxiter, q)
+end
+function inexact_mul(ε::DielectricAdjoint, δρ; tol=0.0, kwargs...)
     δρ = reshape(δρ, size(ε.ρ))
     basis = ε.ham.basis
     δV = apply_kernel(basis, δρ; ε.ρ, ε.q)
     res = apply_χ0(ε.ham, ε.ψ, ε.occupation, ε.εF, ε.eigenvalues, δV;
                    miniter=1, ε.occupation_threshold, tol,
-                   ε.bandtolalg, ε.q, ε.maxiter)
+                   ε.bandtolalg, ε.q, ε.maxiter, kwargs...)
     χ0δV = res.δρ
     Ax = vec(δρ - χ0δV)  # (1 - χ0 K δρ)
     (; Ax, info=(; tol, basis, res...))