GeneralizedSchurDecompositionSolver.jl 2.62 KB
Newer Older
Michel Juillard's avatar
Michel Juillard committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
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