3
3
OpenBLASLUFactorization()
4
4
```
5
5
6
- A direct wrapper over OpenBLAS's LU factorization (`getrf!` and `getrs!`).
7
- This solver makes direct calls to OpenBLAS_jll without going through Julia's
6
+ A direct wrapper over OpenBLAS's LU factorization (`getrf!` and `getrs!`).
7
+ This solver makes direct calls to OpenBLAS_jll without going through Julia's
8
8
libblastrampoline, which can provide performance benefits in certain configurations.
9
9
10
10
## Performance Characteristics
11
11
12
- - Pre-allocates workspace to avoid allocations during solving
13
- - Makes direct `ccall`s to OpenBLAS routines
14
- - Can be faster than `LUFactorization` when OpenBLAS is well-optimized for the hardware
15
- - Supports `Float32`, `Float64`, `ComplexF32`, and `ComplexF64` element types
12
+ - Pre-allocates workspace to avoid allocations during solving
13
+ - Makes direct `ccall`s to OpenBLAS routines
14
+ - Can be faster than `LUFactorization` when OpenBLAS is well-optimized for the hardware
15
+ - Supports `Float32`, `Float64`, `ComplexF32`, and `ComplexF64` element types
16
16
17
17
## When to Use
18
18
19
- - When you want to ensure OpenBLAS is used regardless of the system BLAS configuration
20
- - When benchmarking shows better performance than `LUFactorization` on your specific hardware
21
- - When you need consistent behavior across different systems (always uses OpenBLAS)
19
+ - When you want to ensure OpenBLAS is used regardless of the system BLAS configuration
20
+ - When benchmarking shows better performance than `LUFactorization` on your specific hardware
21
+ - When you need consistent behavior across different systems (always uses OpenBLAS)
22
22
23
23
## Example
24
24
@@ -36,96 +36,96 @@ struct OpenBLASLUFactorization <: AbstractFactorization end
36
36
module OpenBLASLU
37
37
38
38
using LinearAlgebra
39
- using LinearAlgebra: BlasInt, LU, require_one_based_indexing, checksquare
40
- using LinearAlgebra . LAPACK : chkfinite, chkstride1, @blasfunc , chkargsok, chktrans, chklapackerror
39
+ using LinearAlgebra. LAPACK : chkfinite, chkstride1, @blasfunc , chkargsok, chktrans,
40
+ chklapackerror
41
41
using OpenBLAS_jll
42
42
43
43
function getrf! (A:: AbstractMatrix{<:ComplexF64} ;
44
- ipiv = similar (A, BlasInt, min (size (A, 1 ), size (A, 2 ))),
45
- info = Ref {BlasInt} (),
44
+ ipiv = similar (A, LinearAlgebra . BlasInt, min (size (A, 1 ), size (A, 2 ))),
45
+ info = Ref {LinearAlgebra. BlasInt} (),
46
46
check = false )
47
- require_one_based_indexing (A)
47
+ LinearAlgebra . require_one_based_indexing (A)
48
48
check && chkfinite (A)
49
49
chkstride1 (A)
50
50
m, n = size (A)
51
51
lda = max (1 , stride (A, 2 ))
52
52
if isempty (ipiv)
53
- ipiv = similar (A, BlasInt, min (size (A, 1 ), size (A, 2 )))
53
+ ipiv = similar (A, LinearAlgebra . BlasInt, min (size (A, 1 ), size (A, 2 )))
54
54
end
55
55
ccall ((@blasfunc (zgetrf_), OpenBLAS_jll. libopenblas), Cvoid,
56
- (Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64},
57
- Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
56
+ (Ref{LinearAlgebra . BlasInt}, Ref{LinearAlgebra . BlasInt}, Ptr{ComplexF64},
57
+ Ref{LinearAlgebra . BlasInt}, Ptr{LinearAlgebra . BlasInt}, Ptr{LinearAlgebra . BlasInt}),
58
58
m, n, A, lda, ipiv, info)
59
59
chkargsok (info[])
60
60
A, ipiv, info[], info # Error code is stored in LU factorization type
61
61
end
62
62
63
63
function getrf! (A:: AbstractMatrix{<:ComplexF32} ;
64
- ipiv = similar (A, BlasInt, min (size (A, 1 ), size (A, 2 ))),
65
- info = Ref {BlasInt} (),
64
+ ipiv = similar (A, LinearAlgebra . BlasInt, min (size (A, 1 ), size (A, 2 ))),
65
+ info = Ref {LinearAlgebra. BlasInt} (),
66
66
check = false )
67
- require_one_based_indexing (A)
67
+ LinearAlgebra . require_one_based_indexing (A)
68
68
check && chkfinite (A)
69
69
chkstride1 (A)
70
70
m, n = size (A)
71
71
lda = max (1 , stride (A, 2 ))
72
72
if isempty (ipiv)
73
- ipiv = similar (A, BlasInt, min (size (A, 1 ), size (A, 2 )))
73
+ ipiv = similar (A, LinearAlgebra . BlasInt, min (size (A, 1 ), size (A, 2 )))
74
74
end
75
75
ccall ((@blasfunc (cgetrf_), OpenBLAS_jll. libopenblas), Cvoid,
76
- (Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF32},
77
- Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
76
+ (Ref{LinearAlgebra . BlasInt}, Ref{LinearAlgebra . BlasInt}, Ptr{ComplexF32},
77
+ Ref{LinearAlgebra . BlasInt}, Ptr{LinearAlgebra . BlasInt}, Ptr{LinearAlgebra . BlasInt}),
78
78
m, n, A, lda, ipiv, info)
79
79
chkargsok (info[])
80
80
A, ipiv, info[], info # Error code is stored in LU factorization type
81
81
end
82
82
83
83
function getrf! (A:: AbstractMatrix{<:Float64} ;
84
- ipiv = similar (A, BlasInt, min (size (A, 1 ), size (A, 2 ))),
85
- info = Ref {BlasInt} (),
84
+ ipiv = similar (A, LinearAlgebra . BlasInt, min (size (A, 1 ), size (A, 2 ))),
85
+ info = Ref {LinearAlgebra. BlasInt} (),
86
86
check = false )
87
- require_one_based_indexing (A)
87
+ LinearAlgebra . require_one_based_indexing (A)
88
88
check && chkfinite (A)
89
89
chkstride1 (A)
90
90
m, n = size (A)
91
91
lda = max (1 , stride (A, 2 ))
92
92
if isempty (ipiv)
93
- ipiv = similar (A, BlasInt, min (size (A, 1 ), size (A, 2 )))
93
+ ipiv = similar (A, LinearAlgebra . BlasInt, min (size (A, 1 ), size (A, 2 )))
94
94
end
95
95
ccall ((@blasfunc (dgetrf_), OpenBLAS_jll. libopenblas), Cvoid,
96
- (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64},
97
- Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
96
+ (Ref{LinearAlgebra . BlasInt}, Ref{LinearAlgebra . BlasInt}, Ptr{Float64},
97
+ Ref{LinearAlgebra . BlasInt}, Ptr{LinearAlgebra . BlasInt}, Ptr{LinearAlgebra . BlasInt}),
98
98
m, n, A, lda, ipiv, info)
99
99
chkargsok (info[])
100
100
A, ipiv, info[], info # Error code is stored in LU factorization type
101
101
end
102
102
103
103
function getrf! (A:: AbstractMatrix{<:Float32} ;
104
- ipiv = similar (A, BlasInt, min (size (A, 1 ), size (A, 2 ))),
105
- info = Ref {BlasInt} (),
104
+ ipiv = similar (A, LinearAlgebra . BlasInt, min (size (A, 1 ), size (A, 2 ))),
105
+ info = Ref {LinearAlgebra. BlasInt} (),
106
106
check = false )
107
- require_one_based_indexing (A)
107
+ LinearAlgebra . require_one_based_indexing (A)
108
108
check && chkfinite (A)
109
109
chkstride1 (A)
110
110
m, n = size (A)
111
111
lda = max (1 , stride (A, 2 ))
112
112
if isempty (ipiv)
113
- ipiv = similar (A, BlasInt, min (size (A, 1 ), size (A, 2 )))
113
+ ipiv = similar (A, LinearAlgebra . BlasInt, min (size (A, 1 ), size (A, 2 )))
114
114
end
115
115
ccall ((@blasfunc (sgetrf_), OpenBLAS_jll. libopenblas), Cvoid,
116
- (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float32},
117
- Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
116
+ (Ref{LinearAlgebra . BlasInt}, Ref{LinearAlgebra . BlasInt}, Ptr{Float32},
117
+ Ref{LinearAlgebra . BlasInt}, Ptr{LinearAlgebra . BlasInt}, Ptr{LinearAlgebra . BlasInt}),
118
118
m, n, A, lda, ipiv, info)
119
119
chkargsok (info[])
120
120
A, ipiv, info[], info # Error code is stored in LU factorization type
121
121
end
122
122
123
123
function getrs! (trans:: AbstractChar ,
124
124
A:: AbstractMatrix{<:ComplexF64} ,
125
- ipiv:: AbstractVector{BlasInt} ,
125
+ ipiv:: AbstractVector{LinearAlgebra. BlasInt} ,
126
126
B:: AbstractVecOrMat{<:ComplexF64} ;
127
- info = Ref {BlasInt} ())
128
- require_one_based_indexing (A, ipiv, B)
127
+ info = Ref {LinearAlgebra. BlasInt} ())
128
+ LinearAlgebra . require_one_based_indexing (A, ipiv, B)
129
129
LinearAlgebra. LAPACK. chktrans (trans)
130
130
chkstride1 (A, B, ipiv)
131
131
n = LinearAlgebra. checksquare (A)
@@ -137,20 +137,22 @@ function getrs!(trans::AbstractChar,
137
137
end
138
138
nrhs = size (B, 2 )
139
139
ccall ((@blasfunc (zgetrs_), OpenBLAS_jll. libopenblas), Cvoid,
140
- (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt},
141
- Ptr{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, Ptr{BlasInt}, Clong),
140
+ (Ref{UInt8}, Ref{LinearAlgebra. BlasInt}, Ref{LinearAlgebra. BlasInt},
141
+ Ptr{ComplexF64}, Ref{LinearAlgebra. BlasInt},
142
+ Ptr{LinearAlgebra. BlasInt}, Ptr{ComplexF64}, Ref{LinearAlgebra. BlasInt},
143
+ Ptr{LinearAlgebra. BlasInt}, Clong),
142
144
trans, n, size (B, 2 ), A, max (1 , stride (A, 2 )), ipiv, B, max (1 , stride (B, 2 )), info,
143
145
1 )
144
- LinearAlgebra. LAPACK. chklapackerror (BlasInt (info[]))
146
+ LinearAlgebra. LAPACK. chklapackerror (LinearAlgebra . BlasInt (info[]))
145
147
B
146
148
end
147
149
148
150
function getrs! (trans:: AbstractChar ,
149
151
A:: AbstractMatrix{<:ComplexF32} ,
150
- ipiv:: AbstractVector{BlasInt} ,
152
+ ipiv:: AbstractVector{LinearAlgebra. BlasInt} ,
151
153
B:: AbstractVecOrMat{<:ComplexF32} ;
152
- info = Ref {BlasInt} ())
153
- require_one_based_indexing (A, ipiv, B)
154
+ info = Ref {LinearAlgebra. BlasInt} ())
155
+ LinearAlgebra . require_one_based_indexing (A, ipiv, B)
154
156
LinearAlgebra. LAPACK. chktrans (trans)
155
157
chkstride1 (A, B, ipiv)
156
158
n = LinearAlgebra. checksquare (A)
@@ -162,20 +164,22 @@ function getrs!(trans::AbstractChar,
162
164
end
163
165
nrhs = size (B, 2 )
164
166
ccall ((@blasfunc (cgetrs_), OpenBLAS_jll. libopenblas), Cvoid,
165
- (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF32}, Ref{BlasInt},
166
- Ptr{BlasInt}, Ptr{ComplexF32}, Ref{BlasInt}, Ptr{BlasInt}, Clong),
167
+ (Ref{UInt8}, Ref{LinearAlgebra. BlasInt}, Ref{LinearAlgebra. BlasInt},
168
+ Ptr{ComplexF32}, Ref{LinearAlgebra. BlasInt},
169
+ Ptr{LinearAlgebra. BlasInt}, Ptr{ComplexF32}, Ref{LinearAlgebra. BlasInt},
170
+ Ptr{LinearAlgebra. BlasInt}, Clong),
167
171
trans, n, size (B, 2 ), A, max (1 , stride (A, 2 )), ipiv, B, max (1 , stride (B, 2 )), info,
168
172
1 )
169
- LinearAlgebra. LAPACK. chklapackerror (BlasInt (info[]))
173
+ LinearAlgebra. LAPACK. chklapackerror (LinearAlgebra . BlasInt (info[]))
170
174
B
171
175
end
172
176
173
177
function getrs! (trans:: AbstractChar ,
174
178
A:: AbstractMatrix{<:Float64} ,
175
- ipiv:: AbstractVector{BlasInt} ,
179
+ ipiv:: AbstractVector{LinearAlgebra. BlasInt} ,
176
180
B:: AbstractVecOrMat{<:Float64} ;
177
- info = Ref {BlasInt} ())
178
- require_one_based_indexing (A, ipiv, B)
181
+ info = Ref {LinearAlgebra. BlasInt} ())
182
+ LinearAlgebra . require_one_based_indexing (A, ipiv, B)
179
183
LinearAlgebra. LAPACK. chktrans (trans)
180
184
chkstride1 (A, B, ipiv)
181
185
n = LinearAlgebra. checksquare (A)
@@ -187,20 +191,21 @@ function getrs!(trans::AbstractChar,
187
191
end
188
192
nrhs = size (B, 2 )
189
193
ccall ((@blasfunc (dgetrs_), OpenBLAS_jll. libopenblas), Cvoid,
190
- (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
191
- Ptr{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}, Clong),
194
+ (Ref{UInt8}, Ref{LinearAlgebra. BlasInt}, Ref{LinearAlgebra. BlasInt},
195
+ Ptr{Float64}, Ref{LinearAlgebra. BlasInt},
196
+ Ptr{LinearAlgebra. BlasInt}, Ptr{Float64}, Ref{LinearAlgebra. BlasInt}, Ptr{LinearAlgebra. BlasInt}, Clong),
192
197
trans, n, size (B, 2 ), A, max (1 , stride (A, 2 )), ipiv, B, max (1 , stride (B, 2 )), info,
193
198
1 )
194
- LinearAlgebra. LAPACK. chklapackerror (BlasInt (info[]))
199
+ LinearAlgebra. LAPACK. chklapackerror (LinearAlgebra . BlasInt (info[]))
195
200
B
196
201
end
197
202
198
203
function getrs! (trans:: AbstractChar ,
199
204
A:: AbstractMatrix{<:Float32} ,
200
- ipiv:: AbstractVector{BlasInt} ,
205
+ ipiv:: AbstractVector{LinearAlgebra. BlasInt} ,
201
206
B:: AbstractVecOrMat{<:Float32} ;
202
- info = Ref {BlasInt} ())
203
- require_one_based_indexing (A, ipiv, B)
207
+ info = Ref {LinearAlgebra. BlasInt} ())
208
+ LinearAlgebra . require_one_based_indexing (A, ipiv, B)
204
209
LinearAlgebra. LAPACK. chktrans (trans)
205
210
chkstride1 (A, B, ipiv)
206
211
n = LinearAlgebra. checksquare (A)
@@ -212,11 +217,12 @@ function getrs!(trans::AbstractChar,
212
217
end
213
218
nrhs = size (B, 2 )
214
219
ccall ((@blasfunc (sgetrs_), OpenBLAS_jll. libopenblas), Cvoid,
215
- (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float32}, Ref{BlasInt},
216
- Ptr{BlasInt}, Ptr{Float32}, Ref{BlasInt}, Ptr{BlasInt}, Clong),
220
+ (Ref{UInt8}, Ref{LinearAlgebra. BlasInt}, Ref{LinearAlgebra. BlasInt},
221
+ Ptr{Float32}, Ref{LinearAlgebra. BlasInt},
222
+ Ptr{LinearAlgebra. BlasInt}, Ptr{Float32}, Ref{LinearAlgebra. BlasInt}, Ptr{LinearAlgebra. BlasInt}, Clong),
217
223
trans, n, size (B, 2 ), A, max (1 , stride (A, 2 )), ipiv, B, max (1 , stride (B, 2 )), info,
218
224
1 )
219
- LinearAlgebra. LAPACK. chklapackerror (BlasInt (info[]))
225
+ LinearAlgebra. LAPACK. chklapackerror (LinearAlgebra . BlasInt (info[]))
220
226
B
221
227
end
222
228
@@ -227,7 +233,7 @@ default_alias_b(::OpenBLASLUFactorization, ::Any, ::Any) = false
227
233
228
234
const PREALLOCATED_OPENBLAS_LU = begin
229
235
A = rand (0 , 0 )
230
- luinst = ArrayInterface. lu_instance (A), Ref {BlasInt} ()
236
+ luinst = ArrayInterface. lu_instance (A), Ref {LinearAlgebra. BlasInt} ()
231
237
end
232
238
233
239
function LinearSolve. init_cacheval (alg:: OpenBLASLUFactorization , A, b, u, Pl, Pr,
@@ -241,7 +247,7 @@ function LinearSolve.init_cacheval(alg::OpenBLASLUFactorization,
241
247
maxiters:: Int , abstol, reltol, verbose:: LinearVerbosity ,
242
248
assumptions:: OperatorAssumptions )
243
249
A = rand (eltype (A), 0 , 0 )
244
- ArrayInterface. lu_instance (A), Ref {BlasInt} ()
250
+ ArrayInterface. lu_instance (A), Ref {LinearAlgebra. BlasInt} ()
245
251
end
246
252
247
253
function SciMLBase. solve! (cache:: LinearCache , alg:: OpenBLASLUFactorization ;
@@ -251,7 +257,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization;
251
257
if cache. isfresh
252
258
cacheval = @get_cacheval (cache, :OpenBLASLUFactorization )
253
259
res = OpenBLASLU. getrf! (A; ipiv = cacheval[1 ]. ipiv, info = cacheval[2 ])
254
- fact = LU (res[1 : 3 ]. .. ), res[4 ]
260
+ fact = LinearAlgebra . LU (res[1 : 3 ]. .. ), res[4 ]
255
261
cache. cacheval = fact
256
262
257
263
if ! LinearAlgebra. issuccess (fact[1 ])
@@ -262,7 +268,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization;
262
268
end
263
269
264
270
A, info = @get_cacheval (cache, :OpenBLASLUFactorization )
265
- require_one_based_indexing (cache. u, cache. b)
271
+ LinearAlgebra . require_one_based_indexing (cache. u, cache. b)
266
272
m, n = size (A, 1 ), size (A, 2 )
267
273
if m > n
268
274
Bc = copy (cache. b)
@@ -273,5 +279,6 @@ function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization;
273
279
OpenBLASLU. getrs! (' N' , A. factors, A. ipiv, cache. u; info)
274
280
end
275
281
276
- SciMLBase. build_linear_solution (alg, cache. u, nothing , cache; retcode = ReturnCode. Success)
277
- end
282
+ SciMLBase. build_linear_solution (
283
+ alg, cache. u, nothing , cache; retcode = ReturnCode. Success)
284
+ end
0 commit comments