Skip to content

fixed example in readme #21

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Apr 3, 2025
Merged

fixed example in readme #21

Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
6e2ff84
fixed example in readme
ajinkya-k Apr 2, 2025
ee320a8
fixed hermitian example
ajinkya-k Apr 2, 2025
a95f319
using floats everywhere
ajinkya-k Apr 3, 2025
f94b0d2
removed extra spaces
ajinkya-k Apr 3, 2025
b7c1ab7
removed float wrapper
ajinkya-k Apr 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
126 changes: 63 additions & 63 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,49 +13,49 @@ For example, starting with a 6 by 6 matrix whose elements are numbered 1 to 36 i
```julia
julia> using LinearAlgebra, RectangularFullPacked

julia> A = reshape(1:36, (6, 6))
6×6 reshape(::UnitRange{Int64}, 6, 6) with eltype Int64:
1 7 13 19 25 31
2 8 14 20 26 32
3 9 15 21 27 33
4 10 16 22 28 34
5 11 17 23 29 35
6 12 18 24 30 36
julia> A = reshape(1.:36., (6, 6))
6×6 reshape(::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, 6, 6) with eltype Float64:
1.0 7.0 13.0 19.0 25.0 31.0
2.0 8.0 14.0 20.0 26.0 32.0
3.0 9.0 15.0 21.0 27.0 33.0
4.0 10.0 16.0 22.0 28.0 34.0
5.0 11.0 17.0 23.0 29.0 35.0
6.0 12.0 18.0 24.0 30.0 36.0
```
the lower triangular matrix `AL` is constructed by replacing the elements above the diagonal with zero.
```julia
julia> AL = tril!(collect(A))
6×6 Matrix{Int64}:
1 0 0 0 0 0
2 8 0 0 0 0
3 9 15 0 0 0
4 10 16 22 0 0
5 11 17 23 29 0
6 12 18 24 30 36
```
`AL` requires the same amount of storage as does `A` even though there are only 21 potential non-zeros in `AL`.
6×6 Matrix{Float64}:
1.0 0.0 0.0 0.0 0.0 0.0
2.0 8.0 0.0 0.0 0.0 0.0
3.0 9.0 15.0 0.0 0.0 0.0
4.0 10.0 16.0 22.0 0.0 0.0
5.0 11.0 17.0 23.0 29.0 0.0
6.0 12.0 18.0 24.0 30.0 36.0
```
`AL` requires the same amount of storage as does `A` even though there are only 21 potential non-zeros in `AL`.
The RFP version of the lower triangular matrix
```julia
julia> ArfpL = Int.(TriangularRFP(float.(A), :L))
6×6 Matrix{Int64}:
1 0 0 0 0 0
2 8 0 0 0 0
3 9 15 0 0 0
4 10 16 22 0 0
5 11 17 23 29 0
6 12 18 24 30 36
```julia
julia> ArfpL = TriangularRFP(float.(A), :L)
6×6 TriangularRFP{Float64}:
1.0 0.0 0.0 0.0 0.0 0.0
2.0 8.0 0.0 0.0 0.0 0.0
3.0 9.0 15.0 0.0 0.0 0.0
4.0 10.0 16.0 22.0 0.0 0.0
5.0 11.0 17.0 23.0 29.0 0.0
6.0 12.0 18.0 24.0 30.0 36.0
```
provides the same displayed form but the underlying, "parent" array is 7 by 3
```julia
julia> ALparent = Int.(ArfpL.data)
7×3 Matrix{Int64}:
22 23 24
1 29 30
2 8 36
3 9 15
4 10 16
5 11 17
6 12 18
provides the same displayed form but the underlying, "parent" array is 7 by 3
```julia
julia> ALparent = ArfpL.data
7×3 Matrix{Float64}:
22.0 23.0 24.0
1.0 29.0 30.0
2.0 8.0 36.0
3.0 9.0 15.0
4.0 10.0 16.0
5.0 11.0 17.0
6.0 12.0 18.0
```

The three blocks of `AL` are the lower triangle of `AL[1:3, 1:3]`, stored as the lower triangle of `ALparent[2:4, :]`; the square block `AL[4:6, 1:3]` stored in `ALparent[5:7, :]`; and the lower triangle of `AL[4:6, 4:6]` stored transposed in `ALparent[1:3, :]`.
Expand All @@ -64,50 +64,50 @@ For odd values of n, the parent is of size `(n, div(n + 1, 2))` and the non-zero

For example,
```julia
julia> AL = tril!(collect(reshape(1:25, 5, 5)))
5×5 Matrix{Int64}:
1 0 0 0 0
2 7 0 0 0
3 8 13 0 0
4 9 14 19 0
5 10 15 20 25
julia> AL = tril!(collect(reshape(1.:25., 5, 5)))
5×5 Matrix{Float64}:
1.0 0.0 0.0 0.0 0.0
2.0 7.0 0.0 0.0 0.0
3.0 8.0 13.0 0.0 0.0
4.0 9.0 14.0 19.0 0.0
5.0 10.0 15.0 20.0 25.0

julia> ArfpL = Int.(TriangularRFP(float(AL), :L).data)
5×3 Matrix{Int64}:
1 19 20
2 7 25
3 8 13
4 9 14
5 10 15
julia> ArfpL = TriangularRFP(float(AL), :L).data
5×3 Matrix{Float64}:
1.0 19.0 20.0
2.0 7.0 25.0
3.0 8.0 13.0
4.0 9.0 14.0
5.0 10.0 15.0
```

RFP storage is especially useful for large positive definite Hermitian matrices because the Cholesky factor can be evaluated nearly as quickly (by applying Level-3 BLAS to the blocks) as in full storage mode but requiring about half the storage.
RFP storage is especially useful for large positive definite Hermitian matrices because the Cholesky factor can be evaluated nearly as quickly (by applying Level-3 BLAS to the blocks) as in full storage mode but requiring about half the storage.

A trivial example is
```julia
julia> A = [2 1 2; 1 2 0; 1 0 2]
3×3 Matrix{Int64}:
2 1 2
1 2 0
1 0 2
A trivial example is
```julia
julia> A = [2. 1 2; 1 2 0; 1 0 2]
3×3 Matrix{Float64}:
2.0 1.0 2.0
1.0 2.0 0.0
1.0 0.0 2.0

julia> cholesky(Hermitian(A, :L))
julia> cholesky(Hermitian(A, :L))
Cholesky{Float64, Matrix{Float64}}
L factor:
3×3 LowerTriangular{Float64, Matrix{Float64}}:
1.41421 ⋅ ⋅
0.707107 1.22474 ⋅
0.707107 -0.408248 1.1547

julia> ArfpL = HermitianRFP(TriangularRFP(float.(A), :L))
julia> ArfpL = Hermitian(TriangularRFP(float.(A), :L))
3×3 HermitianRFP{Float64}:
2.0 1.0 1.0
1.0 2.0 0.0
1.0 0.0 2.0

julia> cholesky!(ArfpL).data
julia> cholesky!(ArfpL).data
3×2 Matrix{Float64}:
1.41421 1.1547
0.707107 1.22474
0.707107 -0.408248
```
```
Loading