CompatHelper: bump compat for MatrixAlgebraKit to 0.6, (keep existing compat)

Project.toml

@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
+version = "0.10.11"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.10.10"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -41,7 +41,7 @@ GPUArraysCore = "0.1, 0.2"
 LinearAlgebra = "1.10"
 MacroTools = "0.5.13"
 MapBroadcast = "0.1.5"
-MatrixAlgebraKit = "0.5"
+MatrixAlgebraKit = "0.6"
 SparseArraysBase = "0.7.1"
 SplitApplyCombine = "1.2.3"
 TensorAlgebra = "0.3, 0.4"

src/factorizations/orthnull.jl

@@ -1,26 +1,21 @@
-using MatrixAlgebraKit:
-    MatrixAlgebraKit,
-    default_svd_algorithm,
-    left_null!,
-    left_null_svd!,
-    left_orth!,
-    left_polar!,
-    lq_compact!,
-    null_truncation_strategy,
-    qr_compact!,
-    right_null!,
-    right_null_svd!,
-    right_orth!,
-    right_polar!,
-    select_algorithm,
-    svd_compact!
+using MatrixAlgebraKit: MatrixAlgebraKit, AbstractAlgorithm, default_svd_algorithm,
+    left_null!, left_orth!, left_polar!, lq_compact!, null_truncation_strategy, qr_compact!,
+    right_null!, right_orth!, right_polar!, select_algorithm, svd_compact!
 
-function MatrixAlgebraKit.initialize_output(
-        ::typeof(left_orth!), A::AbstractBlockSparseMatrix
-    )
-    return nothing
+using MatrixAlgebraKit: LeftOrthAlgorithm
+for kind in ("polar", "qr", "svd")
+    @eval begin
+        function MatrixAlgebraKit.initialize_output(
+                ::typeof(left_orth!), A::AbstractBlockSparseMatrix,
+                alg::LeftOrthAlgorithm{Symbol($kind)},
+            )
+            return nothing
+        end
+    end
 end
-function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::AbstractBlockSparseMatrix, F)
+function MatrixAlgebraKit.check_input(
+        ::typeof(left_orth!), A::AbstractBlockSparseMatrix, F, alg::AbstractAlgorithm
+    )
     !isnothing(F) && throw(
         ArgumentError(
             "`left_orth!` on block sparse matrices does not support specifying the output"
@@ -29,55 +24,60 @@ function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::AbstractBlockSpar
     return nothing
 end
 
-function MatrixAlgebraKit.left_orth_qr!(A::AbstractBlockSparseMatrix, F, alg)
+using MatrixAlgebraKit: LeftOrthViaQR
+function MatrixAlgebraKit.left_orth!(A::AbstractBlockSparseMatrix, F, alg::LeftOrthViaQR)
     !isnothing(F) && throw(
         ArgumentError(
             "`left_orth!` on block sparse matrices does not support specifying the output"
         ),
     )
-    alg′ = select_algorithm(qr_compact!, A, alg)
+    alg′ = select_algorithm(qr_compact!, A, alg.alg)
     return qr_compact!(A, alg′)
 end
-function MatrixAlgebraKit.left_orth_polar!(A::AbstractBlockSparseMatrix, F, alg)
+using MatrixAlgebraKit: LeftOrthViaPolar
+function MatrixAlgebraKit.left_orth!(A::AbstractBlockSparseMatrix, F, alg::LeftOrthViaPolar)
     !isnothing(F) && throw(
         ArgumentError(
             "`left_orth!` on block sparse matrices does not support specifying the output"
         ),
     )
-    alg′ = select_algorithm(left_polar!, A, alg)
+    alg′ = select_algorithm(left_polar!, A, alg.alg)
     return left_polar!(A, alg′)
 end
-function MatrixAlgebraKit.left_orth_svd!(
-        A::AbstractBlockSparseMatrix, F, alg, trunc::Nothing = nothing
+using MatrixAlgebraKit: LeftOrthViaSVD, does_truncate
+function MatrixAlgebraKit.left_orth!(
+        A::AbstractBlockSparseMatrix, F, alg::LeftOrthViaSVD
     )
     !isnothing(F) && throw(
         ArgumentError(
             "`left_orth!` on block sparse matrices does not support specifying the output"
         ),
     )
-    alg′ = select_algorithm(svd_compact!, A, alg)
-    U, S, Vᴴ = svd_compact!(A, alg′)
-    return U, S * Vᴴ
-end
-function MatrixAlgebraKit.left_orth_svd!(A::AbstractBlockSparseMatrix, F, alg, trunc)
-    !isnothing(F) && throw(
-        ArgumentError(
-            "`left_orth!` on block sparse matrices does not support specifying the output"
-        ),
-    )
-    alg′ = select_algorithm(svd_compact!, A, alg)
-    alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc)
-    U, S, Vᴴ = svd_trunc!(A, alg_trunc)
+    U, S, Vᴴ = if !does_truncate(alg.alg)
+        alg′ = select_algorithm(svd_compact!, A, alg.alg)
+        svd_compact!(A, alg′)
+    else
+        alg′ = select_algorithm(svd_compact!, A, alg.alg)
+        alg_trunc = select_algorithm(svd_trunc!, A, alg′)
+        svd_trunc!(A, alg_trunc)
+    end
     return U, S * Vᴴ
 end
 
-function MatrixAlgebraKit.initialize_output(
-        ::typeof(right_orth!), A::AbstractBlockSparseMatrix
-    )
-    return nothing
+using MatrixAlgebraKit: RightOrthAlgorithm
+for kind in ("lq", "polar", "svd")
+    @eval begin
+        function MatrixAlgebraKit.initialize_output(
+                ::typeof(right_orth!), A::AbstractBlockSparseMatrix,
+                alg::RightOrthAlgorithm{Symbol($kind)},
+            )
+            return nothing
+        end
+    end
 end
 function MatrixAlgebraKit.check_input(
-        ::typeof(right_orth!), A::AbstractBlockSparseMatrix, F::Nothing
+        ::typeof(right_orth!), A::AbstractBlockSparseMatrix, F::Nothing,
+        alg::AbstractAlgorithm,
     )
     !isnothing(F) && throw(
         ArgumentError(
@@ -87,65 +87,63 @@ function MatrixAlgebraKit.check_input(
     return nothing
 end
 
-function MatrixAlgebraKit.right_orth_lq!(A::AbstractBlockSparseMatrix, F, alg)
+using MatrixAlgebraKit: RightOrthViaLQ
+function MatrixAlgebraKit.right_orth!(A::AbstractBlockSparseMatrix, F, alg::RightOrthViaLQ)
     !isnothing(F) && throw(
         ArgumentError(
             "`right_orth!` on block sparse matrices does not support specifying the output"
         ),
     )
-    alg′ = select_algorithm(lq_compact!, A, alg)
+    alg′ = select_algorithm(lq_compact!, A, alg.alg)
     return lq_compact!(A, alg′)
 end
-function MatrixAlgebraKit.right_orth_polar!(A::AbstractBlockSparseMatrix, F, alg)
+using MatrixAlgebraKit: RightOrthViaPolar
+function MatrixAlgebraKit.right_orth!(
+        A::AbstractBlockSparseMatrix, F, alg::RightOrthViaPolar
+    )
     !isnothing(F) && throw(
         ArgumentError(
             "`right_orth!` on block sparse matrices does not support specifying the output"
         ),
     )
-    alg′ = select_algorithm(right_polar!, A, alg)
+    alg′ = select_algorithm(right_polar!, A, alg.alg)
     return right_polar!(A, alg′)
 end
-function MatrixAlgebraKit.right_orth_svd!(
-        A::AbstractBlockSparseMatrix, F, alg, trunc::Nothing = nothing
+using MatrixAlgebraKit: RightOrthViaSVD
+function MatrixAlgebraKit.right_orth!(
+        A::AbstractBlockSparseMatrix, F, alg::RightOrthViaSVD
     )
     !isnothing(F) && throw(
         ArgumentError(
             "`right_orth!` on block sparse matrices does not support specifying the output"
         ),
     )
-    alg′ = select_algorithm(svd_compact!, A, alg)
-    U, S, Vᴴ = svd_compact!(A, alg′)
-    return U * S, Vᴴ
-end
-function MatrixAlgebraKit.right_orth_svd!(A::AbstractBlockSparseMatrix, F, alg, trunc)
-    !isnothing(F) && throw(
-        ArgumentError(
-            "`right_orth!` on block sparse matrices does not support specifying the output"
-        ),
-    )
-    alg′ = select_algorithm(svd_compact!, A, alg)
-    alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc)
-    U, S, Vᴴ = svd_trunc!(A, alg_trunc)
+    U, S, Vᴴ = if !does_truncate(alg.alg)
+        alg′ = select_algorithm(svd_compact!, A, alg.alg)
+        svd_compact!(A, alg′)
+    else
+        alg′ = select_algorithm(svd_compact!, A, alg.alg)
+        alg_trunc = select_algorithm(svd_trunc!, A, alg′)
+        svd_trunc!(A, alg_trunc)
+    end
     return U * S, Vᴴ
 end
 
 function MatrixAlgebraKit.initialize_output(
-        ::typeof(left_null!), A::AbstractBlockSparseMatrix
+        ::typeof(left_null!), A::AbstractBlockSparseMatrix, alg::AbstractAlgorithm
     )
     return nothing
 end
 function MatrixAlgebraKit.check_input(
-        ::typeof(left_null!), A::AbstractBlockSparseMatrix, N::Nothing
+        ::typeof(left_null!), A::AbstractBlockSparseMatrix, N::Nothing,
+        alg::AbstractAlgorithm,
     )
     return nothing
 end
-function MatrixAlgebraKit.left_null_qr!(A::AbstractBlockSparseMatrix, N, alg)
-    return left_null_svd!(A, N, default_svd_algorithm(A))
-end
-function MatrixAlgebraKit.left_null_svd!(
-        A::AbstractBlockSparseMatrix, N, alg, trunc::Nothing
+function MatrixAlgebraKit.left_null!(
+        A::AbstractBlockSparseMatrix, N, alg::AbstractAlgorithm
     )
-    return left_null_svd!(A, N, alg, null_truncation_strategy(; atol = 0, rtol = 0))
+    return error("Not implemented.")
 end
 function MatrixAlgebraKit.truncate(
         ::typeof(left_null!),
@@ -156,16 +154,19 @@ function MatrixAlgebraKit.truncate(
 end
 
 function MatrixAlgebraKit.initialize_output(
-        ::typeof(right_null!), A::AbstractBlockSparseMatrix
+        ::typeof(right_null!), A::AbstractBlockSparseMatrix, alg::AbstractAlgorithm
     )
     return nothing
 end
 function MatrixAlgebraKit.check_input(
-        ::typeof(right_null!), A::AbstractBlockSparseMatrix, N::Nothing
+        ::typeof(right_null!), A::AbstractBlockSparseMatrix, N::Nothing,
+        alg::AbstractAlgorithm,
     )
     return nothing
 end
-function MatrixAlgebraKit.right_null_lq!(A::AbstractBlockSparseMatrix, N, alg)
+function MatrixAlgebraKit.right_null!(
+        A::AbstractBlockSparseMatrix, N, alg::AbstractAlgorithm
+    )
     return error("Not implement.")
 end
 function MatrixAlgebraKit.truncate(

src/factorizations/polar.jl

@@ -27,7 +27,7 @@ end
 function MatrixAlgebraKit.left_polar!(A::AbstractBlockSparseMatrix, alg::PolarViaSVD)
     check_input(left_polar!, A)
     # TODO: Use more in-place operations here, avoid `copy`.
-    U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+    U, S, Vᴴ = svd_compact!(A, alg.svd_alg)
     W = U * Vᴴ
     # TODO: `copy` is required for now because of:
     # https://github.com/ITensor/BlockSparseArrays.jl/issues/24
@@ -38,7 +38,7 @@ end
 function MatrixAlgebraKit.right_polar!(A::AbstractBlockSparseMatrix, alg::PolarViaSVD)
     check_input(right_polar!, A)
     # TODO: Use more in-place operations here, avoid `copy`.
-    U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+    U, S, Vᴴ = svd_compact!(A, alg.svd_alg)
     Wᴴ = U * Vᴴ
     # TODO: `copy` is required for now because of:
     # https://github.com/ITensor/BlockSparseArrays.jl/issues/24

src/factorizations/svd.jl

@@ -1,6 +1,6 @@
 using DiagonalArrays: diagonaltype
-using MatrixAlgebraKit:
-    MatrixAlgebraKit, check_input, default_svd_algorithm, svd_compact!, svd_full!, svd_vals!
+using MatrixAlgebraKit: MatrixAlgebraKit, check_input, default_svd_algorithm, svd_compact!,
+    svd_full!, svd_vals!
 using TypeParameterAccessors: realtype
 
 function MatrixAlgebraKit.default_svd_algorithm(
@@ -239,3 +239,14 @@ function MatrixAlgebraKit.svd_vals!(
     end
     return S
 end
+
+# Computing truncation errors. Currently not supported for block sparse matrices,
+# needs to be implemented.
+struct BlockSparseDiagView{T, A <: AbstractBlockSparseMatrix{T}} <:
+    AbstractBlockSparseVector{T}
+    a::A
+end
+MatrixAlgebraKit.diagview(a::AbstractBlockSparseMatrix) = BlockSparseDiagView(a)
+# TODO: Implement this, and/or develop `BlockSparseDiagView` so that the generic
+# version in MatrixAlgebraKit works.
+MatrixAlgebraKit.truncation_error!(a::BlockSparseDiagView, ind) = nothing

src/factorizations/utility.jl

@@ -1,5 +1,22 @@
+using LinearAlgebra: LinearAlgebra
 using MatrixAlgebraKit: MatrixAlgebraKit
 
+# Friendlier definitions for block sparse matrices. MatrixAlgebraKit.jl v0.6
+# defines them in terms of `diagview`:
+# https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/blob/v0.6.0/src/common/matrixproperties.jl#L40-L64
+# which isn't well supported right now for block sparse matrices.
+# TODO: Improve `diagview` for block sparse matrices and remove these definitions.
+function MatrixAlgebraKit.is_left_isometric(
+        A::AnyAbstractBlockSparseMatrix; isapprox_kwargs...
+    )
+    return isapprox(A' * A, LinearAlgebra.I; isapprox_kwargs...)
+end
+function MatrixAlgebraKit.is_right_isometric(
+        A::AnyAbstractBlockSparseMatrix; isapprox_kwargs...
+    )
+    return MatrixAlgebraKit.is_left_isometric(A'; isapprox_kwargs...)
+end
+
 function infimum(r1::AbstractUnitRange, r2::AbstractUnitRange)
     (isone(first(r1)) && isone(first(r2))) ||
         throw(ArgumentError("infimum only defined for ranges starting at 1"))

test/Project.toml

@@ -31,7 +31,7 @@ DiagonalArrays = "0.3"
 GPUArraysCore = "0.2"
 JLArrays = "0.2, 0.3"
 LinearAlgebra = "1"
-MatrixAlgebraKit = "0.5"
+MatrixAlgebraKit = "0.6"
 Random = "1"
 SafeTestsets = "0.1"
 SparseArraysBase = "0.7"