initial files

Manifest.toml

+# This file is machine-generated - editing it directly is not advised
+
+[[FastLapackInterface]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "a82ae140565456800a09c026f200e7046b247ac7"
+repo-rev = "master"
+repo-url = "https://git.dynare.org/julia-packages/fastlapackinterface.jl.git"
+uuid = "29a986be-02c6-4525-aec4-84b980013641"
+version = "0.1.0"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

Project.toml

+name = "PolynomialMatrixEquations"
+uuid = "4f9d485d-518f-41ed-81c8-372cd804c93b"
+authors = ["michel "]
+version = "0.1.0"
+
+[deps]
+FastLapackInterface = "29a986be-02c6-4525-aec4-84b980013641"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/CyclicReduction.jl

+using FastLapackInterface.LinSolveAlgo: LinSolveWs, linsolve_core!
+
+using LinearAlgebra
+using LinearAlgebra.BLAS: gemm!
+export CyclicReductionWs, cyclic_reduction!, cyclic_reduction_check
+
+mutable struct CyclicReductionWs
+    linsolve_ws::LinSolveWs
+    ahat1::Matrix{Float64}
+    a1copy::Matrix{Float64}
+    m::Matrix{Float64,}
+    m00::SubArray{Float64}
+    m02::SubArray{Float64}
+    m20::SubArray{Float64}
+    m22::SubArray{Float64}
+    m1::Matrix{Float64}
+    m1_a0::SubArray{Float64}
+    m1_a2::SubArray{Float64}
+    m2::Matrix{Float64}
+    m2_a0::SubArray{Float64}
+    m2_a2::SubArray{Float64}
+    info::Int64
+
+    function CyclicReductionWs(n)
+        linsolve_ws = LinSolveWs(n)
+        ahat1 = Matrix{Float64}(undef, n,n)
+	a1copy = Matrix{Float64}(undef, n,n)
+        m = Matrix{Float64}(undef, 2*n,2*n)
+        m00 = view(m,1:n,1:n)
+        m02 = view(m,1:n,n .+ (1:n))
+        m20 = view(m,n .+ (1:n), 1:n)
+        m22 = view(m,n .+ (1:n),n .+(1:n))
+        m1 = Matrix{Float64}(undef, n, 2*n)
+        m1_a0 = view(m1,1:n,1:n)
+        m1_a2 = view(m1,1:n,n .+ (1:n))
+        m2 = Matrix{Float64}(undef, 2*n, n)
+        m2_a0 = view(m2, 1:n, 1:n)
+        m2_a2 = view(m2, n .+ (1:n), 1:n)
+        info = 0
+        new(linsolve_ws,ahat1,a1copy,m, m00, m02, m20, m22, m1, m1_a0, m1_a2, m2, m2_a0, m2_a2, info) 
+    end
+end
+
+"""
+    cyclic_reduction!(x::Array{Float64},a0::Array{Float64},a1::Array{Float64},a2::Array{Float64},ws::CyclicReductionWs, cvg_tol::Float64, max_it::Int64)
+
+Solve the quadratic matrix equation a0 + a1*x + a2*x*x = 0, using the cyclic reduction method from Bini et al. (???).
+
+The solution is returned in matrix x. In case of nonconvergency, x is set to NaN and an error code is returned in ws.info
+* info = 0: return OK
+* info = 1: no stable solution (????)
+* info = 2: multiple stable solutions (????)
+
+# Example
+```meta
+DocTestSetup = quote
+     using CyclicReduction
+     n = 3
+     ws = CyclicReductionWs(n)
+     a0 = [0.5 0 0; 0 0.5 0; 0 0 0];
+     a1 = eye(n)
+     a2 = [0 0 0; 0 0 0; 0 0 0.8]
+     x = zeros(n,n)
+end
+```
+
+```jldoctest
+julia> display(names(CyclicReduction))
+```
+
+```jldoctest
+julia> cyclic_reduction!(x,a0,a1,a2,ws,1e-8,50)
+```
+"""
+function cyclic_reduction!(x::Array{Float64},a0::Array{Float64},a1::Array{Float64},a2::Array{Float64},ws::CyclicReductionWs, cvg_tol::Float64, max_it::Int64)
+    copyto!(x,a0)
+    copyto!(ws.ahat1,1,a1,1,length(a1))
+    @inbounds copyto!(ws.m1_a0, a0)
+    @inbounds copyto!(ws.m1_a2, a2)
+    @inbounds copyto!(ws.m2_a0, a0)
+    @inbounds copyto!(ws.m2_a2, a2)
+    it = 0
+    @inbounds while it < max_it
+        #        ws.m = [a0; a2]*(a1\[a0 a2])
+	copyto!(ws.a1copy,a1)
+        linsolve_core!(ws.a1copy, ws.m1, ws.linsolve_ws)
+        gemm!('N','N',-1.0,ws.m2,ws.m1,0.0,ws.m)
+        @simd for i in eachindex(a1)
+            a1[i] += ws.m02[i] + ws.m20[i]
+        end
+        copyto!(ws.m1_a0, ws.m00)
+        copyto!(ws.m1_a2, ws.m22)
+        copyto!(ws.m2_a0, ws.m00)
+        copyto!(ws.m2_a2, ws.m22)
+        if any(isinf.(ws.m))
+            if norm(ws.m1_a0) < cvg_tol
+                ws.info = 2
+            else
+                ws.info = 1
+            end
+            fill!(x,NaN)
+            return
+        end
+        ws.ahat1 += ws.m20
+        crit = norm(ws.m1_a0,1)
+        if crit < cvg_tol
+	    # keep iterating until condition on a2 is met
+            if norm(ws.m1_a2,1) < cvg_tol
+                break
+            end
+        end
+        it += 1
+    end
+    if it == max_it
+        if norm(ws.m1_a0) < cvg_tol
+            ws.info = 2
+        else
+            ws.info = 1
+        end
+        fill!(x,NaN)
+        return
+    else
+        linsolve_core!(ws.ahat1, x, ws.linsolve_ws)
+        @inbounds lmul!(-1.0,x)
+        ws.info = 0
+    end        
+end
+
+function cyclic_reduction_check(x::Array{Float64,2},a0::Array{Float64,2}, a1::Array{Float64,2}, a2::Array{Float64,2},cvg_tol::Float64)
+    res = a0 + a1*x + a2*x*x
+    if (sum(sum(abs.(res))) > cvg_tol)
+        print("the norm of the residuals, ", res, ", compared to the tolerance criterion ",cvg_tol)
+    end
+    nothing
+end
+

src/GeneralizedSchurDecompositionSolver.jl

+using FastLapackInterface.SchurAlgo: DggesWs, dgges!
+using FastLapackInterface.LinSolveAlgo: LinSolveWs, linsolve_core!
+using LinearAlgebra
+using LinearAlgebra.BLAS
+
+struct GsSolverWs
+    dgges_ws::DggesWs
+    linsolve_ws_11::LinSolveWs
+    linsolve_ws_22::LinSolveWs
+    D11::SubArray{Float64}
+    E11::SubArray{Float64}
+    vsr::Matrix{Float64}
+    Z11::SubArray{Float64}
+    Z12::SubArray{Float64}
+    Z21::SubArray{Float64}
+    Z22::SubArray{Float64}
+    tmp1::Matrix{Float64}
+    tmp2::Matrix{Float64}
+    tmp3::Matrix{Float64}
+    g1::Matrix{Float64}
+    g2::Matrix{Float64}
+    eigval::Vector{ComplexF64}
+    
+    function GsSolverWs(d,e,n1)
+        dgges_ws = DggesWs(Ref{UInt8}('N'), Ref{UInt8}('N'), Ref{UInt8}('N'), e, d)
+        n = size(d,1)
+        n2 = n - n1
+        linsolve_ws_11 = LinSolveWs(n1)
+        linsolve_ws_22 = LinSolveWs(n2)
+        D11 = view(d,1:n1,1:n1)
+        E11 = view(e,1:n1,1:n1)
+        vsr = Matrix{Float64}(undef, n, n)
+        Z11 = view(vsr,1:n1,1:n1)
+        Z12 = view(vsr,1:n1,n1+1:n)
+        Z21 = view(vsr,n1+1:n,1:n1)
+        Z22 = view(vsr,n1+1:n,n1+1:n)
+        tmp1 = Matrix{Float64}(undef, n2, n1)
+        tmp2 = Matrix{Float64}(undef, n1, n1)
+        tmp3 = Matrix{Float64}(undef, n1, n1)
+        g1 = Matrix{Float64}(undef, n1, n1)
+        g2 = Matrix{Float64}(undef, n2, n1)
+        eigval = Vector{ComplexF64}(undef, n)
+        new(dgges_ws,linsolve_ws_11, linsolve_ws_22, D11,E11,vsr,Z11,Z12,Z21,Z22,tmp1,tmp2,tmp3,g1,g2, eigval)
+    end
+end
+
+struct UnstableSystemException <: Exception end
+struct UndeterminateSystemException <: Exception end
+
+#"""
+#    gs_solver!(ws::GsSolverWs,d::Matrix{Float64},e::Matrix{Float64},n1::Int64,qz_criterium)
+#
+#finds the unique stable solution for the following system:
+#
+#```
+#d \left[\begin{array}{c}I\\g_2\end{array}\right]g_1 = e \left[\begin{array}{c}I\\g_2\end{array}\right]
+#```
+#The solution is returned in ``ws.g1`` and ``ws.g2``
+#"""
+function gs_solver!(ws::GsSolverWs,d::Matrix{Float64},e::Matrix{Float64},n1::Int64,qz_criterium::Float64)
+
+    dgges!('N', 'V', e, d, zeros(1,1), ws.vsr, ws.eigval, ws.dgges_ws)
+    nstable = ws.dgges_ws.sdim[]
+    
+    if nstable < n1
+        throw(UnstableSystemException)
+    elseif nstable > n1
+        throw(UndeterminateSystemExcpetion)
+    end
+    ws.g2 .= ws.Z12'
+    linsolve_core!(ws.Z22', ws.g2, ws.linsolve_ws_22)
+    lmul!(-1.0,ws.g2)
+    ws.tmp2 .= ws.Z11'
+    ws.D11 .= view(d,1:nstable,1:nstable)
+    linsolve_core!(ws.D11', ws.tmp2, ws.linsolve_ws_11)
+    ws.E11 .= view(e,1:nstable,1:nstable)
+    ws.tmp3 .= ws.E11'
+    linsolve_core!(ws.Z11', ws.tmp3, ws.linsolve_ws_11)
+    gemm!('T','T',1.0,ws.tmp2,ws.tmp3,0.0,ws.g1)
+end
+

src/PolynomialMatrixEquations.jl

+module PolynomialMatrixEquations
+
+include("CyclicReduction.jl")
+export CyclicReductionWs, cyclic_reduction!, cyclic_reduction_check
+
+include("GeneralizedSchurDecompositionSolver.jl")
+export GsSolverWs, gs_solver!
+
+
+
+end # module

test/testCycleReduction.jl

+using PolynomialMatrixEquations
+using LinearAlgebra
+using Test
+
+n = 3
+ws = CyclicReductionWs(n)
+a0 = [0.5 0 0; 0 0.5 0; 0 0 0];
+a1 = I(n) .+ zeros(n, n)
+a2 = [0 0 0; 0 0 0; 0 0 0.8]
+x = zeros(n,n)
+cyclic_reduction!(x,a0,a1,a2,ws,1e-8,50)
+@test ws.info == 0
+@test a0 + a1*x + a2*x*x ≈ zeros(n, n)
+display(x)
+
+cyclic_reduction_check(x,a0,a1,a2,1e-8)
+a0 = [0.5 0 0; 0 1.1 0; 0 0 0];
+a1 .= I(n) .+ zeros(n, n)
+a2 = [0 0 0; 0 0 0; 0 0 0.8]
+cyclic_reduction!(x,a0,a1,a2,ws,1e-8,50)
+@test ws.info == 1
+
+a0 = [0.5 0 0; 0 0.5 0; 0 0 0];
+a1 = I(n) .+ zeros(n, n)
+a2 = [0 0 0; 0 0 0; 0 0 1.2]
+cyclic_reduction!(x,a0,a1,a2,ws,1e-8,50)
+@test ws.info == 2

test/testGeneralizedSchurDecompositionSolver.jl

+using PolynomialMatrixEquations
+using LinearAlgebra
+using Random
+using Test
+
+Random.seed!(132)
+n = 3
+d_org = randn(n,n)
+e_org = randn(n,n)
+
+n1 = count(abs.(eigen(e, d).values) .< 1+1e-6)
+
+ws = GsSolverWs(d, e, n1)
+
+qz_criterium = 1.0 + 1e-6
+d = copy(d_org)
+e = copy(e_org)
+gs_solver!(ws, d , e, n1 , qz_criterium)
+#@test ws.info == 0
+@test d_org*[I(n1); ws.g2]*ws.g1 ≈ e_org*[I(n1); ws.g2]
+
+#=
+cyclic_reduction_check(x,a0,a1,a2,1e-8)
+a0 = [0.5 0 0; 0 1.1 0; 0 0 0];
+a1 .= I(n) .+ zeros(n, n)
+a2 = [0 0 0; 0 0 0; 0 0 0.8]
+cyclic_reduction!(x,a0,a1,a2,ws,1e-8,50)
+@test ws.info == 1
+
+a0 = [0.5 0 0; 0 0.5 0; 0 0 0];
+a1 = I(n) .+ zeros(n, n)
+a2 = [0 0 0; 0 0 0; 0 0 1.2]
+cyclic_reduction!(x,a0,a1,a2,ws,1e-8,50)
+@test ws.info == 2
+=#