search

JuliaLinearAlgebra / BandedMatrices.jl

A Julia package for representing banded matrices
http://julialinearalgebra.github.io/BandedMatrices.jl/
Other
129 stars 38 forks source link

LU decomposition #272

Open eschnett opened 2 years ago

eschnett commented 2 years ago

I needed to invert large banded matrices, and I found that BandedMatrices.jl didn't support this operation. For example:

julia> using LinearAlgebra
julia> using BandedMatrices

julia> B = brand(Int8, 5, 5, 2, 2)
5×5 BandedMatrix{Int8} with bandwidths (2, 2):
 68  -17   -36     ⋅    ⋅
  2  -88  -118    83    ⋅
 83    8   -47    52   83
  ⋅   -6   -44  -116   98
  ⋅    ⋅    57  -104  -87

julia> lu(B)
BandedMatrices.BandedLU{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}
L factor:
5×5 Matrix{Float64}:
  1.0           0.0           0.0       0.0        0.0
  0.0240964     1.0           0.0       0.0        0.0
  5.7603e-16   -1.33319e-16   1.0       0.0        0.0
 -3.15648e-16   0.0680328    -0.632442  1.0        0.0
  0.819277      0.267077      0.591554  0.0155523  1.0
U factor:
5×5 BandedMatrix{Float64} with bandwidths (2, 4):
 83.0    8.0      -47.0      52.0     83.0
  0.0  -88.1928  -116.867    81.747   -2.0
  0.0    0.0       57.0    -104.0    -87.0
   ⋅     0.0        0.0    -187.335   43.1136
   ⋅      ⋅         0.0       0.0    -16.6712

This uses a full matrix to store L, which is not good for me.

Similarly, when using rational numbers, I found that

julia> B = Rational{BigInt}.(brand(Int8, 5, 5, 2, 2))
5×5 BandedMatrix{Rational{BigInt}} with bandwidths (2, 2):
 -128//1  -64//1   74//1      ⋅       ⋅
   27//1  113//1  115//1   -88//1     ⋅
  -10//1   87//1   51//1  -106//1  -67//1
     ⋅     -6//1  -48//1   -43//1  -47//1
     ⋅       ⋅     29//1   -93//1   71//1

julia> lu(B)
ERROR: UndefVarError: banded_lufact! not defined

I have taken the liberty of implementing the functionality that I need outside of BandedMatrices.jl. See here for the functions I implemented.

The respective functionality is good enough for me, but is not quite generic. For example, I implemented \ to solve a linear system, but not an in-place ldiv!. These functions also don't support a banded RHS for \.

I would be happy to contribute this to BandedMatrices.jl. If you are interested, then please let me know how this could/should be integrated into your package. Otherwise, feel free to copy the code I wrote. This code is inspired (copied/translated) from Wikipedia (see the URL above), LAPACK (dtrsm), and the Julia LinearAlgebra package.

dlfivefifty commented 2 years ago

Yes please! I believe the issue is that LU is currently only implemented in LAPACK so a Julia native implementation would be fantastic

eschnett commented 2 years ago

I've looked at the code for LU decomposition in BandedMatrices.jl, and it's not clear to me which part of the code would need to stay and which part should be rewritten. It seems that some LU decomposition functionality is there, but it's not implemented in an ideal way – e.g. it seems the L factor is full instead of banded.

dlfivefifty commented 2 years ago

I suggest you read the docs on LAPACK's gbtrf:

https://netlib.org/lapack/explore-html/d2/d2d/group__double_g_bcomputational_ga7fc91ba3f250ad3844eba25d59f5d7be.html

Is the type BandedLU necessary? It seems that the Base LU type would suffice.

Yes since gbtrf does not actually permute the rows in memory (as that would destroy banded structure), instead it just records where the pivots were. So it's not the same format as LinearAlgebra.LU.

Should the pivot search remain? A pivot search destroys the band structure. One could remove it, and ask people to convert to a full matrix if they want a pivot search.

Yes because you maintain bandedness if you do what LAPACK does.

I am interested in supporting types that are not supported by LAPACK. Would it make sense to remove the calls to LAPACK, implementing the respective functions in Julia?

No, still use LAPACK for BLAS types. Only use the Julia implementation for general types. You can use BigFloat as a test.

A comment in the function lu! states "l of the bands are ignored and overwritten". That seems weird. The problem seems to be that the LAPACK function dgbsv allows (and here requires) over-allocating storage, which is not supported by this package. If there is no functionality to resize banded matrices without de-allocating storage, then filling the extra bands with zeros might be a good approach?

I don't agree: we do allow overallocating storage, e.g., if the data is a (strided) view of a bigger array.

j-fu commented 1 year ago

Hi, FYI: I did some similar job for dgemm etc in sparspak: https://github.com/PetrKryslUCSD/Sparspak.jl/blob/main/src/Utilities/GenericBlasLapackFragments.jl

This was based on the netlib fortran code.