Skip to content

Commit ccd298e

Browse files
committed
Merge pull request #12659 from mfasi/lapack_trexc
Add wrapper for LAPACK _trexc methods (reordering of the Schur form)
2 parents 9561900 + 6622f74 commit ccd298e

File tree

2 files changed

+73
-6
lines changed

2 files changed

+73
-6
lines changed

base/linalg/lapack.jl

Lines changed: 59 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -3858,10 +3858,37 @@ for (gees, gges, elty, relty) in
38583858
end
38593859
end
38603860
# Reorder Schur forms
3861-
for (trsen, tgsen, elty) in
3862-
((:dtrsen_, :dtgsen_, :Float64),
3863-
(:strsen_, :stgsen_, :Float32))
3861+
for (trexc, trsen, tgsen, elty) in
3862+
((:dtrexc_, :dtrsen_, :dtgsen_, :Float64),
3863+
(:strexc_, :strsen_, :stgsen_, :Float32))
38643864
@eval begin
3865+
trexc!(ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) = trexc!('V', ifst, ilst, T, Q)
3866+
function trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
3867+
# * .. Scalar Arguments ..
3868+
# CHARACTER COMPQ
3869+
# INTEGER IFST, ILST, INFO, LDQ, LDT, N
3870+
# * ..
3871+
# * .. Array Arguments ..
3872+
# DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
3873+
chkstride1(T, Q)
3874+
n = chksquare(T)
3875+
ldt = max(1, stride(T, 2))
3876+
ldq = max(1, stride(Q, 2))
3877+
work = Array($elty, n)
3878+
info = Array(BlasInt, 1)
3879+
3880+
ccall(($(blasfunc(trexc)), liblapack), Void,
3881+
(Ptr{UInt8}, Ptr{BlasInt},
3882+
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
3883+
Ptr{BlasInt}, Ptr{BlasInt},
3884+
Ptr{$elty}, Ptr{BlasInt}),
3885+
&compq, &n,
3886+
T, &ldt, Q, &ldq,
3887+
&ifst, &ilst,
3888+
work, info)
3889+
@lapackerror
3890+
T, Q
3891+
end
38653892
function trsen!(select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
38663893
# * .. Scalar Arguments ..
38673894
# CHARACTER COMPQ, JOB
@@ -3978,10 +4005,36 @@ for (trsen, tgsen, elty) in
39784005
end
39794006
end
39804007

3981-
for (trsen, tgsen, elty) in
3982-
((:ztrsen_, :ztgsen_, :Complex128),
3983-
(:ctrsen_, :ctgsen_, :Complex64))
4008+
for (trexc, trsen, tgsen, elty) in
4009+
((:ztrexc_, :ztrsen_, :ztgsen_, :Complex128),
4010+
(:ctrexc_, :ctrsen_, :ctgsen_, :Complex64))
39844011
@eval begin
4012+
trexc!(ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) = trexc!('V', ifst, ilst, T, Q)
4013+
function trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
4014+
# * .. Scalar Arguments ..
4015+
# CHARACTER COMPQ
4016+
# INTEGER IFST, ILST, INFO, LDQ, LDT, N
4017+
# * ..
4018+
# * .. Array Arguments ..
4019+
# DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
4020+
chkstride1(T, Q)
4021+
n = chksquare(T)
4022+
ldt = max(1, stride(T, 2))
4023+
ldq = max(1, stride(Q, 2))
4024+
info = Array(BlasInt, 1)
4025+
4026+
ccall(($(blasfunc(trexc)), liblapack), Void,
4027+
(Ptr{UInt8}, Ptr{BlasInt},
4028+
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
4029+
Ptr{BlasInt}, Ptr{BlasInt},
4030+
Ptr{BlasInt}),
4031+
&compq, &n,
4032+
T, &ldt, Q, &ldq,
4033+
&ifst, &ilst,
4034+
info)
4035+
@lapackerror
4036+
T, Q
4037+
end
39854038
function trsen!(select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
39864039
# * .. Scalar Arguments ..
39874040
# CHARACTER COMPQ, JOB

test/linalg/lapack.jl

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -455,6 +455,20 @@ for elty in (Float32, Float64, Complex64, Complex128)
455455
@test_throws DimensionMismatch LAPACK.laic1!(1,rand(elty,10),real(rand(elty)),rand(elty,11),rand(elty))
456456
end
457457

458+
#trexc
459+
for elty in (Float32, Float64, Complex64, Complex128)
460+
for c in ('V', 'N')
461+
A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4])
462+
T,Q,d = schur(A)
463+
Base.LinAlg.LAPACK.trexc!(c,LinAlg.BlasInt(1),LinAlg.BlasInt(2),T,Q)
464+
@test d[1] T[2,2]
465+
@test d[2] T[1,1]
466+
if (c == 'V')
467+
@test Q * T * Q' A
468+
end
469+
end
470+
end
471+
458472
# Test our own linear algebra functionaly against LAPACK
459473
for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
460474
for nn in (5,10,15)

0 commit comments

Comments
 (0)