SchurAlgo.jl 6.56 KB
Newer Older
1
module SchurAlgo
Michel Juillard's avatar
Michel Juillard committed
2
# general Schur decomposition with reordering
3
# adapted from ./base/linalg/lapack.jl
Michel Juillard's avatar
Michel Juillard committed
4
5
6
7

include("exceptions.jl")

import Base: USE_BLAS64, LAPACKException
8
9
10
import LinearAlgebra: BlasInt, BlasFloat, checksquare, chkstride1
import LinearAlgebra.BLAS: @blasfunc, libblas
import LinearAlgebra.LAPACK: liblapack, chklapackerror
Michel Juillard's avatar
Michel Juillard committed
11

12
13
export DgeesWS, dgees!, DggesWS, dgges!

Michel Juillard's avatar
Michel Juillard committed
14
15
16
const criterium = 1+1e-6

function mycompare{T}(wr_::Ptr{T}, wi_::Ptr{T})
17
18
19
  wr = unsafe_load(wr_)
  wi = unsafe_load(wi_)
  return convert(Cint, ((wr*wr + wi*wi) < criterium) ? 1 : 0)
Michel Juillard's avatar
Michel Juillard committed
20
21
22
end

function mycompare{T}(alphar_::Ptr{T}, alphai_::Ptr{T}, beta_::Ptr{T})
23
24
25
26
  alphar = unsafe_load(alphar_)
  alphai = unsafe_load(alphai_)
  beta = unsafe_load(beta_)
  return convert(Cint, ((alphar*alphar + alphai*alphai) < criterium*beta*beta) ? 1 : 0)
Michel Juillard's avatar
Michel Juillard committed
27
28
end

29
30
const mycompare_c = @cfunction(mycompare, Cint, (Ptr{Cdouble}, Ptr{Cdouble}))
const mycompare_g_c = @cfunction(mycompare, Cint, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}))
Michel Juillard's avatar
Michel Juillard committed
31

32
mutable struct DgeesWS
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
  jobvs::Ref{UInt8}
  sdim::Ref{BlasInt}
  wr::Vector{Float64}
  wi::Vector{Float64}
  ldvs::Ref{BlasInt}
  vs::Matrix{Float64}
  work::Vector{Float64}
  lwork::Ref{BlasInt}
  bwork::Vector{Int64}
  eigen_values::Vector{Complex{Float64}}
  info::Ref{BlasInt}

  function DgeesWS(jobvs::Ref{UInt8}, A::StridedMatrix{Float64}, sdim::Ref{BlasInt},
                    wr::Vector{Float64}, wi::Vector{Float64}, ldvs::Ref{BlasInt}, vs::Matrix{Float64},
                    work::Vector{Float64}, lwork::Ref{BlasInt}, bwork::Vector{Int64},
                    eigen_values::Vector{Complex{Float64}}, info::Ref{BlasInt})
      n = Ref{BlasInt}(size(A,1))
      RldA = Ref{BlasInt}(max(1,stride(A,2)))
      Rsort = Ref{UInt8}('N')
      ccall((@blasfunc(dgees_), liblapack), Cvoid,
            (Ref{UInt8}, Ref{UInt8}, Ptr{Cvoid},
             Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
             Ptr{BlasInt}, Ptr{Float64},
             Ptr{Float64}, Ptr{Float64}, Ref{BlasInt},
             Ptr{Float64}, Ref{BlasInt}, Ptr{Int64},
             Ref{BlasInt}),
            jobvs, Rsort, mycompare_c,
            n, A, RldA,
            sdim, wr, wi,
            vs, ldvs,
            work, lwork, bwork,
            info)
      chklapackerror(info[])
      lwork = Ref{BlasInt}(real(work[1]))
      work = Vector{Float64}(lwork[])
      new(jobvs, sdim, wr, wi, ldvs, vs, work, lwork, bwork, eigen_values, info)
  end
Michel Juillard's avatar
Michel Juillard committed
70
71
72
end

function DgeesWS(A::StridedMatrix{Float64})
73
74
75
76
77
78
79
80
81
82
83
84
85
86
  chkstride1(A)
  n, = checksquare(A)
  jobvs = Ref{UInt8}('V')
  sdim = Ref{BlasInt}(0)
  wr = Vector{Float64}(n)
  wi = Vector{Float64}(n)
  ldvs = Ref{BlasInt}(jobvs[] == UInt32('V') ? n : 1)
  vs = Matrix{Float64}(ldvs[], n)
  work = Vector{Float64}(1)
  lwork = Ref{BlasInt}(-1)
  bwork = Vector{Int64}(n)
  eigen_values = Vector{Complex{Float64}}(n)
  info = Ref{BlasInt}(0)
  DgeesWS(jobvs, A, sdim, wr, wi, ldvs, vs, work, lwork, bwork, eigen_values, info)
Michel Juillard's avatar
Michel Juillard committed
87
end
88
89

function DgeesWS(n::Int64)
90
91
  A = zeros(n,n)
  DgeesWS(A)
92
93
end

Michel Juillard's avatar
Michel Juillard committed
94
function dgees!(ws::DgeesWS,A::StridedMatrix{Float64})
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
  n = Ref{BlasInt}(size(A,1))
  RldA = Ref{BlasInt}(max(1,stride(A,2)))
  sort = Ref{UInt8}('S')
  ccall((@blasfunc(dgees_), liblapack), Cvoid,
        (Ref{UInt8}, Ref{UInt8}, Ptr{Cvoid},
         Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
         Ref{BlasInt}, Ptr{Float64}, Ptr{Float64},
         Ptr{Float64}, Ref{BlasInt},
         Ptr{Float64}, Ref{BlasInt}, Ptr{Int64},
         Ref{BlasInt}),
        ws.jobvs, sort, mycompare_c,
        n, A, RldA,
        ws.sdim, ws.wr, ws.wi,
        ws.vs, ws.ldvs,
        ws.work, ws.lwork, ws.bwork,
        ws.info)
  ws.eigen_values = complex.(ws.wr, ws.wi)
  chklapackerror(ws.info[])
Michel Juillard's avatar
Michel Juillard committed
113
114
end

115
mutable struct DggesWS
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
  alphar::Vector{Float64}
  alphai::Vector{Float64}
  beta::Vector{Float64}
  lwork::BlasInt
  work::Vector{Float64}
  bwork::Vector{Int64}
  sdim::BlasInt

  function DggesWS(jobvsl::Ref{UInt8}, jobvsr::Ref{UInt8}, sort::Ref{UInt8}, A::StridedMatrix{Float64}, B::StridedMatrix{Float64})
      chkstride1(A, B)
      n, m = checksquare(A, B)
      if n != m
          throw(DimensionMismatch("Dimensions of A, ($n,$n), and B, ($m,$m), must match"))
      end
      n = BlasInt(size(A,1))
      alphar = Vector{Float64}(n)
      alphai = Vector{Float64}(n)
      beta = Vector{Float64}(n)
      bwork = Vector{Int64}(n)
      ldvsl = BlasInt(1)
      ldvsr = BlasInt(1)
      sdim = BlasInt(0)
      lwork = BlasInt(-1)
      work = Vector{Float64}(1)
      sdim = BlasInt(0)
      info = BlasInt(0)
      ccall((@blasfunc(dgges_), liblapack), Cvoid,
            (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ptr{Cvoid},
             Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64},
             Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ptr{Float64},
             Ptr{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64},
             Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Int64},
             Ref{BlasInt}),
            jobvsl, jobvsr, sort, mycompare_c,
            n, A, max(1,stride(A, 2)), B,
            max(1,stride(B, 2)), sdim, alphar, alphai,
            beta, C_NULL, ldvsl, C_NULL,
            ldvsr, work, lwork, bwork,
            info)
      chklapackerror(info)
      println(real(work[1]))
      lwork = BlasInt(real(work[1]))
      work = Vector{Float64}(lwork)
      new(alphar,alphai,beta,lwork,work,bwork,sdim)
  end
Michel Juillard's avatar
Michel Juillard committed
161
162
end

Michel Juillard's avatar
Michel Juillard committed
163
164
165
166
function DggesWS(A::StridedMatrix{Float64}, B::StridedMatrix{Float64})
    DggesWS(Ref{UInt8}('N'), Ref{UInt8}('N'), Ref{UInt8}('N'), A, B)
end

Michel Juillard's avatar
Michel Juillard committed
167
168
169
170
171
172
173
174
175
function dgges!(jobvsl::Char, jobvsr::Char, A::StridedMatrix{Float64}, B::StridedMatrix{Float64},
                vsl::Matrix{Float64}, vsr::Matrix{Float64}, eigval::Array{Complex64,1},
                ws::DggesWS)
    n = size(A,1)
    ldvsl = jobvsl == 'V' ? n : 1
    ldvsr = jobvsr == 'V' ? n : 1
    sort = 'S'
    sdim = Ref{BlasInt}(0)
    info = Ref{BlasInt}(0)
176
177
    ccall((@blasfunc(dgges_), liblapack), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ptr{Cvoid},
Michel Juillard's avatar
Michel Juillard committed
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
           Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64},
           Ref{BlasInt}, Ptr{BlasInt}, Ptr{Float64}, Ptr{Float64},
           Ptr{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64},
           Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Int64},
           Ref{Int64}),
          jobvsl, jobvsr, sort, mycompare_c,
          n, A, max(1,stride(A, 2)), B,
          max(1,stride(B, 2)), sdim, ws.alphar, ws.alphai,
          ws.beta, vsl, ldvsl, vsr,
          ldvsr, ws.work, ws.lwork, ws.bwork,
          info)
    ws.sdim = sdim[]
    if info[]
        throw(DggesException(info[]))
    end
    for i in 1:n
        eigval[i] = complex(ws.alphar[i],ws.alphai[i])/ws.beta[i]
    end
end

end