Revert changes to SVD

Author: mtfishman

Project.toml

@@ -22,9 +22,11 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
 
 [weakdeps]
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
+TensorProducts = "decf83d6-1968-43f4-96dc-fdb3fe15fc6d"
 
 [extensions]
 BlockSparseArraysTensorAlgebraExt = "TensorAlgebra"
+BlockSparseArraysTensorProductsExt = "TensorProducts"
 
 [compat]
 Adapt = "4.1.1"
@@ -43,6 +45,7 @@ MatrixAlgebraKit = "0.2.2"
 SparseArraysBase = "0.7.1"
 SplitApplyCombine = "1.2.3"
 TensorAlgebra = "0.3.2"
+TensorProducts = "0.1.7"
 Test = "1.10"
 TypeParameterAccessors = "0.4.1"
 julia = "1.10"

ext/BlockSparseArraysTensorProductsExt/BlockSparseArraysTensorProductsExt.jl

@@ -0,0 +1,14 @@
+module BlockSparseArraysTensorProductsExt
+
+using BlockSparseArrays: BlockUnitRange, blockrange, eachblockaxis
+using TensorProducts: TensorProducts, tensor_product
+# TODO: Dispatch on `FusionStyle` to allow different kinds of products,
+# for example to allow merging common symmetry sectors.
+function TensorProducts.tensor_product(a1::BlockUnitRange, a2::BlockUnitRange)
+  new_blockaxes = vec(
+    map(splat(tensor_product), Iterators.product(eachblockaxis(a1), eachblockaxis(a2)))
+  )
+  return blockrange(new_blockaxes)
+end
+
+end

src/abstractblocksparsearray/abstractblocksparsearray.jl

@@ -93,7 +93,7 @@ function Base.setindex!(
   # https://github.com/JuliaLang/julia/pull/52487).
   # TODO: Delete once we drop support for Julia v1.10.
   aI = @view! a[I...]
-  copyto!(aI, value)
+  aI .= value
   return a
 end
 

src/factorizations/svd.jl

@@ -49,22 +49,21 @@ function MatrixAlgebraKit.initialize_output(
   ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
 )
   bm, bn = blocksize(A)
-  (bmn, mindim) = findmin((bm, bn))
+  bmn = min(bm, bn)
 
   brows = eachblockaxis(axes(A, 1))
   bcols = eachblockaxis(axes(A, 2))
   u_axes = similar(brows, bmn)
-  v_axes = similar(bcols, bmn)
+  v_axes = similar(brows, bmn)
 
   # fill in values for blocks that are present
   bIs = collect(eachblockstoredindex(A))
   browIs = Int.(first.(Tuple.(bIs)))
   bcolIs = Int.(last.(Tuple.(bIs)))
   for bI in eachblockstoredindex(A)
     row, col = Int.(Tuple(bI))
-    dim = (row, col)[mindim]
-    u_axes[dim] = infimum(brows[row], bcols[col])
-    v_axes[dim] = infimum(bcols[col], brows[row])
+    u_axes[col] = infimum(brows[row], bcols[col])
+    v_axes[col] = infimum(bcols[col], brows[row])
   end
 
   # fill in values for blocks that aren't present, pairing them in order of occurence
@@ -84,10 +83,9 @@ function MatrixAlgebraKit.initialize_output(
   # allocate output
   for bI in eachblockstoredindex(A)
     brow, bcol = Tuple(bI)
-    bdim = (brow, bcol)[mindim]
     block = @view!(A[bI])
     block_alg = block_algorithm(alg, block)
-    U[brow, bdim], S[bdim, bdim], Vt[bdim, bcol] = MatrixAlgebraKit.initialize_output(
+    U[brow, bcol], S[bcol, bcol], Vt[bcol, bcol] = MatrixAlgebraKit.initialize_output(
       svd_compact!, block, block_alg
     )
   end