Merge pull request #17 from louisponet/main

Resizing support

Project.toml

@@ -1,7 +1,7 @@
 name = "FastLapackInterface"
 uuid = "29a986be-02c6-4525-aec4-84b980013641"
 authors = ["Louis Ponet, Michel Juillard"]
-version = "1.1.0"
+version = "1.2.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

README.md

@@ -24,4 +24,4 @@ decompose!(ws, A)
 A similar API exists for the above decompositions. For more information and examples please see the [documentation](https://dynarejulia.github.io/FastLapackInterface.jl/dev/).
 
 ## Package version
--   v1.x: works only with Julia >= 1.7
+-   v1.x: works only with Julia >= 1.6.3

benchmark/bench_qr.jl

@@ -90,21 +90,21 @@ for n in sizes
     suite["QRWYWs"]["geqrt!"]["LAPACK"]["$n"]    = @benchmarkable bench_geqrt!($As, $T)
 end
 
-###### QRpWs
+###### QRPivotedWs
 
-suite["QRpWs"] = BenchmarkGroup()
-suite["QRpWs"]["creation"] = BenchmarkGroup()
+suite["QRPivotedWs"] = BenchmarkGroup()
+suite["QRPivotedWs"]["creation"] = BenchmarkGroup()
 for n in sizes
-    suite["QRpWs"]["creation"]["$n"] = @benchmarkable QRWYWs(A) setup = (A = rand($n, $n))
+    suite["QRPivotedWs"]["creation"]["$n"] = @benchmarkable QRWYWs(A) setup = (A = rand($n, $n))
 end
 
 for n in sizes
-    suite["QRpWs"]["creation"]["$n"] = @benchmarkable QRWYWs(A) setup = (A = rand($n, $n))
+    suite["QRPivotedWs"]["creation"]["$n"] = @benchmarkable QRWYWs(A) setup = (A = rand($n, $n))
 end
 
-suite["QRpWs"]["geqp3!"]              = BenchmarkGroup()
-suite["QRpWs"]["geqp3!"]["workspace"] = BenchmarkGroup()
-suite["QRpWs"]["geqp3!"]["LAPACK"]    = BenchmarkGroup()
+suite["QRPivotedWs"]["geqp3!"]              = BenchmarkGroup()
+suite["QRPivotedWs"]["geqp3!"]["workspace"] = BenchmarkGroup()
+suite["QRPivotedWs"]["geqp3!"]["LAPACK"]    = BenchmarkGroup()
 
 function bench_geqp3!(As, ws)
     for A in As
@@ -121,10 +121,10 @@ for n in sizes
     As   = [rand(n, n) for i in 1:vector_length]
     τ    = zeros(n)
     lpvt = zeros(Int, n)
-    ws   = QRpWs(As[1])
+    ws   = QRPivotedWs(As[1])
 
-    suite["QRpWs"]["geqp3!"]["workspace"]["$n"] = @benchmarkable bench_geqp3!($As, $ws)
-    suite["QRpWs"]["geqp3!"]["LAPACK"]["$n"]    = @benchmarkable bench_geqp3!($As, $lpvt, $τ)
+    suite["QRPivotedWs"]["geqp3!"]["workspace"]["$n"] = @benchmarkable bench_geqp3!($As, $ws)
+    suite["QRPivotedWs"]["geqp3!"]["LAPACK"]["$n"]    = @benchmarkable bench_geqp3!($As, $lpvt, $τ)
 end
 
 end

docs/src/LAPACK.md

@@ -1,4 +1,5 @@
 # [LAPACK](@id LAPACK)
+This section details the [`LAPACK functions`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.LAPACK) that are supported for use with various [`Workspaces`](@ref WorkSpaces). Each function has a `resize` keyword argument that is `true` by default, allowing for automatic resizing of the workspaces to accomodate larger Matrices or different features than they were originally constructed for.
 
 ## Unified Interface
 After having created the [`Workspace`](@ref WorkSpaces) that corresponds to the targeted factorization or decomposition,
@@ -14,7 +15,7 @@ factorize!
 LinearAlgebra.LAPACK.geqrf!(::QRWs, ::AbstractMatrix)
 LinearAlgebra.LAPACK.ormqr!(::QRWs, ::AbstractChar, ::AbstractChar, ::AbstractMatrix, ::AbstractVecOrMat)
 LinearAlgebra.LAPACK.geqrt!(::QRWYWs, ::AbstractMatrix)
-LinearAlgebra.LAPACK.geqp3!(::QRpWs, ::AbstractMatrix)
+LinearAlgebra.LAPACK.geqp3!(::QRPivotedWs, ::AbstractMatrix)
 ```
 
 ## Schur

docs/src/workspaces.md

@@ -8,13 +8,17 @@ This can then be used to perform the decompositions without extra allocations.
 ```@docs
 Workspace
 ```
+Each [`Workspace`](@ref) also has a function to [`resize!`](@ref) to allow for its use with larger matrices or with more features (e.g. the computation of left eigenvectors and right eigenvectors using [`EigenWs`](@ref)).
+```@docs
+resize!(::Workspace, ::AbstractMatrix; kwargs...)
+```
 
 ## [QR](@id QR-id)
 
 ```@docs
 QRWs
 QRWYWs
-QRpWs
+QRPivotedWs
 ```
 
 ## [Schur](@id Schur-id)

src/FastLapackInterface.jl

@@ -20,7 +20,7 @@ abstract type Workspace end
 include("lu.jl")
 export LUWs
 include("qr.jl")
-export QRWs, QRWYWs, QRpWs
+export QRWs, QRWYWs, QRPivotedWs
 include("schur.jl")
 export SchurWs, GeneralizedSchurWs
 include("eigen.jl")

src/bunch_kaufman.jl

@@ -79,35 +79,44 @@ for (sytrfs,  elty) in
     (((:dsytrf_,:dsytrf_rook_),:Float64),
      ((:ssytrf_,:ssytrf_rook_),:Float32),
      ((:zsytrf_,:zsytrf_rook_, :zhetrf_,:zhetrf_rook_),:ComplexF64),
-     ((:csytrf_,:csytrf_rook_, :chetrf_,:chetrf_rook_),:ComplexF32))  
-    @eval function BunchKaufmanWs(A::AbstractMatrix{$elty})
-        chkstride1(A)
-        n = checksquare(A)
-        ipiv  = similar(A, BlasInt, n)
-        if n == 0
-            return BunchKaufmanWs($elty[], ipiv)
+     ((:csytrf_,:csytrf_rook_, :chetrf_,:chetrf_rook_),:ComplexF32))
+
+    @eval begin
+        function Base.resize!(ws::BunchKaufmanWs, A::AbstractMatrix{$elty})
+            chkstride1(A)
+            n = checksquare(A)
+            if n == 0
+                return ws
+            end
+            resize!(ws.ipiv, n)
+            info  = Ref{BlasInt}()
+            ccall((@blasfunc($(sytrfs[1])), liblapack), Cvoid,
+                  (Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
+                   Ptr{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Clong),
+                  'U', n, A, stride(A,2), ws.ipiv, ws.work, -1, info, 1)
+            chkargsok(info[])
+            resize!(ws.work, BlasInt(real(ws.work[1])))
+            return ws
+        end
+        function BunchKaufmanWs(A::AbstractMatrix{$elty})
+            return resize!(BunchKaufmanWs(Vector{$elty}(undef, 1), BlasInt[]), A)
         end
-        work  = Vector{$elty}(undef, 1)
-        lwork = BlasInt(-1)
-        info  = Ref{BlasInt}()
-        ccall((@blasfunc($(sytrfs[1])), liblapack), Cvoid,
-              (Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
-               Ptr{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Clong),
-              'U', n, A, stride(A,2), ipiv, work, lwork, info, 1)
-        chkargsok(info[])
-        resize!(work, BlasInt(real(work[1])))
-        return BunchKaufmanWs(work, ipiv)
     end
     for (sytrf, fn) in zip(sytrfs, (:sytrf!, :sytrf_rook!, :hetrf!, :hetrf_rook!))
-        @eval function $fn(ws::BunchKaufmanWs{$elty}, uplo::AbstractChar, A::AbstractMatrix{$elty})
+        @eval function $fn(ws::BunchKaufmanWs{$elty}, uplo::AbstractChar, A::AbstractMatrix{$elty}; resize=true)
             chkstride1(A)
             n = checksquare(A)
-            chkuplo(uplo)
-            lipiv = length(ws.ipiv)
-            @assert n <= lipiv "Workspace was allocated for matrices of maximum size ($lipiv, $lipiv)." 
             if n == 0
                 return A, ws.ipiv, zero(BlasInt)
             end
+            chkuplo(uplo)
+            if n > length(ws.ipiv)
+                if resize
+                    resize!(ws, A)
+                else
+                    throw(ArgumentError("Workspace is too small, use resize!(ws, A)."))
+                end
+            end
             info  = Ref{BlasInt}()
             ccall((@blasfunc($sytrf), liblapack), Cvoid,
                   (Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
@@ -120,10 +129,11 @@ for (sytrfs,  elty) in
 end
 
 """
-    sytrf!(ws, uplo, A) -> (A, ws.ipiv, info)
+    sytrf!(ws, uplo, A; resize=true) -> (A, ws.ipiv, info)
 
 Computes the Bunch-Kaufman factorization of a symmetric matrix `A`,
 using previously allocated workspace `ws`.
+If the workspace was too small and `resize==true` it will automatically resized.
 If `uplo = U`, the upper half of `A` is stored. If `uplo = L`, the lower
 half is stored.
 
@@ -132,25 +142,25 @@ the error code `info` which is a non-negative integer. If `info` is positive
 the matrix is singular and the diagonal part of the factorization is exactly
 zero at position `info`.
 """
-sytrf!(ws::BunchKaufmanWs, uplo::AbstractChar, A::AbstractMatrix)
+sytrf!(ws::BunchKaufmanWs, uplo::AbstractChar, A::AbstractMatrix; kwargs...)
 
 """
-    sytrf_rook!(ws, uplo, A) -> (A, ws.ipiv, info)
+    sytrf_rook!(ws, uplo, A; resize=true) -> (A, ws.ipiv, info)
 
 Similar to [`sytrf!`](@ref) but using the bounded ("rook") diagonal pivoting method.
 """
-sytrf_rook!(ws::BunchKaufmanWs, uplo::AbstractChar, A::AbstractMatrix)
+sytrf_rook!(ws::BunchKaufmanWs, uplo::AbstractChar, A::AbstractMatrix; kwargs...)
 
 """
-    hetrf!(ws, uplo, A) -> (A, ws.ipiv, info)
+    hetrf!(ws, uplo, A; resize=true) -> (A, ws.ipiv, info)
 
 Similar as [`sytrf!`](@ref) but for Hermitian matrices.
 """
-hetrf!(ws::BunchKaufmanWs, uplo::AbstractChar, A::AbstractMatrix)
+hetrf!(ws::BunchKaufmanWs, uplo::AbstractChar, A::AbstractMatrix; kwargs...)
 
 """
-    hetrf_rook!(ws, uplo, A) -> (A, ws.ipiv, info)
+    hetrf_rook!(ws, uplo, A; resize=true) -> (A, ws.ipiv, info)
 
 Similar to [`hetrf!`](@ref) but using the bounded ("rook") diagonal pivoting method.
 """
-hetrf_rook!(ws::BunchKaufmanWs, uplo::AbstractChar, A::AbstractMatrix)
+hetrf_rook!(ws::BunchKaufmanWs, uplo::AbstractChar, A::AbstractMatrix; kwargs...)

src/cholesky.jl

@@ -45,17 +45,30 @@ for (pstrf, elty, rtyp) in
      (:zpstrf_,:ComplexF64,:Float64),
      (:cpstrf_,:ComplexF32,:Float32))
     @eval begin
-        function CholeskyPivotedWs(A::AbstractMatrix{$elty})
+        function Base.resize!(ws::CholeskyPivotedWs, A::AbstractMatrix{$elty})
             n = checksquare(A)
-            return CholeskyPivotedWs(Vector{$rtyp}(undef, 2n), similar(A, BlasInt, n))
+            resize!(ws.work, 2n)
+            resize!(ws.piv, n)
+            return ws
+        end
+        function CholeskyPivotedWs(A::AbstractMatrix{$elty})
+            return resize!(CholeskyPivotedWs($rtyp[], BlasInt[]), A)
         end
 
-        function pstrf!(ws::CholeskyPivotedWs, uplo::AbstractChar, A::AbstractMatrix{$elty}, tol::Real)
+        function pstrf!(ws::CholeskyPivotedWs, uplo::AbstractChar, A::AbstractMatrix{$elty}, tol::Real; resize=true)
             chkstride1(A)
             n = checksquare(A)
             chkuplo(uplo)
             rank = Ref{BlasInt}()
             info = Ref{BlasInt}()
+            if length(ws.piv) < n
+                if resize
+                    resize!(ws, A)
+                else
+                    throw(ArgumentError("Workspace is too small, use resize!(ws, A)."))
+                end
+            end
+
             ccall((@blasfunc($pstrf), liblapack), Cvoid,
                   (Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt},
                    Ptr{BlasInt}, Ref{$rtyp}, Ptr{$rtyp}, Ptr{BlasInt}, Clong),
@@ -67,16 +80,17 @@ for (pstrf, elty, rtyp) in
 end
 
 """
-    pstrf!(ws, uplo, A, tol) -> (A, ws.piv, rank, info)
+    pstrf!(ws, uplo, A, tol; resize=true) -> (A, ws.piv, rank, info)
 
 Computes the (upper if `uplo = U`, lower if `uplo = L`) pivoted Cholesky
 decomposition of positive-definite matrix `A` with a user-set tolerance
 `tol`, using a preallocated [`CholeskyPivotedWs`](@ref).
+If the workspace was too small and `resize==true` it will be automatically resized.
 `A` is overwritten by its Cholesky decomposition.
 
 Returns `A`, the pivots `piv`, the rank of `A`, and an `info` code. If `info = 0`,
 the factorization succeeded. If `info = i > 0 `, then `A` is indefinite or
 rank-deficient.
 """
-pstrf!(ws::CholeskyPivotedWs, uplo::AbstractChar, A::AbstractMatrix, tol::Real)
+pstrf!(ws::CholeskyPivotedWs, uplo::AbstractChar, A::AbstractMatrix, tol::Real; kwargs...)