More matrix factorizations (#52)

Author: mtfishman

Project.toml

@@ -1,7 +1,7 @@
 name = "TensorAlgebra"
 uuid = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.2.9"
+version = "0.2.10"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/MatrixAlgebra.jl

@@ -0,0 +1,136 @@
+module MatrixAlgebra
+
+export eigen,
+  eigen!,
+  eigvals,
+  eigvals!,
+  factorize,
+  factorize!,
+  lq,
+  lq!,
+  orth,
+  orth!,
+  polar,
+  polar!,
+  qr,
+  qr!,
+  svd,
+  svd!,
+  svdvals,
+  svdvals!
+
+using LinearAlgebra: LinearAlgebra
+using MatrixAlgebraKit
+
+for (f, f_full, f_compact) in (
+  (:qr, :qr_full, :qr_compact),
+  (:qr!, :qr_full!, :qr_compact!),
+  (:lq, :lq_full, :lq_compact),
+  (:lq!, :lq_full!, :lq_compact!),
+)
+  @eval begin
+    function $f(A::AbstractMatrix; full::Bool=false, kwargs...)
+      f = full ? $f_full : $f_compact
+      return f(A; kwargs...)
+    end
+  end
+end
+
+for (eigen, eigh_full, eig_full, eigh_trunc, eig_trunc) in (
+  (:eigen, :eigh_full, :eig_full, :eigh_trunc, :eig_trunc),
+  (:eigen!, :eigh_full!, :eig_full!, :eigh_trunc!, :eig_trunc!),
+)
+  @eval begin
+    function $eigen(A::AbstractMatrix; trunc=nothing, ishermitian=nothing, kwargs...)
+      ishermitian = @something ishermitian LinearAlgebra.ishermitian(A)
+      f = if !isnothing(trunc)
+        ishermitian ? $eigh_trunc : $eig_trunc
+      else
+        ishermitian ? $eigh_full : $eig_full
+      end
+      return f(A; kwargs...)
+    end
+  end
+end
+
+for (eigvals, eigh_vals, eig_vals) in
+    ((:eigvals, :eigh_vals, :eig_vals), (:eigvals!, :eigh_vals!, :eig_vals!))
+  @eval begin
+    function $eigvals(A::AbstractMatrix; ishermitian=nothing, kwargs...)
+      ishermitian = @something ishermitian LinearAlgebra.ishermitian(A)
+      f = (ishermitian ? $eigh_vals : $eig_vals)
+      return f(A; kwargs...)
+    end
+  end
+end
+
+for (svd, svd_trunc, svd_full, svd_compact) in (
+  (:svd, :svd_trunc, :svd_full, :svd_compact),
+  (:svd!, :svd_trunc!, :svd_full!, :svd_compact!),
+)
+  @eval begin
+    function $svd(A::AbstractMatrix; full::Bool=false, trunc=nothing, kwargs...)
+      return if !isnothing(trunc)
+        @assert !full "Specified both full and truncation, currently not supported"
+        $svd_trunc(A; trunc, kwargs...)
+      else
+        (full ? $svd_full : $svd_compact)(A; kwargs...)
+      end
+    end
+  end
+end
+
+for (svdvals, svd_vals) in ((:svdvals, :svd_vals), (:svdvals!, :svd_vals!))
+  @eval begin
+    function $svdvals(A::AbstractMatrix; ishermitian=nothing, kwargs...)
+      return $svd_vals(A; kwargs...)
+    end
+  end
+end
+
+for (polar, left_polar, right_polar) in
+    ((:polar, :left_polar, :right_polar), (:polar!, :left_polar!, :right_polar!))
+  @eval begin
+    function $polar(A::AbstractMatrix; side=:left, kwargs...)
+      f = if side == :left
+        $left_polar
+      elseif side == :right
+        $right_polar
+      else
+        throw(ArgumentError("`side=$side` not supported."))
+      end
+      return f(A; kwargs...)
+    end
+  end
+end
+
+for (orth, left_orth, right_orth) in
+    ((:orth, :left_orth, :right_orth), (:orth!, :left_orth!, :right_orth!))
+  @eval begin
+    function $orth(A::AbstractMatrix; side=:left, kwargs...)
+      f = if side == :left
+        $left_orth
+      elseif side == :right
+        $right_orth
+      else
+        throw(ArgumentError("`side=$side` not supported."))
+      end
+      return f(A; kwargs...)
+    end
+  end
+end
+
+for (factorize, orth_f) in ((:factorize, :(MatrixAlgebra.orth)), (:factorize!, :orth!))
+  @eval begin
+    function $factorize(A::AbstractMatrix; orth=:left, kwargs...)
+      f = if orth in (:left, :right)
+        $orth_f
+      else
+        throw(ArgumentError("`orth=$orth` not supported."))
+      end
+      return f(A; side=orth, kwargs...)
+    end
+  end
+end
+
+end

src/TensorAlgebra.jl

@@ -1,7 +1,24 @@
 module TensorAlgebra
 
-export contract, contract!, eigen, eigvals, lq, left_null, qr, right_null, svd, svdvals
+export contract,
+  contract!,
+  eigen,
+  eigvals,
+  factorize,
+  left_null,
+  left_orth,
+  left_polar,
+  lq,
+  qr,
+  right_null,
+  right_orth,
+  right_polar,
+  orth,
+  polar,
+  svd,
+  svdvals
 
+include("MatrixAlgebra.jl")
 include("blockedtuple.jl")
 include("blockedpermutation.jl")
 include("BaseExtensions/BaseExtensions.jl")