@@ -22,37 +22,13 @@ function lu(A::AbstractMatrix, pivot = Val(true), thread = Val(true); kwargs...)
22
22
return lu! (copy (A), normalize_pivot (pivot), thread; kwargs... )
23
23
end
24
24
25
- struct NotIPIV <: AbstractVector{BlasInt}
26
- len:: Int
27
- end
28
- Base. size (A:: NotIPIV ) = (A. len,)
29
- Base. getindex (:: NotIPIV , i:: Int ) = i
30
- Base. view (:: NotIPIV , r:: AbstractUnitRange ) = NotIPIV (length (r))
31
- init_pivot (:: Val{false} , minmn) = NotIPIV (minmn)
32
- init_pivot (:: Val{true} , minmn) = Vector {BlasInt} (undef, minmn)
33
-
34
- if isdefined (LinearAlgebra, :_ipiv_cols! )
35
- function LinearAlgebra. _ipiv_cols! (:: LU{<:Any, <:Any, NotIPIV} , :: OrdinalRange ,
36
- B:: StridedVecOrMat )
37
- return B
38
- end
39
- end
40
- if isdefined (LinearAlgebra, :_ipiv_rows! )
41
- function LinearAlgebra. _ipiv_rows! (:: LU{<:Any, <:Any, NotIPIV} , :: OrdinalRange ,
42
- B:: StridedVecOrMat )
43
- return B
44
- end
45
- end
46
-
47
25
function lu! (A, pivot = Val (true ), thread = Val (true ); check = true , kwargs... )
48
26
m, n = size (A)
49
27
minmn = min (m, n)
50
- # we want the type on both branches to match. When pivot = Val(false), we construct
51
- # a `NotIPIV`, which `LinearAlgebra.generic_lufact!` does not.
52
- F = if pivot === Val (true ) && minmn < 10 # avx introduces small performance degradation
28
+ F = if minmn < 10 # avx introduces small performance degradation
53
29
LinearAlgebra. generic_lufact! (A, to_stdlib_pivot (pivot); check = check)
54
30
else
55
- lu! (A, init_pivot (pivot , minmn), normalize_pivot (pivot), thread; check = check,
31
+ lu! (A, Vector {BlasInt} (undef , minmn), normalize_pivot (pivot), thread; check = check,
56
32
kwargs... )
57
33
end
58
34
return F
@@ -68,8 +44,6 @@ pick_threshold() = LoopVectorization.register_size() == 64 ? 48 : 40
68
44
recurse (:: StridedArray ) = true
69
45
recurse (_) = false
70
46
71
- _ptrarray (ipiv) = PtrArray (ipiv)
72
- _ptrarray (ipiv:: NotIPIV ) = ipiv
73
47
function lu! (A:: AbstractMatrix{T} , ipiv:: AbstractVector{<:Integer} ,
74
48
pivot = Val (true ), thread = Val (true );
75
49
check:: Bool = true ,
@@ -84,7 +58,7 @@ function lu!(A::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer},
84
58
if T <: Union{Float32, Float64}
85
59
GC. @preserve ipiv A begin info = recurse! (view (PtrArray (A), axes (A)... ), pivot,
86
60
m, n, mnmin,
87
- _ptrarray (ipiv), info, blocksize,
61
+ PtrArray (ipiv), info, blocksize,
88
62
thread) end
89
63
else
90
64
info = recurse! (A, pivot, m, n, mnmin, ipiv, info, blocksize, thread)
116
90
# [AL AR]
117
91
AL = @view A[:, 1 : m]
118
92
AR = @view A[:, (m + 1 ): n]
119
- apply_permutation! (ipiv, AR, Val {Thread} ( ))
120
- ldiv! (_unit_lower_triangular (AL), AR, Val {Thread} ( ))
93
+ apply_permutation! (ipiv, AR, Val (Thread ))
94
+ ldiv! (_unit_lower_triangular (AL), AR, Val (Thread ))
121
95
end
122
96
info
123
97
end
@@ -213,10 +187,8 @@ function reckernel!(A::AbstractMatrix{T}, pivot::Val{Pivot}, m, n, ipiv, info, b
213
187
Pivot && apply_permutation! (P2, A21, thread)
214
188
215
189
info != previnfo && (info += n1)
216
- if Pivot
217
- @turbo warn_check_args= false for i in 1 : n2
218
- P2[i] += n1
219
- end
190
+ @turbo warn_check_args= false for i in 1 : n2
191
+ P2[i] += n1
220
192
end
221
193
return info
222
194
end # inbounds
@@ -262,8 +234,8 @@ function _generic_lufact!(A, ::Val{Pivot}, ipiv, info) where {Pivot}
262
234
amax = absi
263
235
end
264
236
end
265
- ipiv[k] = kp
266
237
end
238
+ ipiv[k] = kp
267
239
if ! iszero (A[kp, k])
268
240
if k != kp
269
241
# Interchange
0 commit comments