@@ -33,64 +33,67 @@ _lu(A::StaticMatrix{1,N,T}, ::Type{Val{Pivot}}) where {N,T,Pivot} =
33
33
(SMatrix {1,1,T} (one (T)), A, SVector {1,Int} (1 ))
34
34
35
35
function _lu (A:: StaticMatrix{M,1} , :: Type{Val{Pivot}} ) where {M,Pivot}
36
- kp = 1
37
- if Pivot
38
- amax = abs (A[1 ,1 ])
39
- for i = 2 : M
40
- absi = abs (A[i,1 ])
41
- if absi > amax
42
- kp = i
43
- amax = absi
36
+ @inbounds begin
37
+ kp = 1
38
+ if Pivot
39
+ amax = abs (A[1 ,1 ])
40
+ for i = 2 : M
41
+ absi = abs (A[i,1 ])
42
+ if absi > amax
43
+ kp = i
44
+ amax = absi
45
+ end
44
46
end
45
47
end
48
+ ps = tailindices (Val{M})
49
+ if kp != 1
50
+ ps = setindex (ps, 1 , kp- 1 )
51
+ end
52
+ U = SMatrix {1,1} (A[kp,1 ])
53
+ # Scale first column
54
+ Akkinv = inv (A[kp,1 ])
55
+ Ls = A[ps,1 ] * Akkinv
56
+ if ! isfinite (Akkinv)
57
+ Ls = zeros (Ls)
58
+ end
59
+ L = [SVector {1} (one (eltype (Ls))); Ls]
60
+ p = [SVector {1,Int} (kp); ps]
46
61
end
47
- ps = tailindices (Val{M})
48
- if kp != 1
49
- ps = setindex (ps, 1 , kp- 1 )
50
- end
51
- U = SMatrix {1,1} (A[kp,1 ])
52
- # Scale first column
53
- Akkinv = inv (A[kp,1 ])
54
- Ls = A[ps,1 ] * Akkinv
55
- if ! isfinite (Akkinv)
56
- Ls = zeros (Ls)
57
- end
58
- L = [SVector {1} (one (eltype (Ls))); Ls]
59
- p = [SVector {1,Int} (kp); ps]
60
62
return (SMatrix {M,1} (L), U, p)
61
63
end
62
64
63
65
function _lu (A:: StaticMatrix{M,N,T} , :: Type{Val{Pivot}} ) where {M,N,T,Pivot}
64
- kp = 1
65
- if Pivot
66
- amax = abs (A[1 ,1 ])
67
- for i = 2 : M
68
- absi = abs (A[i,1 ])
69
- if absi > amax
70
- kp = i
71
- amax = absi
66
+ @inbounds begin
67
+ kp = 1
68
+ if Pivot
69
+ amax = abs (A[1 ,1 ])
70
+ for i = 2 : M
71
+ absi = abs (A[i,1 ])
72
+ if absi > amax
73
+ kp = i
74
+ amax = absi
75
+ end
72
76
end
73
77
end
74
- end
75
- ps = tailindices (Val{M})
76
- if kp != 1
77
- ps = setindex (ps, 1 , kp- 1 )
78
- end
79
- Ufirst = SMatrix {1,N} (A[kp,:])
80
- # Scale first column
81
- Akkinv = inv (A[kp,1 ])
82
- Ls = A[ps,1 ] * Akkinv
83
- if ! isfinite (Akkinv)
84
- Ls = zeros (Ls)
85
- end
86
-
87
- # Update the rest
88
- Arest = A[ps,tailindices (Val{N})] - Ls* Ufirst[:,tailindices (Val{N})]
89
- Lrest, Urest, prest = _lu (Arest, Val{Pivot})
90
- p = [SVector {1,Int} (kp); ps[prest]]
91
- L = [[SVector {1} (one (eltype (Ls))); Ls[prest]] [zeros (SMatrix {1} (Lrest[1 ,:])); Lrest]]
92
- U = [Ufirst; [zeros (Urest[:,1 ]) Urest]]
78
+ ps = tailindices (Val{M})
79
+ if kp != 1
80
+ ps = setindex (ps, 1 , kp- 1 )
81
+ end
82
+ Ufirst = SMatrix {1,N} (A[kp,:])
83
+ # Scale first column
84
+ Akkinv = inv (A[kp,1 ])
85
+ Ls = A[ps,1 ] * Akkinv
86
+ if ! isfinite (Akkinv)
87
+ Ls = zeros (Ls)
88
+ end
93
89
90
+ # Update the rest
91
+ Arest = A[ps,tailindices (Val{N})] - Ls* Ufirst[:,tailindices (Val{N})]
92
+ Lrest, Urest, prest = _lu (Arest, Val{Pivot})
93
+ p = [SVector {1,Int} (kp); ps[prest]]
94
+ L = [[SVector {1} (one (eltype (Ls))); Ls[prest]] [zeros (SMatrix {1} (Lrest[1 ,:])); Lrest]]
95
+ U = [Ufirst; [zeros (Urest[:,1 ]) Urest]]
96
+ end
94
97
return (L, U, p)
95
98
end
96
99
0 commit comments