From 0138dfb1d13d96dcf6998bd033d1f92cadcc7af7 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Sun, 24 Sep 2023 12:42:21 -0500 Subject: [PATCH] Update Krylov processes to return more arguments --- docs/src/processes.md | 54 +++++++++++++++++------------------ src/krylov_processes.jl | 63 +++++++++++++++++++++++++---------------- test/test_processes.jl | 62 +++++++++++++++++++++++++++++----------- 3 files changed, 111 insertions(+), 68 deletions(-) diff --git a/docs/src/processes.md b/docs/src/processes.md index 2da69d04e..7452b0577 100644 --- a/docs/src/processes.md +++ b/docs/src/processes.md @@ -60,18 +60,18 @@ For matrices $C \in \mathbb{C}^{n \times n} \enspace$ and $\enspace T \in \mathb \left\{\sum_{i=0}^{k-1} C^i T \, \Omega_i \, \middle \vert \, \Omega_i \in \mathbb{C}^{p \times p},~0 \le i \le k-1 \right\}. ``` -## Hermitian Lanczos +## [Hermitian Lanczos](@id hermitian-lanczos) ![hermitian_lanczos](./graphics/hermitian_lanczos.png) After $k$ iterations of the Hermitian Lanczos process, the situation may be summarized as ```math \begin{align*} - A V_k &= V_k T_k + \beta_{k+1,k} v_{k+1} e_k^T = V_{k+1} T_{k+1,k}, \\ + A V_k &= V_k T_k + \beta_{k+1,k} v_{k+1} e_k^T = V_{k+1} T_{k+1,k}, \\ V_k^H V_k &= I_k, \end{align*} ``` -where $V_k$ is an orthonormal basis of the Krylov subspace $\mathcal{K}_k (A,b)$, +where $V_k$ is an orthonormal basis of the Krylov subspace $\mathcal{K}_k(A,b)$, ```math T_k = \begin{bmatrix} @@ -89,22 +89,22 @@ T_{k+1,k} = ``` Note that depending on how we normalize the vectors that compose $V_k$, $T_{k+1,k}$ can be a real tridiagonal matrix even if $A$ is a complex matrix. -The function [`hermitian_lanczos`](@ref hermitian_lanczos) returns $V_{k+1}$ and $T_{k+1,k}$. +The function [`hermitian_lanczos`](@ref hermitian_lanczos(::Any, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $\beta_1$ and $T_{k+1,k}$. Related methods: [`SYMMLQ`](@ref symmlq), [`CG`](@ref cg), [`CR`](@ref cr), [`CAR`](@ref car), [`MINRES`](@ref minres), [`MINRES-QLP`](@ref minres_qlp), [`MINARES`](@ref minares), [`CGLS`](@ref cgls), [`CRLS`](@ref crls), [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`CG-LANCZOS`](@ref cg_lanczos) and [`CG-LANCZOS-SHIFT`](@ref cg_lanczos_shift). ```@docs -hermitian_lanczos +hermitian_lanczos(::Any, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat) ``` -## Non-Hermitian Lanczos +## [Non-Hermitian Lanczos](@id nonhermitian-lanczos) ![nonhermitian_lanczos](./graphics/nonhermitian_lanczos.png) After $k$ iterations of the non-Hermitian Lanczos process (also named the Lanczos biorthogonalization process), the situation may be summarized as ```math \begin{align*} - A V_k &= V_k T_k + \beta_{k+1} v_{k+1} e_k^T = V_{k+1} T_{k+1,k}, \\ + A V_k &= V_k T_k + \beta_{k+1} v_{k+1} e_k^T = V_{k+1} T_{k+1,k}, \\ A^H U_k &= U_k T_k^H + \bar{\gamma}_{k+1} u_{k+1} e_k^T = U_{k+1} T_{k,k+1}^H, \\ V_k^H U_k &= U_k^H V_k = I_k, \end{align*} @@ -131,7 +131,7 @@ T_{k,k+1} = \end{bmatrix}. ``` -The function [`nonhermitian_lanczos`](@ref nonhermitian_lanczos) returns $V_{k+1}$, $T_{k+1,k}$, $U_{k+1}$ and $T_{k,k+1}^H$. +The function [`nonhermitian_lanczos`](@ref nonhermitian_lanczos(::Any, ::AbstractVector{FC}, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $\beta_1$, $T_{k+1,k}$, $U_{k+1}$, $\bar{\gamma}_1$ and $T_{k,k+1}^H$. Related methods: [`BiLQ`](@ref bilq), [`QMR`](@ref qmr), [`BiLQR`](@ref bilqr), [`CGS`](@ref cgs) and [`BICGSTAB`](@ref bicgstab). @@ -140,7 +140,7 @@ Related methods: [`BiLQ`](@ref bilq), [`QMR`](@ref qmr), [`BiLQR`](@ref bilqr), With these scaling factors, the non-Hermitian Lanczos process coincides with the Hermitian Lanczos process when $A = A^H$ and $b = c$. ```@docs -nonhermitian_lanczos +nonhermitian_lanczos(::Any, ::AbstractVector{FC}, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat) ``` The non-Hermitian Lanczos process can be also implemented without $A^H$ (transpose-free variant). @@ -158,7 +158,7 @@ The polynomials are defined from the recursions: \end{align*} ``` -Because $\alpha_k = u_k^H A v_k$ and $(\bar{\gamma}_{k+1} u_{k+1})^H (\beta_{k+1} v_{k+1}) = \gamma_{k+1} \beta_{k+1}$, we can define the coefficients of $T_{k+1,k}$ and $T_{k,k+1}^H$ as follows: +Because $\alpha_k = u_k^H A v_k$ and $(\bar{\gamma}_{k+1} u_{k+1})^H (\beta_{k+1} v_{k+1}) = \gamma_{k+1} \beta_{k+1}$, we can determine the coefficients of $T_{k+1,k}$ and $T_{k,k+1}^H$ as follows: ```math \begin{align*} \alpha_k &= \dfrac{1}{\gamma_k \beta_k} \langle~Q_{k-1}(A^H) c \, , \, A P_{k-1}(A) b~\rangle \\ @@ -169,7 +169,7 @@ Because $\alpha_k = u_k^H A v_k$ and $(\bar{\gamma}_{k+1} u_{k+1})^H (\beta_{k+1 \end{align*} ``` -## Arnoldi +## [Arnoldi](@id arnoldi) ![arnoldi](./graphics/arnoldi.png) @@ -197,7 +197,7 @@ H_{k+1,k} = \end{bmatrix}. ``` -The function [`arnoldi`](@ref arnoldi) returns $V_{k+1}$ and $H_{k+1,k}$. +The function [`arnoldi`](@ref arnoldi(::Any, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $\beta$ and $H_{k+1,k}$. Related methods: [`DIOM`](@ref diom), [`FOM`](@ref fom), [`DQGMRES`](@ref dqgmres), [`GMRES`](@ref gmres) and [`FGMRES`](@ref fgmres). @@ -205,17 +205,17 @@ Related methods: [`DIOM`](@ref diom), [`FOM`](@ref fom), [`DQGMRES`](@ref dqgmre The Arnoldi process coincides with the Hermitian Lanczos process when $A$ is Hermitian. ```@docs -arnoldi +arnoldi(::Any, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat) ``` -## Golub-Kahan +## [Golub-Kahan](@id golub-kahan) ![golub_kahan](./graphics/golub_kahan.png) After $k$ iterations of the Golub-Kahan bidiagonalization process, the situation may be summarized as ```math \begin{align*} - A V_k &= U_{k+1} B_k, \\ + A V_k &= U_{k+1} B_k, \\ A^H U_{k+1} &= V_k B_k^H + \bar{\alpha}_{k+1} v_{k+1} e_{k+1}^T = V_{k+1} L_{k+1}^H, \\ V_k^H V_k &= U_k^H U_k = I_k, \end{align*} @@ -246,7 +246,7 @@ B_k = ``` Note that depending on how we normalize the vectors that compose $V_k$ and $U_k$, $L_k$ can be a real bidiagonal matrix even if $A$ is a complex matrix. -The function [`golub_kahan`](@ref golub_kahan) returns $V_{k+1}$, $U_{k+1}$ and $L_{k+1}$. +The function [`golub_kahan`](@ref golub_kahan(::Any, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $U_{k+1}$, $\beta_1$ and $L_{k+1}$. Related methods: [`LNLQ`](@ref lnlq), [`CRAIG`](@ref craig), [`CRAIGMR`](@ref craigmr), [`LSLQ`](@ref lslq), [`LSQR`](@ref lsqr) and [`LSMR`](@ref lsmr). @@ -255,17 +255,17 @@ Related methods: [`LNLQ`](@ref lnlq), [`CRAIG`](@ref craig), [`CRAIGMR`](@ref cr It is also related to the Hermitian Lanczos process applied to $\begin{bmatrix} 0 & A \\ A^H & 0 \end{bmatrix}$ with initial vector $\begin{bmatrix} b \\ 0 \end{bmatrix}$. ```@docs -golub_kahan +golub_kahan(::Any, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat) ``` -## Saunders-Simon-Yip +## [Saunders-Simon-Yip](@id saunders-simon-yip) ![saunders_simon_yip](./graphics/saunders_simon_yip.png) After $k$ iterations of the Saunders-Simon-Yip process (also named the orthogonal tridiagonalization process), the situation may be summarized as ```math \begin{align*} - A U_k &= V_k T_k + \beta_{k+1} v_{k+1} e_k^T = V_{k+1} T_{k+1,k}, \\ + A U_k &= V_k T_k + \beta_{k+1} v_{k+1} e_k^T = V_{k+1} T_{k+1,k}, \\ A^H V_k &= U_k T_k^H + \bar{\gamma}_{k+1} u_{k+1} e_k^T = U_{k+1} T_{k,k+1}^H, \\ V_k^H V_k &= U_k^H U_k = I_k, \end{align*} @@ -292,18 +292,18 @@ T_{k,k+1} = \end{bmatrix}. ``` -The function [`saunders_simon_yip`](@ref saunders_simon_yip) returns $V_{k+1}$, $T_{k+1,k}$, $U_{k+1}$ and $T_{k,k+1}^H$. +The function [`saunders_simon_yip`](@ref saunders_simon_yip(::Any, ::AbstractVector{FC}, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $\beta_1$, $T_{k+1,k}$, $U_{k+1}$, $\bar{\gamma}_1$ and $T_{k,k+1}^H$. Related methods: [`USYMLQ`](@ref usymlq), [`USYMQR`](@ref usymqr), [`TriLQR`](@ref trilqr), [`TriCG`](@ref tricg) and [`TriMR`](@ref trimr). -```@docs -saunders_simon_yip -``` - !!! note The Saunders-Simon-Yip is equivalent to the block-Lanczos process applied to $\begin{bmatrix} 0 & A \\ A^H & 0 \end{bmatrix}$ with initial matrix $\begin{bmatrix} b & 0 \\ 0 & c \end{bmatrix}$. -## Montoison-Orban +```@docs +saunders_simon_yip(::Any, ::AbstractVector{FC}, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat) +``` + +## [Montoison-Orban](@id montoison-orban) ![montoison_orban](./graphics/montoison_orban.png) @@ -347,7 +347,7 @@ F_{k+1,k} = \end{bmatrix}. ``` -The function [`montoison_orban`](@ref montoison_orban) returns $V_{k+1}$, $H_{k+1,k}$, $U_{k+1}$ and $F_{k+1,k}$. +The function [`montoison_orban`](@ref montoison_orban(::Any, ::Any, ::AbstractVector{FC}, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $\beta$, $H_{k+1,k}$, $U_{k+1}$, $\gamma$ and $F_{k+1,k}$. Related methods: [`GPMR`](@ref gpmr). @@ -356,5 +356,5 @@ Related methods: [`GPMR`](@ref gpmr). It also coincides with the Saunders-Simon-Yip process when $B = A^H$. ```@docs -montoison_orban +montoison_orban(::Any, ::Any, ::AbstractVector{FC}, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat) ``` diff --git a/src/krylov_processes.jl b/src/krylov_processes.jl index d482e04eb..ebcf7f513 100644 --- a/src/krylov_processes.jl +++ b/src/krylov_processes.jl @@ -1,7 +1,7 @@ export hermitian_lanczos, nonhermitian_lanczos, arnoldi, golub_kahan, saunders_simon_yip, montoison_orban """ - V, T = hermitian_lanczos(A, b, k) + V, β, T = hermitian_lanczos(A, b, k) #### Input arguments @@ -12,6 +12,7 @@ export hermitian_lanczos, nonhermitian_lanczos, arnoldi, golub_kahan, saunders_s #### Output arguments * `V`: a dense n × (k+1) matrix; +* `β`: a coefficient such that βv₁ = b; * `T`: a sparse (k+1) × k tridiagonal matrix. #### Reference @@ -41,6 +42,7 @@ function hermitian_lanczos(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOr end end + local β₁ V = M(undef, n, k+1) T = SparseMatrixCSC(k+1, k, colptr, rowval, nzval) @@ -49,8 +51,8 @@ function hermitian_lanczos(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOr vᵢ = view(V,:,i) vᵢ₊₁ = q = view(V,:,i+1) if i == 1 - βᵢ = @knrm2(n, b) - vᵢ .= b ./ βᵢ + β₁ = @knrm2(n, b) + vᵢ .= b ./ β₁ end mul!(q, A, vᵢ) if i ≥ 2 @@ -67,11 +69,11 @@ function hermitian_lanczos(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOr vᵢ₊₁ .= q ./ βᵢ₊₁ pαᵢ = pαᵢ + 3 end - return V, T + return V, β₁, T end """ - V, T, U, Tᴴ = nonhermitian_lanczos(A, b, c, k) + V, β, T, U, γᴴ, Tᴴ = nonhermitian_lanczos(A, b, c, k) #### Input arguments @@ -83,8 +85,10 @@ end #### Output arguments * `V`: a dense n × (k+1) matrix; +* `β`: a coefficient such that βv₁ = b; * `T`: a sparse (k+1) × k tridiagonal matrix; * `U`: a dense n × (k+1) matrix; +* `γᴴ`: a coefficient such that γᴴu₁ = c; * `Tᴴ`: a sparse (k+1) × k tridiagonal matrix. #### Reference @@ -115,6 +119,7 @@ function nonhermitian_lanczos(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k end end + local β₁, γ₁ᴴ V = M(undef, n, k+1) U = M(undef, n, k+1) T = SparseMatrixCSC(k+1, k, colptr, rowval, nzval_T) @@ -128,10 +133,10 @@ function nonhermitian_lanczos(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k uᵢ₊₁ = p = view(U,:,i+1) if i == 1 cᴴb = @kdot(n, c, b) - βᵢ = √(abs(cᴴb)) - γᵢ = cᴴb / βᵢ - vᵢ .= b ./ βᵢ - uᵢ .= c ./ conj(γᵢ) + β₁ = √(abs(cᴴb)) + γ₁ᴴ = conj(cᴴb / β₁) + vᵢ .= b ./ β₁ + uᵢ .= c ./ γ₁ᴴ end mul!(q, A , vᵢ) mul!(p, Aᴴ, uᵢ) @@ -161,11 +166,11 @@ function nonhermitian_lanczos(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k end pαᵢ = pαᵢ + 3 end - return V, T, U, Tᴴ + return V, β₁, T, U, γ₁ᴴ, Tᴴ end """ - V, H = arnoldi(A, b, k; reorthogonalization=false) + V, β, H = arnoldi(A, b, k; reorthogonalization=false) #### Input arguments @@ -180,6 +185,7 @@ end #### Output arguments * `V`: a dense n × (k+1) matrix; +* `β`: a coefficient such that βv₁ = b; * `H`: a dense (k+1) × k upper Hessenberg matrix. #### Reference @@ -191,6 +197,7 @@ function arnoldi(A, b::AbstractVector{FC}, k::Int; reorthogonalization::Bool=fal S = ktypeof(b) M = vector_to_matrix(S) + local β V = M(undef, n, k+1) H = zeros(FC, k+1, k) @@ -218,11 +225,11 @@ function arnoldi(A, b::AbstractVector{FC}, k::Int; reorthogonalization::Bool=fal H[j+1,j] = @knrm2(n, q) vⱼ₊₁ .= q ./ H[j+1,j] end - return V, H + return V, β, H end """ - V, U, L = golub_kahan(A, b, k) + V, U, β, L = golub_kahan(A, b, k) #### Input arguments @@ -234,6 +241,7 @@ end * `V`: a dense n × (k+1) matrix; * `U`: a dense m × (k+1) matrix; +* `β`: a coefficient such that βu₁ = b; * `L`: a sparse (k+1) × (k+1) lower bidiagonal matrix. #### References @@ -262,6 +270,7 @@ function golub_kahan(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOrComple rowval[2k+1] = k+1 colptr[k+2] = 2k+2 + local β₁ V = M(undef, n, k+1) U = M(undef, m, k+1) L = SparseMatrixCSC(k+1, k+1, colptr, rowval, nzval) @@ -274,8 +283,8 @@ function golub_kahan(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOrComple vᵢ₊₁ = p = view(V,:,i+1) if i == 1 wᵢ = vᵢ - βᵢ = @knrm2(m, b) - uᵢ .= b ./ βᵢ + β₁ = @knrm2(m, b) + uᵢ .= b ./ β₁ mul!(wᵢ, Aᴴ, uᵢ) αᵢ = @knrm2(n, wᵢ) nzval[pαᵢ] = αᵢ # Lᵢ.ᵢ = αᵢ @@ -294,11 +303,11 @@ function golub_kahan(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOrComple nzval[pαᵢ+2] = αᵢ₊₁ # Lᵢ₊₁.ᵢ₊₁ = αᵢ₊₁ pαᵢ = pαᵢ + 2 end - return V, U, L + return V, U, β₁, L end """ - V, T, U, Tᴴ = saunders_simon_yip(A, b, c, k) + V, β, T, U, γᴴ, Tᴴ = saunders_simon_yip(A, b, c, k) #### Input arguments @@ -310,8 +319,10 @@ end #### Output arguments * `V`: a dense m × (k+1) matrix; +* `β`: a coefficient such that βv₁ = b; * `T`: a sparse (k+1) × k tridiagonal matrix; * `U`: a dense n × (k+1) matrix; +* `γᴴ`: a coefficient such that γᴴu₁ = c; * `Tᴴ`: a sparse (k+1) × k tridiagonal matrix. #### Reference @@ -342,6 +353,7 @@ function saunders_simon_yip(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k:: end end + local β₁, γ₁ᴴ V = M(undef, m, k+1) U = M(undef, n, k+1) T = SparseMatrixCSC(k+1, k, colptr, rowval, nzval_T) @@ -354,10 +366,10 @@ function saunders_simon_yip(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k:: vᵢ₊₁ = q = view(V,:,i+1) uᵢ₊₁ = p = view(U,:,i+1) if i == 1 - β = @knrm2(m, b) - γ = @knrm2(n, c) - vᵢ .= b ./ β - uᵢ .= c ./ γ + β₁ = @knrm2(m, b) + γ₁ᴴ = @knrm2(n, c) + vᵢ .= b ./ β₁ + uᵢ .= c ./ γ₁ᴴ end mul!(q, A , uᵢ) mul!(p, Aᴴ, vᵢ) @@ -386,11 +398,11 @@ function saunders_simon_yip(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k:: end pαᵢ = pαᵢ + 3 end - return V, T, U, Tᴴ + return V, β₁, T, U, γ₁ᴴ, Tᴴ end """ - V, H, U, F = montoison_orban(A, B, b, c, k; reorthogonalization=false) + V, β, H, U, γ, F = montoison_orban(A, B, b, c, k; reorthogonalization=false) #### Input arguments @@ -407,8 +419,10 @@ end #### Output arguments * `V`: a dense m × (k+1) matrix; +* `β`: a coefficient such that βv₁ = b; * `H`: a dense (k+1) × k upper Hessenberg matrix; * `U`: a dense n × (k+1) matrix; +* `γ`: a coefficient such that γu₁ = c; * `F`: a dense (k+1) × k upper Hessenberg matrix. #### Reference @@ -420,6 +434,7 @@ function montoison_orban(A, B, b::AbstractVector{FC}, c::AbstractVector{FC}, k:: S = ktypeof(b) M = vector_to_matrix(S) + local β, γ V = M(undef, m, k+1) U = M(undef, n, k+1) H = zeros(FC, k+1, k) @@ -463,5 +478,5 @@ function montoison_orban(A, B, b::AbstractVector{FC}, c::AbstractVector{FC}, k:: F[j+1,j] = @knrm2(n, p) uⱼ₊₁ .= p ./ F[j+1,j] end - return V, H, U, F + return V, β, H, U, γ, F end diff --git a/test/test_processes.jl b/test/test_processes.jl index 49b2bcd07..237f70ce5 100644 --- a/test/test_processes.jl +++ b/test/test_processes.jl @@ -19,6 +19,7 @@ end m = 250 n = 500 k = 20 + s = 5 for FC in (Float64, ComplexF64) R = real(FC) @@ -29,9 +30,13 @@ end @testset "Data Type: $FC" begin @testset "Hermitian Lanczos" begin - A, b = symmetric_indefinite(n, FC=FC) - V, T = hermitian_lanczos(A, b, k) + A = rand(FC, n, n) + A = A' * A + b = rand(FC, n) + V, β₁, T = hermitian_lanczos(A, b, k) + @test norm(V[:,1:s]' * V[:,1:s] - I) ≤ 1e-4 + @test β₁ * V[:,1] ≈ b @test A * V[:,1:k] ≈ V * T storage_hermitian_lanczos_bytes(n, k) = 4k * nbits_I + (3k-1) * nbits_R + n*(k+1) * nbits_FC @@ -42,10 +47,15 @@ end end @testset "Non-Hermitian Lanczos" begin - A, b = nonsymmetric_definite(n, FC=FC) - c = -b - V, T, U, Tᴴ = nonhermitian_lanczos(A, b, c, k) - + A = rand(FC, n, n) + b = rand(FC, n) + c = rand(FC, n) + V, β₁, T, U, γ₁ᴴ, Tᴴ = nonhermitian_lanczos(A, b, c, k) + + @test norm(V[:,1:s]' * U[:,1:s] - I) ≤ 1e-4 + @test norm(U[:,1:s]' * V[:,1:s] - I) ≤ 1e-4 + @test β₁ * V[:,1] ≈ b + @test γ₁ᴴ * U[:,1] ≈ c @test T[1:k,1:k] ≈ Tᴴ[1:k,1:k]' @test A * V[:,1:k] ≈ V * T @test A' * U[:,1:k] ≈ U * Tᴴ @@ -58,10 +68,13 @@ end end @testset "Arnoldi" begin + A = rand(FC, n, n) + b = rand(FC, n) for reorthogonalization in (false, true) - A, b = nonsymmetric_indefinite(n, FC=FC) - V, H = arnoldi(A, b, k; reorthogonalization) + V, β, H = arnoldi(A, b, k; reorthogonalization) + @test norm(V[:,1:s]' * V[:,1:s] - I) ≤ 1e-4 + @test β * V[:,1] ≈ b @test A * V[:,1:k] ≈ V * H function storage_arnoldi_bytes(n, k) @@ -75,10 +88,14 @@ end end @testset "Golub-Kahan" begin - A, b = under_consistent(m, n, FC=FC) - V, U, L = golub_kahan(A, b, k) + A = rand(FC, m, n) + b = rand(FC, m) + V, U, β₁, L = golub_kahan(A, b, k) B = L[1:k+1,1:k] + @test norm(V[:,1:s]' * V[:,1:s] - I) ≤ 1e-4 + @test norm(U[:,1:s]' * U[:,1:s] - I) ≤ 1e-4 + @test β₁ * U[:,1] ≈ b @test A * V[:,1:k] ≈ U * B @test A' * U ≈ V * L' @test A' * A * V[:,1:k] ≈ V * L' * B @@ -92,10 +109,15 @@ end end @testset "Saunders-Simon-Yip" begin - A, b = under_consistent(m, n, FC=FC) - _, c = over_consistent(n, m, FC=FC) - V, T, U, Tᴴ = saunders_simon_yip(A, b, c, k) - + A = rand(FC, m, n) + b = rand(FC, m) + c = rand(FC, n) + V, β₁, T, U, γ₁ᴴ, Tᴴ = saunders_simon_yip(A, b, c, k) + + @test norm(V[:,1:s]' * V[:,1:s] - I) ≤ 1e-4 + @test norm(U[:,1:s]' * U[:,1:s] - I) ≤ 1e-4 + @test β₁ * V[:,1] ≈ b + @test γ₁ᴴ * U[:,1] ≈ c @test T[1:k,1:k] ≈ Tᴴ[1:k,1:k]' @test A * U[:,1:k] ≈ V * T @test A' * V[:,1:k] ≈ U * Tᴴ @@ -118,11 +140,17 @@ end end @testset "Montoison-Orban" begin + A = rand(FC, m, n) + B = rand(FC, n, m) + b = rand(FC, m) + c = rand(FC, n) for reorthogonalization in (false, true) - A, b = under_consistent(m, n, FC=FC) - B, c = over_consistent(n, m, FC=FC) - V, H, U, F = montoison_orban(A, B, b, c, k; reorthogonalization) + V, β, H, U, γ, F = montoison_orban(A, B, b, c, k; reorthogonalization) + @test norm(V[:,1:s]' * V[:,1:s] - I) ≤ 1e-4 + @test norm(U[:,1:s]' * U[:,1:s] - I) ≤ 1e-4 + @test β * V[:,1] ≈ b + @test γ * U[:,1] ≈ c @test A * U[:,1:k] ≈ V * H @test B * V[:,1:k] ≈ U * F @test B * A * U[:,1:k-1] ≈ U * F * H[1:k,1:k-1]