Skip to content

Commit e0624d4

Browse files
MichelJuillardMichelJuillard
committed
works with both Julia v1.6 and v1.7
1 parent b7cbb38 commit e0624d4

7 files changed

Lines changed: 67 additions & 104 deletions

File tree

‎src/FastLapackInterface.jl‎

Lines changed: 20 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,26 @@
11
module FastLapackInterface
22

3+
import Base.strides
4+
5+
using LinearAlgebra
6+
using LinearAlgebra: BlasInt, BlasFloat, checksquare, chkstride1
7+
using LinearAlgebra.BLAS: @blasfunc
8+
using LinearAlgebra.LAPACK: chklapackerror
9+
10+
11+
@static if VERSION < v"1.7"
12+
using LinearAlgebra.LAPACK: liblapack
13+
else
14+
const liblapack = "libblastrampoline"
15+
end
16+
317
include("LinSolveAlgo.jl")
18+
export LinSolveWs, linsolve_core!, linsolve_core_no_lu!, lu!
419
include("QrAlgo.jl")
20+
export QrWs, QrpWs, geqrf_core!, geqp3!, ormqr_core!
521
include("SchurAlgo.jl")
22+
export DgeesWs, dgees!, DggesWs, dgges!
23+
24+
end #module
25+
#import LinearAlgebra: USE_BLAS64, LAPACKException
626

7-
end # module

‎src/LinSolveAlgo.jl‎

Lines changed: 2 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,3 @@
1-
module LinSolveAlgo
2-
3-
import Base.strides
4-
5-
const libblastrampoline = "libblastrampoline"
6-
7-
using LinearAlgebra
8-
import LinearAlgebra.BlasInt
9-
import LinearAlgebra.BLAS.@blasfunc
10-
import LinearAlgebra.LAPACK: chklapackerror
11-
12-
export LinSolveWs, linsolve_core!, linsolve_core_no_lu!, lu!
13-
141
struct LinSolveWs{T<:Number,U<:Integer}
152
lu::Vector{T}
163
ipiv::Vector{BlasInt}
@@ -42,7 +29,7 @@ for (getrf, getrs, elty) in (
4229
lda = Ref{BlasInt}(max(1, mm))
4330
info = Ref{BlasInt}(0)
4431
ccall(
45-
(@blasfunc($getrf), libblastrampoline),
32+
(@blasfunc($getrf), liblapack),
4633
Cvoid,
4734
(
4835
Ref{BlasInt},
@@ -88,7 +75,7 @@ for (getrf, getrs, elty) in (
8875
ldb = Ref{BlasInt}(max(1, stride(b, 2)))
8976
info = Ref{BlasInt}(0)
9077
ccall(
91-
(@blasfunc($getrs), libblastrampoline),
78+
(@blasfunc($getrs), liblapack),
9279
Cvoid,
9380
(
9481
Ref{UInt8},
@@ -129,5 +116,3 @@ end
129116
# for its parent as the transposition is done at the
130117
# solution stage
131118
lu!(a::Adjoint, ws::LinSolveWs) = lu!(a.parent, ws)
132-
133-
end

‎src/QrAlgo.jl‎

Lines changed: 8 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,3 @@
1-
module QrAlgo
2-
3-
const libblastrampoline = "libblastrampoline"
4-
5-
using LinearAlgebra
6-
import LinearAlgebra: BlasInt
7-
import LinearAlgebra.BLAS: @blasfunc
8-
import LinearAlgebra.LAPACK: chklapackerror
9-
10-
export QrWs, QrpWs, geqrf_core!, geqp3!, ormqr_core!
11-
121
abstract type QR end
132

143
struct QrWs{T<:Number} <: QR
@@ -37,7 +26,7 @@ for (geqrf, ormqr, elty) in (
3726
lwork = Ref{BlasInt}(-1)
3827
info = Ref{BlasInt}(0)
3928
ccall(
40-
(@blasfunc($geqrf), libblastrampoline),
29+
(@blasfunc($geqrf), liblapack),
4130
Nothing,
4231
(
4332
Ref{BlasInt},
@@ -70,7 +59,7 @@ for (geqrf, ormqr, elty) in (
7059
n = Ref{BlasInt}(nn)
7160
RldA = Ref{BlasInt}(max(1, stride(A, 2)))
7261
ccall(
73-
(@blasfunc($geqrf), libblastrampoline),
62+
(@blasfunc($geqrf), liblapack),
7463
Nothing,
7564
(
7665
Ref{BlasInt},
@@ -116,7 +105,7 @@ for (geqrf, ormqr, elty) in ((:dgeqrf_, :dormqr_, :Float64), (:sgeqrf_, :sormqr_
116105
RldA = Ref{BlasInt}(max(1, stride(A, 2)))
117106
RldC = Ref{BlasInt}(max(1, stride(C, 2)))
118107
ccall(
119-
(@blasfunc($ormqr), libblastrampoline),
108+
(@blasfunc($ormqr), liblapack),
120109
Nothing,
121110
(
122111
Ref{UInt8},
@@ -167,7 +156,7 @@ for (geqrf, ormqr, elty) in ((:dgeqrf_, :dormqr_, :Float64), (:sgeqrf_, :sormqr_
167156
RldA = Ref{BlasInt}(max(1, stride(A.parent, 2)))
168157
RldC = Ref{BlasInt}(max(1, stride(C, 2)))
169158
ccall(
170-
(@blasfunc($ormqr), libblastrampoline),
159+
(@blasfunc($ormqr), liblapack),
171160
Nothing,
172161
(
173162
Ref{UInt8},
@@ -221,7 +210,7 @@ for (geqrf, ormqr, elty) in
221210
k = Ref{BlasInt}(length(ws.tau))
222211
RldA = Ref{BlasInt}(max(1,stride(A,2)))
223212
RldC = Ref{BlasInt}(max(1,stride(C,2)))
224-
ccall((@blasfunc($ormqr), libblastrampoline), Nothing,
213+
ccall((@blasfunc($ormqr), liblapack), Nothing,
225214
(Ref{UInt8},Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ref{BlasInt},Ptr{$elty},Ref{BlasInt},
226215
Ptr{$elty},Ptr{$elty},Ref{BlasInt},Ptr{$elty},Ref{BlasInt},Ref{BlasInt}),
227216
side, 'N', m, n, k, A, RldA, ws.tau, C, RldC, ws.work, ws.lwork, ws.info)
@@ -244,7 +233,7 @@ for (geqrf, ormqr, elty) in
244233
k = Ref{BlasInt}(length(ws.tau))
245234
RldA = Ref{BlasInt}(max(1,stride(A.parent,2)))
246235
RldC = Ref{BlasInt}(max(1,stride(C,2)))
247-
ccall((@blasfunc($ormqr), libblastrampoline), Nothing,
236+
ccall((@blasfunc($ormqr), liblapack), Nothing,
248237
(Ref{UInt8},Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ref{BlasInt},Ptr{$elty},Ref{BlasInt},
249238
Ptr{$elty},Ptr{$elty},Ref{BlasInt},Ptr{$elty},Ref{BlasInt},Ref{BlasInt}),
250239
side, $transchar, m, n, k, A.parent, RldA,
@@ -282,7 +271,7 @@ for (geqp3, elty) in (
282271
lwork = BlasInt(-1)
283272
info = Ref{BlasInt}()
284273
ccall(
285-
(@blasfunc($geqp3), libblastrampoline),
274+
(@blasfunc($geqp3), liblapack),
286275
Nothing,
287276
(
288277
Ref{BlasInt},
@@ -315,7 +304,7 @@ for (geqp3, elty) in (
315304
m, n = size(A)
316305
RldA = BlasInt(max(1, stride(A, 2)))
317306
ccall(
318-
(@blasfunc($geqp3), libblastrampoline),
307+
(@blasfunc($geqp3), liblapack),
319308
Nothing,
320309
(
321310
Ref{BlasInt},
@@ -347,4 +336,3 @@ for (geqp3, elty) in (
347336
end
348337
end
349338

350-
end

‎src/SchurAlgo.jl‎

Lines changed: 5 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1,29 +1,11 @@
1-
module SchurAlgo
21
# general Schur decomposition with reordering
32
# adataped from ./base/linalg/lapack.jl
43

54
include("exceptions.jl")
65

7-
const libblastrampoline = "libblastrampoline"
8-
9-
import LinearAlgebra: USE_BLAS64, LAPACKException
10-
import LinearAlgebra: BlasInt, BlasFloat, checksquare, chkstride1
11-
import LinearAlgebra.BLAS: @blasfunc, libblas
12-
import LinearAlgebra.LAPACK: chklapackerror
13-
import Base: has_offset_axes
14-
15-
export DgeesWs, dgees!, DggesWs, dgges!
166

177
const criterium = 1 + 1e-6
188

19-
#=
20-
function mycompare(wr_, wi_)::Cint
21-
wr = unsafe_load(wr_)
22-
wi = unsafe_load(wi_)
23-
return convert(Cint, ((wr*wr + wi*wi) < criterium) ? 1 : 0)
24-
end
25-
=#
26-
279
function mycompare(alphar_::Ptr{T}, alphai_::Ptr{T}, beta_::Ptr{T})::Cint where {T}
2810
alphar = unsafe_load(alphar_)
2911
alphai = unsafe_load(alphai_)
@@ -64,9 +46,8 @@ struct DgeesWs
6446
n = Ref{BlasInt}(size(A, 1))
6547
RldA = Ref{BlasInt}(max(1, stride(A, 2)))
6648
Rsort = Ref{UInt8}('N')
67-
# mycompare_c = @cfunction(mycompare, Cint, (Ptr{Cdouble}, Ptr{Cdouble}))
6849
ccall(
69-
(@blasfunc(dgees_), libblastrampoline),
50+
(@blasfunc(dgees_), liblapack),
7051
Nothing,
7152
(
7253
Ref{UInt8},
@@ -138,7 +119,7 @@ function dgees!(ws::DgeesWs, A::StridedMatrix{Float64})
138119
RldA = Ref{BlasInt}(max(1, stride(A, 2)))
139120
myfunc::Function = make_select_function(>=, 1.0)
140121
ccall(
141-
(@blasfunc(dgees_), libblastrampoline),
122+
(@blasfunc(dgees_), liblapack),
142123
Cvoid,
143124
(
144125
Ref{UInt8},
@@ -192,7 +173,7 @@ function dgees!(ws::DgeesWs, A::StridedMatrix{Float64}, op, crit)
192173
myfunc::Function = make_select_function(op, crit)
193174
mycompare_c = @cfunction($myfunc, Cint, (Ptr{Cdouble}, Ptr{Cdouble}))
194175
ccall(
195-
(@blasfunc(dgees_), libblastrampoline),
176+
(@blasfunc(dgees_), liblapack),
196177
Cvoid,
197178
(
198179
Ref{UInt8},
@@ -267,7 +248,7 @@ mutable struct DggesWs
267248
mycompare_g_c =
268249
@cfunction(mycompare, Cint, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}))
269250
ccall(
270-
(@blasfunc(dgges_), libblastrampoline),
251+
(@blasfunc(dgges_), liblapack),
271252
Nothing,
272253
(
273254
Ref{UInt8},
@@ -343,7 +324,7 @@ function dgges!(
343324
info = Ref{BlasInt}(0)
344325
mycompare_g_c = @cfunction(mycompare, Cint, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}))
345326
ccall(
346-
(@blasfunc(dgges_), libblastrampoline),
327+
(@blasfunc(dgges_), liblapack),
347328
Nothing,
348329
(
349330
Ref{UInt8},
@@ -399,4 +380,3 @@ function dgges!(
399380
end
400381
end
401382

402-
end

‎test/LinSolveAlgo_test.jl‎

Lines changed: 24 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,3 @@
1-
include("../src/LinSolveAlgo.jl")
2-
3-
using LinearAlgebra
4-
using Test
5-
61
n = 3
72
m = 2
83

@@ -15,119 +10,119 @@ for elty in (Float32, Float64, ComplexF32, ComplexF64)
1510
A = copy(A0)
1611

1712
# Full matrix
18-
linws = LinSolveAlgo.LinSolveWs{elty,Int64}(n)
13+
linws = FastLapackInterface.LinSolveWs{elty,Int64}(n)
1914

20-
LinSolveAlgo.lu!(A, linws)
15+
FastLapackInterface.lu!(A, linws)
2116
@test A == A0
2217
F = lu(A)
2318
@test UpperTriangular(reshape(linws.lu, n, n)) ≈ F.U
2419

25-
LinSolveAlgo.lu!(A', linws)
20+
FastLapackInterface.lu!(A', linws)
2621
@test UpperTriangular(reshape(linws.lu, n, n)) ≈ F.U
2722

2823
B = copy(B0)
29-
LinSolveAlgo.linsolve_core!(A, B, linws)
24+
FastLapackInterface.linsolve_core!(A, B, linws)
3025
@test A == A0
3126
@test B ≈ A \ B0
3227

3328
copy!(B, B1)
34-
LinSolveAlgo.linsolve_core_no_lu!(A, B, linws)
29+
FastLapackInterface.linsolve_core_no_lu!(A, B, linws)
3530
@test A == A0
3631
@test B ≈ A \ B1
3732

3833
copy!(B, B0)
39-
LinSolveAlgo.linsolve_core!(A', B, linws)
34+
FastLapackInterface.linsolve_core!(A', B, linws)
4035
@test B ≈ A' \ B0
4136

4237
copy!(B, B1)
43-
LinSolveAlgo.linsolve_core_no_lu!(A', B, linws)
38+
FastLapackInterface.linsolve_core_no_lu!(A', B, linws)
4439
@test A == A0
4540
@test B ≈ A' \ B1
4641

4742
# view of a matrix in upper left corner
48-
linws1 = LinSolveAlgo.LinSolveWs{elty,Int64}(n - 1)
43+
linws1 = FastLapackInterface.LinSolveWs{elty,Int64}(n - 1)
4944
C = view(A, 1:n-1, 1:n-1)
5045
D = view(B, 1:n-1, 1:m-1)
5146
D0 = copy(D)
5247
D1 = copy(D)
53-
LinSolveAlgo.linsolve_core!(C, D, linws1)
48+
FastLapackInterface.linsolve_core!(C, D, linws1)
5449
@test C == view(A, 1:n-1, 1:n-1)
5550
@test D ≈ C \ D0
5651

5752
D = copy(D1)
58-
LinSolveAlgo.linsolve_core_no_lu!(C, D, linws1)
53+
FastLapackInterface.linsolve_core_no_lu!(C, D, linws1)
5954
@test C == view(A, 1:n-1, 1:n-1)
6055
@test D ≈ C \ D1
6156

6257
D = copy(D0)
63-
LinSolveAlgo.linsolve_core!(C', D, linws1)
58+
FastLapackInterface.linsolve_core!(C', D, linws1)
6459
@test C == view(A, 1:n-1, 1:n-1)
6560
@test D ≈ C' \ D0
6661

6762
D = copy(D1)
68-
LinSolveAlgo.linsolve_core_no_lu!(C', D, linws1)
63+
FastLapackInterface.linsolve_core_no_lu!(C', D, linws1)
6964
@test C == view(A, 1:n-1, 1:n-1)
7065
@test D ≈ C' \ D1
7166

7267
# view of a matrix in lower left corner
73-
linws1 = LinSolveAlgo.LinSolveWs{elty,Int64}(n - 1)
68+
linws1 = FastLapackInterface.LinSolveWs{elty,Int64}(n - 1)
7469
C = view(A, 2:n, 1:n-1)
7570
C1 = copy(C)
76-
LinSolveAlgo.lu!(C1, linws1)
71+
FastLapackInterface.lu!(C1, linws1)
7772
F = LinearAlgebra.lu!(C)
7873
@test triu(reshape(linws1.lu, n - 1, n - 1)) ≈ F.U
7974
@test tril(reshape(linws1.lu, n - 1, n - 1), -1) ≈ tril(F.L, -1)
8075

8176
D0 = view(B0, 2:n, 1:m-1)
8277
D1 = view(B1, 2:n, 1:m-1)
8378
D = copy(D0)
84-
LinSolveAlgo.linsolve_core!(C, D, linws1)
79+
FastLapackInterface.linsolve_core!(C, D, linws1)
8580
@test C == view(A, 2:n, 1:n-1)
8681
@test D ≈ C \ D0
8782

8883
D = copy(D1)
89-
LinSolveAlgo.linsolve_core_no_lu!(C, D, linws1)
84+
FastLapackInterface.linsolve_core_no_lu!(C, D, linws1)
9085
@test C == view(A, 2:n, 1:n-1)
9186
@test D ≈ C \ D1
9287

9388
D = copy(D0)
94-
LinSolveAlgo.linsolve_core!(C', D, linws1)
89+
FastLapackInterface.linsolve_core!(C', D, linws1)
9590
@test C == view(A, 2:n, 1:n-1)
9691
@test D ≈ C' \ D0
9792

9893
D = copy(D1)
99-
LinSolveAlgo.linsolve_core_no_lu!(C', D, linws1)
94+
FastLapackInterface.linsolve_core_no_lu!(C', D, linws1)
10095
@test C == view(A, 2:n, 1:n-1)
10196
@test D ≈ C' \ D1
10297

10398
# using too big a work space
104-
linws1 = LinSolveAlgo.LinSolveWs{elty,Int64}(n)
99+
linws1 = FastLapackInterface.LinSolveWs{elty,Int64}(n)
105100
C = view(A, 2:n, 1:n-1)
106101
C1 = copy(C)
107-
LinSolveAlgo.lu!(C1, linws1)
102+
FastLapackInterface.lu!(C1, linws1)
108103
F = LinearAlgebra.lu!(C)
109104
@test triu(reshape(linws1.lu[1:(n-1)^2], n - 1, n - 1)) ≈ F.U
110105
@test tril(reshape(linws1.lu[1:(n-1)^2], n - 1, n - 1), -1) ≈ tril(F.L, -1)
111106

112107
D0 = view(B0, 2:n, 1:m-1)
113108
D1 = view(B1, 2:n, 1:m-1)
114109
D = copy(D0)
115-
LinSolveAlgo.linsolve_core!(C, D, linws1)
110+
FastLapackInterface.linsolve_core!(C, D, linws1)
116111
@test C == view(A, 2:n, 1:n-1)
117112
@test D ≈ C \ D0
118113

119114
D = copy(D1)
120-
LinSolveAlgo.linsolve_core_no_lu!(C, D, linws1)
115+
FastLapackInterface.linsolve_core_no_lu!(C, D, linws1)
121116
@test C == view(A, 2:n, 1:n-1)
122117
@test D ≈ C \ D1
123118

124119
D = copy(D0)
125-
LinSolveAlgo.linsolve_core!(C', D, linws1)
120+
FastLapackInterface.linsolve_core!(C', D, linws1)
126121
@test C == view(A, 2:n, 1:n-1)
127122
@test D ≈ C' \ D0
128123

129124
D = copy(D1)
130-
LinSolveAlgo.linsolve_core_no_lu!(C', D, linws1)
125+
FastLapackInterface.linsolve_core_no_lu!(C', D, linws1)
131126
@test C == view(A, 2:n, 1:n-1)
132127
@test D ≈ C' \ D1
133128

0 commit comments

Comments
 (0)