Expose algorithm kwargs for Grassmann methods

src/TensorKitManifolds.jl

@@ -28,6 +28,9 @@ checkbase(x, y, z, args...) = checkbase(checkbase(x, y), z, args...)
 # the machine epsilon for the elements of an object X, name inspired from eltype
 scalareps(X) = eps(real(scalartype(X)))
 
+# default SVD algorithm used in the algorithms
+default_svd_alg(::AbstractTensorMap) = TensorKit.SVD()
+
 function isisometry(W::AbstractTensorMap; tol=10 * scalareps(W))
     WdW = W' * W
     s = zero(float(real(scalartype(W))))
@@ -61,7 +64,7 @@ end
 
 struct PolarNewton <: TensorKit.OrthogonalFactorizationAlgorithm
 end
-function projectisometric!(W::AbstractTensorMap; alg=Polar())
+function projectisometric!(W::AbstractTensorMap; alg=default_svd_alg(W))
     if alg isa TensorKit.Polar || alg isa TensorKit.SDD
         foreach(blocks(W)) do (c, b)
             return _polarsdd!(b)
@@ -98,7 +101,7 @@ projecthermitian(W::AbstractTensorMap) = projecthermitian!(copy(W))
 projectantihermitian(W::AbstractTensorMap) = projectantihermitian!(copy(W))
 
 function projectisometric(W::AbstractTensorMap;
-                          alg::TensorKit.OrthogonalFactorizationAlgorithm=Polar())
+                          alg=default_svd_alg(W))
     return projectisometric!(copy(W); alg=alg)
 end
 function projectcomplement(X::AbstractTensorMap, W::AbstractTensorMap,

src/grassmann.jl

@@ -7,7 +7,8 @@ module Grassmann
 using TensorKit
 using TensorKit: similarstoragetype, SectorDict
 using ..TensorKitManifolds: projecthermitian!, projectantihermitian!,
-                            projectisometric!, projectcomplement!, PolarNewton
+                            projectisometric!, projectcomplement!, PolarNewton,
+                            default_svd_alg
 import ..TensorKitManifolds: base, checkbase, inner, retract, transport, transport!
 
 # special type to store tangent vectors using Z
@@ -56,7 +57,7 @@ function Base.getproperty(Δ::GrassmannTangent, sym::Symbol)
     elseif sym ∈ (:U, :S, :V)
         v = Base.getfield(Δ, sym)
         v !== nothing && return v
-        U, S, V, = tsvd(Δ.Z)
+        U, S, V, = tsvd(Δ.Z; alg=default_svd_alg(Δ.Z))
         Base.setfield!(Δ, :U, U)
         Base.setfield!(Δ, :S, S)
         Base.setfield!(Δ, :V, V)
@@ -191,15 +192,15 @@ for the isometries `U`, `V`, and `Y`, and the diagonal matrix `S`, and returning
 `Z = U * S * V` and `Y`.
 """
 function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg=nothing)
-    space(Wold) == space(Wnew) || throw(SectorMismatch())
+    space(Wold) == space(Wnew) || throw(SpaceMismatch())
     WodWn = Wold' * Wnew # V' * cos(S) * V * Y
     Wneworth = Wnew - Wold * WodWn
-    Vd, cS, VY = tsvd!(WodWn)
+    Vd, cS, VY = tsvd!(WodWn; alg=default_svd_alg(WodWn))
     Scmplx = acos(cS)
     # acos always returns a complex TensorMap. We cast back to real if possible.
     S = scalartype(WodWn) <: Real && isreal(sectortype(Scmplx)) ? real(Scmplx) : Scmplx
     UsS = Wneworth * VY' # U * sin(S) # should be in polar decomposition form
-    U = projectisometric!(UsS; alg=Polar())
+    U = projectisometric!(UsS)
     Y = Vd * VY
     V = Vd'
     Z = Grassmann.GrassmannTangent(Wold, U * S * V)
@@ -213,9 +214,9 @@ Return the unitary Y such that V*Y and W are "in the same Grassmann gauge" (tech
 from fibre bundles: in the same section), such that they can be related by a Grassmann
 retraction.
 """
-function relativegauge(W::AbstractTensorMap, V::AbstractTensorMap)
-    space(W) == space(V) || throw(SectorMismatch())
-    return projectisometric!(V' * W; alg=Polar())
+function relativegauge(W::AbstractTensorMap, V::AbstractTensorMap; alg=nothing)
+    space(W) == space(V) || throw(SpaceMismatch())
+    return projectisometric!(V' * W)
 end
 
 function transport!(Θ::GrassmannTangent, W::AbstractTensorMap, Δ::GrassmannTangent, α, W′;