search

RalphAS / PeriodicSchurDecompositions.jl

Julia package for periodic Schur decompositions of matrix products
Other
8 stars 0 forks source link

ordschur! fails on BigFloat data #7

Closed andreasvarga closed 1 year ago

andreasvarga commented 1 year ago

I tried to use pschur and ordschur! on extended precision data. The following code fails on BigFloat data:

using PeriodicSchurDecompositions
using JLD
hpd = load("test.jld","hpd");
hpds = [BigFloat.(hpd[i]) for i in 1:200];
PSFs = PeriodicSchurDecompositions.pschur(hpds,:L);
select = abs.(PSFs.values) .< 1;
PeriodicSchurDecompositions.ordschur!(PSFs,select);

ERROR: MethodError: no method matching rmul!(::SubArray{BigFloat, 2, Matrix{BigFloat}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}, ::LinearAlgebra.HessenbergQ{BigFloat, Matrix{BigFloat}, Vector{BigFloat}, false})
Closest candidates are:
  rmul!(::StridedMatrix{T} where T, ::LinearAlgebra.UpperTriangular{<:Any, <:LinearAlgebra.Adjoint}) at C:\Users\Andreas\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\triangular.jl:1108
  rmul!(::StridedMatrix{T} where T, ::LinearAlgebra.UnitUpperTriangular{<:Any, <:LinearAlgebra.Adjoint}) at C:\Users\Andreas\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\triangular.jl:1126
  rmul!(::StridedMatrix{T} where T, ::LinearAlgebra.LowerTriangular{<:Any, <:LinearAlgebra.Adjoint}) at C:\Users\Andreas\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\triangular.jl:1144
  ...
Stacktrace:
 [1] _swapadjqr!(T1::Matrix{BigFloat}, Ts::Vector{Matrix{BigFloat}}, Zs::Vector{Matrix{BigFloat}}, i1::Int64, p1::Int64, p2::Int64; sylcheck::Bool)
   @ PeriodicSchurDecompositions C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\sylswap.jl:138
 [2] _swapadjqr!
   @ C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\sylswap.jl:11 [inlined]
 [3] _swapschur!
   @ C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\rordschur.jl:253 [inlined]
 [4] _moveblock!(P::PeriodicSchur{BigFloat, Matrix{BigFloat}, Matrix{BigFloat}, Matrix{BigFloat}}, jsrc::Int64, jdest::Int64, wantZ::Bool, Q::Vector{Matrix{BigFloat}})
   @ PeriodicSchurDecompositions C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\rordschur.jl:181
 [5] ordschur!(P::PeriodicSchur{BigFloat, Matrix{BigFloat}, Matrix{BigFloat}, Matrix{BigFloat}}, select::BitVector; wantZ::Bool, Z::Nothing)
   @ PeriodicSchurDecompositions C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\rordschur.jl:98
 [6] ordschur!(P::PeriodicSchur{BigFloat, Matrix{BigFloat}, Matrix{BigFloat}, Matrix{BigFloat}}, select::BitVector)
   @ PeriodicSchurDecompositions C:\Users\Andreas\.julia\packages\PeriodicSchurDecompositions\AD2eu\src\rordschur.jl:6
 [7] top-level scope
   @ REPL[26]:1
RalphAS commented 1 year ago

The required method is implemented in the GenericLinearAlgebra package. I'm using a naive Sylvester solver, so the reordering takes a couple of minutes for BigFloat. If Float128 (from Quadmath.jl) is adequate, that is significantly faster.