Skip to content

Commit 895d8ab

Browse files
YingboMaYingboMa
committed
init
0 parents  commit 895d8ab

File tree

8 files changed

+254
-0
lines changed

8 files changed

+254
-0
lines changed

‎.codecov.yml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
comment: false

‎.gitignore

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
*.jl.cov
2+
*.jl.*.cov
3+
*.jl.mem

‎.travis.yml

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
# Documentation: http://docs.travis-ci.com/user/languages/julia/
2+
language: julia
3+
os:
4+
- linux
5+
julia:
6+
- 1.0
7+
- nightly
8+
#matrix:
9+
# allow_failures:
10+
# - julia: nightly
11+
notifications:
12+
email: false
13+
after_success:
14+
# push coverage results to Coveralls
15+
- julia -e 'cd(Pkg.dir("RecursiveFactorization")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
16+
# push coverage results to Codecov
17+
- julia -e 'cd(Pkg.dir("RecursiveFactorization")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'

‎LICENSE.md

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
The PuMaS.jl package is licensed under the MIT "Expat" License:
2+
3+
> Copyright (c) 2019: Yingbo Ma.
4+
>
5+
>
6+
> Permission is hereby granted, free of charge, to any person obtaining a copy
7+
>
8+
> of this software and associated documentation files (the "Software"), to deal
9+
>
10+
> in the Software without restriction, including without limitation the rights
11+
>
12+
> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13+
>
14+
> copies of the Software, and to permit persons to whom the Software is
15+
>
16+
> furnished to do so, subject to the following conditions:
17+
>
18+
>
19+
>
20+
> The above copyright notice and this permission notice shall be included in all
21+
>
22+
> copies or substantial portions of the Software.
23+
>
24+
>
25+
>
26+
> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
27+
>
28+
> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
29+
>
30+
> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
31+
>
32+
> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
33+
>
34+
> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
35+
>
36+
> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
37+
>
38+
> SOFTWARE.
39+
>
40+
>

‎README.md

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
# RecursiveFactorization
2+
3+
---
4+
5+
`RecursiveFactorization.jl` is a package that collects various recursive matrix factorization algorithms.
6+
7+
#### Implemented Algorithms:
8+
9+
- Sivan Toledo's recursive left-looking LU algorithm. DOI: [10.1137/S0895479896297744](https://epubs.siam.org/doi/10.1137/S0895479896297744)

‎src/RecursiveFactorization.jl

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
module RecursiveFactorization
2+
3+
include("./lu.jl")
4+
5+
end # module

‎src/lu.jl

Lines changed: 150 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,150 @@
1+
using LinearAlgebra: BlasInt, LU, UnitLowerTriangular, ldiv!, BLAS, checknonsingular
2+
3+
lu(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}} = Val(true);
4+
check::Bool = true, blocksize::Integer = 16) = lu!(copy(A), pivot;
5+
check = check, blocksize = blocksize)
6+
7+
lu!(A, pivot::Union{Val{false}, Val{true}} = Val(true);
8+
check::Bool = true, blocksize::Integer = 16) = lu!(copy(A), Vector{BlasInt}(undef, min(size(A)...)), pivot;
9+
check = check, blocksize = blocksize)
10+
11+
function lu!(A::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer},
12+
pivot::Union{Val{false}, Val{true}} = Val(true);
13+
check::Bool=true, blocksize::Integer=16) where T
14+
info = Ref(zero(BlasInt))
15+
m, n = size(A)
16+
mnmin = min(m, n)
17+
reckernel!(A, pivot, m, mnmin, ipiv, info, blocksize)
18+
check && checknonsingular(info[])
19+
LU{T, typeof(A)}(A, ipiv, info[])
20+
end
21+
22+
nsplit(::Type{Float64}, n) = n >= 16 ? ((n + 8) ÷ 16) * 8 : n ÷ 2
23+
nsplit(::Type{Float32}, n) = n >= 32 ? ((n + 16) ÷ 32) * 16 : n ÷ 2
24+
nsplit(::Type{ComplexF64}, n) = n >= 8 ? ((n + 4) ÷ 8) * 4 : n ÷ 2
25+
nsplit(::Type{ComplexF32}, n) = n >= 16 ? ((n + 8) ÷ 16) * 8 : n ÷ 2
26+
27+
Base.@propagate_inbounds function apply_permutation!(P, A)
28+
for i in axes(P, 1)
29+
i′ = P[i]
30+
i′ == i && continue
31+
@simd for j in axes(A, 2)
32+
A[i, j], A[i′, j] = A[i′, j], A[i, j]
33+
end
34+
end
35+
nothing
36+
end
37+
38+
function reckernel!(A::AbstractMatrix{T}, pivot::Val{Pivot}, m, n, ipiv, info, blocksize)::Nothing where {T,Pivot}
39+
@inbounds begin
40+
if n <= max(blocksize, 1)
41+
_generic_lufact!(A, pivot, ipiv, info)
42+
return nothing
43+
end
44+
n1 = nsplit(T, n)
45+
n2 = n - n1
46+
m2 = m - n1
47+
48+
# ======================================== #
49+
# Now, our LU process looks like this
50+
# [ P1 ] [ A11 A21 ] [ L11 0 ] [ U11 U12 ]
51+
# [ ] [ ] = [ ] [ ]
52+
# [ P2 ] [ A21 A22 ] [ L21 I ] [ 0 A′22 ]
53+
# ======================================== #
54+
55+
# ======================================== #
56+
# Partition the matrix A
57+
# [AL AR]
58+
AL = @view A[:, 1:n1]
59+
AR = @view A[:, n1+1:n]
60+
# AL AR
61+
# [A11 A12]
62+
# [A21 A22]
63+
A11 = @view A[1:n1, 1:n1]
64+
A12 = @view A[1:n1, n1+1:n]
65+
A21 = @view A[n1+1:m, 1:n1]
66+
A22 = @view A[n1+1:m, n1+1:n]
67+
# [P1]
68+
# [P2]
69+
P1 = @view ipiv[1:n1]
70+
P2 = @view ipiv[n1+1:n]
71+
# ========================================
72+
73+
# [ A11 ] [ L11 ]
74+
# P [ ] = [ ] U11
75+
# [ A21 ] [ L21 ]
76+
reckernel!(AL, pivot, m, n1, P1, info, blocksize)
77+
# [ A12 ] [ P1 ] [ A12 ]
78+
# [ ] <- [ ] [ ]
79+
# [ A22 ] [ 0 ] [ A22 ]
80+
Pivot && apply_permutation!(P1, AR)
81+
# A12 = L11 U12 => U12 = L11 \ A12
82+
ldiv!(UnitLowerTriangular(A11), A12)
83+
# Schur complement:
84+
# We have A22 = L21 U12 + A′22, hence
85+
# A′22 = A22 - L21 U12
86+
BLAS.gemm!('N', 'N', -one(T), A21, A12, one(T), A22)
87+
# record info
88+
previnfo = info[]
89+
# P2 A22 = L22 U22
90+
reckernel!(A22, pivot, m2, n2, P2, info, blocksize)
91+
# A21 <- P2 A21
92+
Pivot && apply_permutation!(P2, A21)
93+
94+
info[] != previnfo && (info[] += n1)
95+
@simd for i in 1:n2
96+
P2[i] += n1
97+
end
98+
return nothing
99+
end # inbounds
100+
end
101+
102+
#=
103+
Modified from https://github.com/JuliaLang/julia/blob/b56a9f07948255dfbe804eef25bdbada06ec2a57/stdlib/LinearAlgebra/src/lu.jl
104+
License is MIT: https://julialang.org/license
105+
=#
106+
function _generic_lufact!(A::StridedMatrix{T}, ::Val{Pivot}, ipiv, info) where {Pivot,T}
107+
m, n = size(A)
108+
minmn = length(ipiv)
109+
@inbounds begin
110+
for k = 1:minmn
111+
# find index max
112+
kp = k
113+
if Pivot
114+
amax = abs(zero(T))
115+
for i = k:m
116+
absi = abs(A[i,k])
117+
if absi > amax
118+
kp = i
119+
amax = absi
120+
end
121+
end
122+
end
123+
ipiv[k] = kp
124+
if !iszero(A[kp,k])
125+
if k != kp
126+
# Interchange
127+
@simd for i = 1:n
128+
tmp = A[k,i]
129+
A[k,i] = A[kp,i]
130+
A[kp,i] = tmp
131+
end
132+
end
133+
# Scale first column
134+
Akkinv = inv(A[k,k])
135+
@simd for i = k+1:m
136+
A[i,k] *= Akkinv
137+
end
138+
elseif info[] == 0
139+
info[] = k
140+
end
141+
# Update the rest
142+
for j = k+1:n
143+
@simd for i = k+1:m
144+
A[i,j] -= A[i,k]*A[k,j]
145+
end
146+
end
147+
end
148+
end
149+
return nothing
150+
end

‎test/runtests.jl

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
using Test
2+
import RecursiveFactorization
3+
import LinearAlgebra
4+
5+
baselu = LinearAlgebra.lu
6+
mylu = RecursiveFactorization.lu
7+
8+
function testlu(MF, BF)
9+
@test MF.info == BF.info
10+
@test MF.L ≈ BF.L
11+
@test MF.U ≈ BF.U
12+
@test MF.P ≈ BF.P
13+
nothing
14+
end
15+
16+
@testset "Test LU factorization" begin
17+
for p in (Val(true), Val(false))
18+
A = rand(100, 100)
19+
MF = mylu(A, p)
20+
BF = baselu(A, p)
21+
testlu(MF, BF)
22+
for i in 1:100
23+
A[:, i] .= 0
24+
MF = mylu(A, p, check=false)
25+
BF = baselu(A, p, check=false)
26+
testlu(MF, BF)
27+
end
28+
end
29+
end

0 commit comments

Comments
 (0)