diff --git a/src/block_krylov_solvers.jl b/src/block_krylov_solvers.jl index a20c967e5..070070ac5 100644 --- a/src/block_krylov_solvers.jl +++ b/src/block_krylov_solvers.jl @@ -31,6 +31,7 @@ mutable struct BlockMinresSolver{T,FC,SV,SM} <: BlockKrylovSolver{T,FC,SV,SM} C :: SM D :: SM τ :: SV + tmp :: SM warm_start :: Bool stats :: SimpleStats{T} end @@ -46,8 +47,9 @@ function BlockMinresSolver(m, n, p, SV, SM) C = SM(undef, p, p) D = SM(undef, 2p, p) τ = SV(undef, p) + tmp = C isa Matrix ? SM(undef, 0, 0) : SM(undef, p, p) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") - solver = BlockMinresSolver{T,FC,SV,SM}(m, n, p, ΔX, X, W, P, Q, C, D, τ, false, stats) + solver = BlockMinresSolver{T,FC,SV,SM}(m, n, p, ΔX, X, W, P, Q, C, D, τ, tmp, false, stats) return solver end diff --git a/src/block_minres.jl b/src/block_minres.jl index 004472003..ad143515c 100644 --- a/src/block_minres.jl +++ b/src/block_minres.jl @@ -23,12 +23,12 @@ Solve the Hermitian linear system AX = B of size n with p right-hand sides using #### Input arguments -* `A`: a linear operator that models a Hermitian matrix of dimension n; -* `B`: a matrix of size n × p. +* `A`: a linear operator that models a Hermitian matrix of dimension `n`; +* `B`: a matrix of size `n × p`. #### Optional argument -* `X0`: a matrix of size n × p that represents an initial guess of the solution X. +* `X0`: a matrix of size `n × p` that represents an initial guess of the solution `X`. #### Keyword arguments @@ -45,7 +45,7 @@ Solve the Hermitian linear system AX = B of size n with p right-hand sides using #### Output arguments -* `X`: a dense matrix of size n × p; +* `X`: a dense matrix of size `n × p`; * `stats`: statistics collected on the run in a [`SimpleStats`](@ref) structure. """ function block_minres end @@ -83,27 +83,6 @@ optargs_block_minres = (:X0,) kwargs_block_minres = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) @eval begin - function block_minres($(def_args_block_minres...), $(def_optargs_block_minres...); $(def_kwargs_block_minres...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}} - start_time = time_ns() - solver = BlockMinresSolver(A, B) - warm_start!(solver, $(optargs_block_minres...)) - elapsed_time = ktimer(start_time) - timemax -= elapsed_time - block_minres!(solver, $(args_block_minres...); $(kwargs_block_minres...)) - solver.stats.timer += elapsed_time - return solver.X, solver.stats - end - - function block_minres($(def_args_block_minres...); $(def_kwargs_block_minres...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}} - start_time = time_ns() - solver = BlockMinresSolver(A, B) - elapsed_time = ktimer(start_time) - timemax -= elapsed_time - block_minres!(solver, $(args_block_minres...); $(kwargs_block_minres...)) - solver.stats.timer += elapsed_time - return solver.X, solver.stats - end - function block_minres!(solver :: BlockMinresSolver{T,FC,SV,SM}, $(def_args_block_minres...); $(def_kwargs_block_minres...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, SV <: AbstractVector{FC}, SM <: AbstractMatrix{FC}} # Timer @@ -174,39 +153,87 @@ kwargs_block_minres = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :his # Initial Ψ₁ and V₁ copyto!(V, R₀) - householder!(V, Z, τ) + if C isa Matrix + householder!(V, Z, τ) + else + householder!(V, Z, τ, solver.tmp) + end # Continue the block-Lanczos process. mul!(W, A, V) # Q ← AVₖ for i = 1 : inner_iter - mul!(Ω, V', W) # Ωₖ = Vₖᴴ * Q - # (iter ≥ 2) && mul!(Q, ...) # Q ← Q - βₖ * Vₖ₋₁ * Ψₖᴴ - mul!(Q, V, R, α, β) # Q = Q - Vₖ * Ωₖ + mul!(Ω, V', W) # Ωₖ = Vₖᴴ * Q + (iter ≥ 2) && mul!(Q, Vold, Ψ') # Q ← Q - Vₖ₋₁ * Ψₖᴴ + mul!(Q, V, R, α, β) # Q = Q - Vₖ * Ωₖ end # Vₖ₊₁ and Ψₖ₊₁ are stored in Q and C. - householder!(Q, C, τ) + if C isa Matrix + householder!(Q, C, τ) + else + householder!(Q, C, τ, solver.tmp) + end # Update the QR factorization of Tₖ₊₁.ₖ. - # Apply previous Householder reflections Ωᵢ. - for i = 1 : inner_iter-1 - D1 .= R[nr+i] - D2 .= R[nr+i+1] - @kormqr!('L', trans, H[i], τ[i], D) - R[nr+i] .= D1 - R[nr+i+1] .= D2 + # Apply previous Householder reflections Θₖ₋₂. + if k ≥ 3 + D1 .= Rₖ₋₂.ₖ + D2 .= Rₖ₋₁.ₖ + kormqr!('L', trans, H[i-2], τ[i-2], D) + Rₖ₋₂.ₖ .= D1 + Rₖ₋₁.ₖ .= D2 end - # Compute and apply current Householder reflection Ωₖ. - H[inner_iter][1:p,:] .= R[nr+inner_iter] + # Apply previous Householder reflections Θₖ₋₁. + if k ≥ 2 + D1 .= Rₖ₋₁.ₖ + D2 .= Rₖ.ₖ + kormqr!('L', trans, H[i-1], τ[i-1], D) + Rₖ₋₁.ₖ .= D1 + Rₖ.ₖ .= D2 + end + + # Compute and apply current Householder reflection θₖ. + H[inner_iter][1:p,:] .= Rₖ.ₖ H[inner_iter][p+1:2p,:] .= C - householder!(H[inner_iter], R[nr+inner_iter], τ[inner_iter], compact=true) + if C isa Matrix + householder!(H[i], Rₖ.ₖ, τ[i], compact=true) + else + householder!(H[i], Rₖ.ₖ, τ[i], solver.tmp, compact=true) + end # Update Zₖ = (Qₖ)ᴴΓE₁ = (Λ₁, ..., Λₖ, Λbarₖ₊₁) - D1 .= Z[inner_iter] + D1 .= Λbarₖ D2 .= zero(FC) - @kormqr!('L', trans, H[inner_iter], τ[inner_iter], D) - Z[inner_iter] .= D1 + kormqr!('L', trans, H[i], τ[i], D) + Λₖ .= D1 + + # Compute the directions Wₖ, the last columns of Wₖ = Vₖ(Rₖ)⁻¹ ⟷ (Rₖ)ᵀ(Wₖ)ᵀ = (Vₖ)ᵀ + # R₁₁w₁ = v₁ + if iter == 1 + # wₖ = wₖ₋₁ + # kaxpy!(n, one(FC), uₖ, wₖ) + # wₖ .= wₖ ./ δₖ + end + # R₂₂w₂ = (v₂ - R₂₁w₁) + if iter == 2 + # wₖ = wₖ₋₂ + # kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ) + # kaxpy!(n, one(FC), uₖ, wₖ) + # wₖ .= wₖ ./ δₖ + end + # Rₖₖwₖ = (vₖ - Rₖₖ₋₁wₖ₋₁ - Rₖₖ₋₂wₖ₋₂) + if iter ≥ 3 + # kscal!(n, -ϵₖ₋₂, wₖ₋₂) + # wₖ = wₖ₋₂ + # kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ) + # kaxpy!(n, one(FC), uₖ, wₖ) + # wₖ .= wₖ ./ δₖ + end + + # Update Xₖ = VₖYₖ = WₖZₖ + # Xₖ = Xₖ₋₁ + Λₖ * wₖ + mul!(X, Λₖ, W[i], γ, β) # Update residual norm estimate. # ‖ M(B - AXₖ) ‖_F = ‖Λbarₖ₊₁‖_F diff --git a/src/krylov_solve.jl b/src/krylov_solve.jl index 5cf8b274a..d928b6450 100644 --- a/src/krylov_solve.jl +++ b/src/krylov_solve.jl @@ -48,7 +48,6 @@ for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs (:TrilqrSolver , :trilqr , args_trilqr , def_args_trilqr , optargs_trilqr , def_optargs_trilqr , kwargs_trilqr , def_kwargs_trilqr ) (:CrmrSolver , :crmr , args_crmr , def_args_crmr , () , () , kwargs_crmr , def_kwargs_crmr ) (:CgSolver , :cg , args_cg , def_args_cg , optargs_cg , def_optargs_cg , kwargs_cg , def_kwargs_cg ) - (:CgLanczosShiftSolver, :cg_lanczos_shift, args_cg_lanczos_shift, def_args_cg_lanczos_shift, () , () , kwargs_cg_lanczos_shift, def_kwargs_cg_lanczos_shift) (:CglsSolver , :cgls , args_cgls , def_args_cgls , () , () , kwargs_cgls , def_kwargs_cgls ) (:CgLanczosSolver , :cg_lanczos , args_cg_lanczos , def_args_cg_lanczos , optargs_cg_lanczos, def_optargs_cg_lanczos, kwargs_cg_lanczos , def_kwargs_cg_lanczos ) (:BilqSolver , :bilq , args_bilq , def_args_bilq , optargs_bilq , def_optargs_bilq , kwargs_bilq , def_kwargs_bilq ) @@ -58,6 +57,8 @@ for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs (:FgmresSolver , :fgmres , args_fgmres , def_args_fgmres , optargs_fgmres , def_optargs_fgmres , kwargs_fgmres , def_kwargs_fgmres ) (:FomSolver , :fom , args_fom , def_args_fom , optargs_fom , def_optargs_fom , kwargs_fom , def_kwargs_fom ) (:GpmrSolver , :gpmr , args_gpmr , def_args_gpmr , optargs_gpmr , def_optargs_gpmr , kwargs_gpmr , def_kwargs_gpmr ) + (:CgLanczosShiftSolver , :cg_lanczos_shift , args_cg_lanczos_shift , def_args_cg_lanczos_shift , (), (), kwargs_cg_lanczos_shift , def_kwargs_cg_lanczos_shift ) + (:CglsLanczosShiftSolver, :cgls_lanczos_shift, args_cgls_lanczos_shift, def_args_cgls_lanczos_shift, (), (), kwargs_cgls_lanczos_shift, def_kwargs_cgls_lanczos_shift) ] # Create the symbol for the in-place method krylov! = Symbol(krylov, :!)