diff --git a/NDTensors/src/DiagonalArrays/README.md b/NDTensors/src/DiagonalArrays/README.md
index 33996cf52d..1b9355b0a9 100644
--- a/NDTensors/src/DiagonalArrays/README.md
+++ b/NDTensors/src/DiagonalArrays/README.md
@@ -5,15 +5,11 @@ A Julia `DiagonalArray` type.
 ````julia
 using NDTensors.DiagonalArrays:
   DiagonalArray,
-  densearray,
-  diagview,
-  diaglength,
-  getdiagindex,
-  setdiagindex!,
-  setdiag!,
-  diagcopyto!
-
-d = DiagonalArray([1., 2, 3], 3, 4, 5)
+  DiagIndex,
+  DiagIndices,
+  densearray
+
+d = DiagonalArray([1.0, 2, 3], 3, 4, 5)
 @show d[1, 1, 1] == 1
 @show d[2, 2, 2] == 2
 @show d[1, 2, 1] == 0
@@ -21,20 +17,20 @@ d = DiagonalArray([1., 2, 3], 3, 4, 5)
 d[2, 2, 2] = 22
 @show d[2, 2, 2] == 22
 
-@show diaglength(d) == 3
+@show length(d[DiagIndices()]) == 3
 @show densearray(d) == d
-@show getdiagindex(d, 2) == d[2, 2, 2]
+@show d[DiagIndex(2)] == d[2, 2, 2]
 
-setdiagindex!(d, 222, 2)
+d[DiagIndex(2)] = 222
 @show d[2, 2, 2] == 222
 
 a = randn(3, 4, 5)
 new_diag = randn(3)
-setdiag!(a, new_diag)
-diagcopyto!(d, a)
+a[DiagIndices()] = new_diag
+d[DiagIndices()] = a[DiagIndices()]
 
-@show diagview(a) == new_diag
-@show diagview(d) == new_diag
+@show a[DiagIndices()] == new_diag
+@show d[DiagIndices()] == new_diag
 ````
 
 You can generate this README with:
diff --git a/NDTensors/src/DiagonalArrays/examples/README.jl b/NDTensors/src/DiagonalArrays/examples/README.jl
index ede04d6c2e..f5e621385e 100644
--- a/NDTensors/src/DiagonalArrays/examples/README.jl
+++ b/NDTensors/src/DiagonalArrays/examples/README.jl
@@ -2,15 +2,7 @@
 #
 # A Julia `DiagonalArray` type.
 
-using NDTensors.DiagonalArrays:
-  DiagonalArray,
-  densearray,
-  diagview,
-  diaglength,
-  getdiagindex,
-  setdiagindex!,
-  setdiag!,
-  diagcopyto!
+using NDTensors.DiagonalArrays: DiagonalArray, DiagIndex, DiagIndices, densearray
 
 d = DiagonalArray([1.0, 2, 3], 3, 4, 5)
 @show d[1, 1, 1] == 1
@@ -20,20 +12,20 @@ d = DiagonalArray([1.0, 2, 3], 3, 4, 5)
 d[2, 2, 2] = 22
 @show d[2, 2, 2] == 22
 
-@show diaglength(d) == 3
+@show length(d[DiagIndices()]) == 3
 @show densearray(d) == d
-@show getdiagindex(d, 2) == d[2, 2, 2]
+@show d[DiagIndex(2)] == d[2, 2, 2]
 
-setdiagindex!(d, 222, 2)
+d[DiagIndex(2)] = 222
 @show d[2, 2, 2] == 222
 
 a = randn(3, 4, 5)
 new_diag = randn(3)
-setdiag!(a, new_diag)
-diagcopyto!(d, a)
+a[DiagIndices()] = new_diag
+d[DiagIndices()] = a[DiagIndices()]
 
-@show diagview(a) == new_diag
-@show diagview(d) == new_diag
+@show a[DiagIndices()] == new_diag
+@show d[DiagIndices()] == new_diag
 
 # You can generate this README with:
 # ```julia
diff --git a/NDTensors/src/DiagonalArrays/src/DiagonalArrays.jl b/NDTensors/src/DiagonalArrays/src/DiagonalArrays.jl
index 45319ad8ce..c8bb9411d1 100644
--- a/NDTensors/src/DiagonalArrays/src/DiagonalArrays.jl
+++ b/NDTensors/src/DiagonalArrays/src/DiagonalArrays.jl
@@ -3,7 +3,7 @@ module DiagonalArrays
 using Compat # allequal
 using LinearAlgebra
 
-export DiagonalArray
+export DiagonalArray, DiagonalMatrix, DiagonalVector, DiagIndex, DiagIndices, densearray
 
 include("diagview.jl")
 
@@ -19,6 +19,9 @@ struct DiagonalArray{T,N,Diag<:AbstractVector{T},Zero} <: AbstractArray{T,N}
   zero::Zero
 end
 
+const DiagonalVector{T,Diag,Zero} = DiagonalArray{T,1,Diag,Zero}
+const DiagonalMatrix{T,Diag,Zero} = DiagonalArray{T,2,Diag,Zero}
+
 function DiagonalArray{T,N}(
   diag::AbstractVector{T}, d::Tuple{Vararg{Int,N}}, zero=DefaultZero()
 ) where {T,N}
@@ -53,6 +56,25 @@ function DiagonalArray(diag::AbstractVector{T}, d::Vararg{Int,N}) where {T,N}
   return DiagonalArray{T,N}(diag, d)
 end
 
+default_size(diag::AbstractVector, n) = ntuple(Returns(length(diag)), n)
+
+# Infer size from diagonal
+function DiagonalArray{T,N}(diag::AbstractVector) where {T,N}
+  return DiagonalArray{T,N}(diag, default_size(diag, N))
+end
+
+function DiagonalArray{<:Any,N}(diag::AbstractVector{T}) where {T,N}
+  return DiagonalArray{T,N}(diag)
+end
+
+function DiagonalMatrix(diag::AbstractVector)
+  return DiagonalArray{<:Any,2}(diag)
+end
+
+function DiagonalVector(diag::AbstractVector)
+  return DiagonalArray{<:Any,1}(diag)
+end
+
 # undef
 function DiagonalArray{T,N}(::UndefInitializer, d::Tuple{Vararg{Int,N}}) where {T,N}
   return DiagonalArray{T,N}(Vector{T}(undef, minimum(d)), d)
diff --git a/NDTensors/src/DiagonalArrays/src/diagview.jl b/NDTensors/src/DiagonalArrays/src/diagview.jl
index 21f7ac8a06..b3e2547362 100644
--- a/NDTensors/src/DiagonalArrays/src/diagview.jl
+++ b/NDTensors/src/DiagonalArrays/src/diagview.jl
@@ -18,6 +18,19 @@ function setdiagindex!(a::AbstractArray, v, i::Integer)
   return a
 end
 
+struct DiagIndex
+  I::Int
+end
+
+function Base.getindex(a::AbstractArray, i::DiagIndex)
+  return getdiagindex(a, i.I)
+end
+
+function Base.setindex!(a::AbstractArray, v, i::DiagIndex)
+  setdiagindex!(a, v, i.I)
+  return a
+end
+
 function setdiag!(a::AbstractArray, v)
   copyto!(diagview(a), v)
   return a
@@ -28,27 +41,42 @@ function diaglength(a::AbstractArray)
   return minimum(size(a))
 end
 
-function diagstride(A::AbstractArray)
+function diagstride(a::AbstractArray)
   s = 1
   p = 1
-  for i in 1:(ndims(A) - 1)
-    p *= size(A, i)
+  for i in 1:(ndims(a) - 1)
+    p *= size(a, i)
     s += p
   end
   return s
 end
 
-function diagindices(A::AbstractArray)
-  diaglength = minimum(size(A))
-  maxdiag = LinearIndices(A)[CartesianIndex(ntuple(Returns(diaglength), ndims(A)))]
-  return 1:diagstride(A):maxdiag
+function diagindices(a::AbstractArray)
+  diaglength = minimum(size(a))
+  maxdiag = LinearIndices(a)[CartesianIndex(ntuple(Returns(diaglength), ndims(a)))]
+  return 1:diagstride(a):maxdiag
 end
 
-function diagview(A::AbstractArray)
-  return @view A[diagindices(A)]
+function diagindices(a::AbstractArray{<:Any,0})
+  return Base.OneTo(1)
+end
+
+function diagview(a::AbstractArray)
+  return @view a[diagindices(a)]
 end
 
 function diagcopyto!(dest::AbstractArray, src::AbstractArray)
   copyto!(diagview(dest), diagview(src))
   return dest
 end
+
+struct DiagIndices end
+
+function Base.getindex(a::AbstractArray, ::DiagIndices)
+  return diagview(a)
+end
+
+function Base.setindex!(a::AbstractArray, v, ::DiagIndices)
+  setdiag!(a, v)
+  return a
+end
diff --git a/NDTensors/src/NDTensors.jl b/NDTensors/src/NDTensors.jl
index 4fe31337f6..2028b9446d 100644
--- a/NDTensors/src/NDTensors.jl
+++ b/NDTensors/src/NDTensors.jl
@@ -48,6 +48,7 @@ include("exports.jl")
 #####################################
 # General functionality
 #
+include("default_kwargs.jl")
 include("algorithm.jl")
 include("aliasstyle.jl")
 include("abstractarray/set_types.jl")
@@ -148,6 +149,8 @@ include("arraystorage/arraystorage/tensor/eigen.jl")
 include("arraystorage/arraystorage/tensor/svd.jl")
 
 # DiagonalArray storage
+include("arraystorage/diagonalarray/storage/contract.jl")
+
 include("arraystorage/diagonalarray/tensor/contract.jl")
 
 # BlockSparseArray storage
diff --git a/NDTensors/src/arraystorage/arraystorage/storage/arraystorage.jl b/NDTensors/src/arraystorage/arraystorage/storage/arraystorage.jl
index 996494b68a..52f2beeed5 100644
--- a/NDTensors/src/arraystorage/arraystorage/storage/arraystorage.jl
+++ b/NDTensors/src/arraystorage/arraystorage/storage/arraystorage.jl
@@ -6,6 +6,7 @@ const ArrayStorage{T,N} = Union{
   SubArray{T,N},
   PermutedDimsArray{T,N},
   StridedView{T,N},
+  DiagonalArray{T,N},
   BlockSparseArray{T,N},
 }
 
@@ -28,3 +29,8 @@ const MatrixOrArrayStorage{T} = Union{MatrixStorage{T},ArrayStorage{T}}
 function to_arraystorage(x::DenseTensor)
   return tensor(reshape(data(x), size(x)), inds(x))
 end
+
+# TODO: Delete once `Diag` is removed.
+function to_arraystorage(x::DiagTensor)
+  return tensor(DiagonalArray(data(x), size(x)), inds(x))
+end
diff --git a/NDTensors/src/arraystorage/arraystorage/storage/contract.jl b/NDTensors/src/arraystorage/arraystorage/storage/contract.jl
index 0770325f9b..39c9696b48 100644
--- a/NDTensors/src/arraystorage/arraystorage/storage/contract.jl
+++ b/NDTensors/src/arraystorage/arraystorage/storage/contract.jl
@@ -4,10 +4,11 @@ function contract(
   labels1,
   array2::MatrixOrArrayStorage,
   labels2,
-  labelsR=contract_labels(labels1, labels2),
+  labelsR=contract_labels(labels1, labels2);
+  kwargs...,
 )
   output_array = contraction_output(array1, labels1, array2, labels2, labelsR)
-  contract!(output_array, labelsR, array1, labels1, array2, labels2)
+  contract!(output_array, labelsR, array1, labels1, array2, labels2; kwargs...)
   return output_array
 end
 
diff --git a/NDTensors/src/arraystorage/arraystorage/tensor/eigen.jl b/NDTensors/src/arraystorage/arraystorage/tensor/eigen.jl
index d3147dcce6..e55038f7d8 100644
--- a/NDTensors/src/arraystorage/arraystorage/tensor/eigen.jl
+++ b/NDTensors/src/arraystorage/arraystorage/tensor/eigen.jl
@@ -5,10 +5,10 @@
 function eigen(
   T::Hermitian{<:Any,<:ArrayStorageTensor};
   maxdim=nothing,
-  mindim=1,
+  mindim=nothing,
   cutoff=nothing,
-  use_absolute_cutoff=false,
-  use_relative_cutoff=true,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
   # These are getting passed erroneously.
   # TODO: Make sure they don't get passed down
   # to here.
@@ -20,12 +20,6 @@ function eigen(
   ortho=nothing,
   svd_alg=nothing,
 )
-  truncate = !isnothing(maxdim) || !isnothing(cutoff)
-  # TODO: Define `default_maxdim(T)`.
-  maxdim = isnothing(maxdim) ? minimum(dims(T)) : maxdim
-  # TODO: Define `default_cutoff(T)`.
-  cutoff = isnothing(cutoff) ? zero(eltype(T)) : cutoff
-
   matrixT = matrix(T)
   ## TODO Here I am calling parent to ensure that the correct `any` function
   ## is envoked for non-cpu matrices
@@ -45,7 +39,7 @@ function eigen(
   DM = DM[p]
   VM = VM[:, p]
 
-  if truncate
+  if any(!isnothing, (maxdim, cutoff))
     DM, truncerr, _ = truncate!!(
       DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff
     )
@@ -68,7 +62,6 @@ function eigen(
   Vinds = indstype((dag(ind(T, 2)), dag(r)))
   Dinds = indstype((l, dag(r)))
   V = tensor(VM, Vinds)
-  # TODO: Replace with `DiagonalArray`.
-  D = tensor(Diag(DM), Dinds)
+  D = tensor(DiagonalMatrix(DM), Dinds)
   return D, V, spec
 end
diff --git a/NDTensors/src/arraystorage/arraystorage/tensor/qr.jl b/NDTensors/src/arraystorage/arraystorage/tensor/qr.jl
index 04ae4501fc..7ee216341f 100644
--- a/NDTensors/src/arraystorage/arraystorage/tensor/qr.jl
+++ b/NDTensors/src/arraystorage/arraystorage/tensor/qr.jl
@@ -1,4 +1,5 @@
-function qr(A::ArrayStorageTensor)
+function qr(A::ArrayStorageTensor; positive=false)
+  positive && error("Not implemented")
   Q, R = qr(storage(A))
   Q = convert(typeof(R), Q)
   i, j = inds(A)
diff --git a/NDTensors/src/arraystorage/arraystorage/tensor/svd.jl b/NDTensors/src/arraystorage/arraystorage/tensor/svd.jl
index f00543959e..69f446ac53 100644
--- a/NDTensors/src/arraystorage/arraystorage/tensor/svd.jl
+++ b/NDTensors/src/arraystorage/arraystorage/tensor/svd.jl
@@ -1,76 +1,72 @@
-# TODO: Rewrite this function to be more modern:
-# 1. Output `Spectrum` as a keyword argument that gets overwritten.
-# 2. Dispatch on `alg`.
-# 3. Make this into two layers, one that handles indices and one that works with `Matrix`.
+backup_svd_alg(::Algorithm"divide_and_conquer") = Algorithm"qr_iteration"()
+backup_svd_alg(::Algorithm"qr_iteration") = Algorithm"recursive"()
+
+function svd(alg::Algorithm"divide_and_conquer", a::ArrayStorage)
+  USV = svd_catch_error(a; alg=LinearAlgebra.DivideAndConquer())
+  if isnothing(USV)
+    return svd(backup_svd_alg(alg), a)
+  end
+  return USV
+end
+
+function svd(alg::Algorithm"qr_iteration", a::ArrayStorage)
+  USV = svd_catch_error(a; alg=LinearAlgebra.QRIteration())
+  if isnothing(USV)
+    return svd(backup_svd_alg(alg), a)
+  end
+  return USV
+end
+
+function svd(alg::Algorithm"recursive", a::ArrayStorage)
+  return svd_recursive(a)
+end
+
+function svd(::Algorithm"QRAlgorithm", a::ArrayStorage)
+  return error("Not implemented yet")
+end
+
+function svd(::Algorithm"JacobiAlgorithm", a::ArrayStorage)
+  return error("Not implemented yet")
+end
+
+function svd(alg::Algorithm, a::ArrayStorage)
+  return error(
+    "svd algorithm $alg is not currently supported. Please see the documentation for currently supported algorithms.",
+  )
+end
+
 """
-    svd(T::ArrayStorageTensor{<:Number,2}; kwargs...)
+    tsvd(a::ArrayStorage{<:Number,2}; kwargs...)
 
 svd of an order-2 DenseTensor
 """
-function svd(
-  T::ArrayStorageTensor;
+function tsvd(
+  a::ArrayStorage;
+  mindim=nothing,
   maxdim=nothing,
-  mindim=1,
   cutoff=nothing,
-  alg="divide_and_conquer", # TODO: Define `default_alg(T)`
-  use_absolute_cutoff=false,
-  use_relative_cutoff=true,
-  # These are getting passed erroneously.
-  # TODO: Make sure they don't get passed down
-  # to here.
-  which_decomp=nothing,
-  tags=nothing,
-  eigen_perturbation=nothing,
-  normalize=nothing,
+  alg=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
+  # Only used by BlockSparse svd
+  min_blockdim=nothing,
 )
-  truncate = !isnothing(maxdim) || !isnothing(cutoff)
-  # TODO: Define `default_maxdim(T)`.
-  maxdim = isnothing(maxdim) ? minimum(dims(T)) : maxdim
-  # TODO: Define `default_cutoff(T)`.
-  cutoff = isnothing(cutoff) ? zero(eltype(T)) : cutoff
-
-  # TODO: Dispatch on `Algorithm(alg)`.
-  if alg == "divide_and_conquer"
-    MUSV = svd_catch_error(matrix(T); alg=LinearAlgebra.DivideAndConquer())
-    if isnothing(MUSV)
-      # If "divide_and_conquer" fails, try "qr_iteration"
-      alg = "qr_iteration"
-      MUSV = svd_catch_error(matrix(T); alg=LinearAlgebra.QRIteration())
-      if isnothing(MUSV)
-        # If "qr_iteration" fails, try "recursive"
-        alg = "recursive"
-        MUSV = svd_recursive(matrix(T))
-      end
-    end
-  elseif alg == "qr_iteration"
-    MUSV = svd_catch_error(matrix(T); alg=LinearAlgebra.QRIteration())
-    if isnothing(MUSV)
-      # If "qr_iteration" fails, try "recursive"
-      alg = "recursive"
-      MUSV = svd_recursive(matrix(T))
-    end
-  elseif alg == "recursive"
-    MUSV = svd_recursive(matrix(T))
-  elseif alg == "QRAlgorithm" || alg == "JacobiAlgorithm"
-    MUSV = svd_catch_error(matrix(T); alg=alg)
-  else
-    error(
-      "svd algorithm $alg is not currently supported. Please see the documentation for currently supported algorithms.",
-    )
-  end
-  if isnothing(MUSV)
-    if any(isnan, T)
+  alg = replace_nothing(alg, default_svd_alg(a))
+  USV = svd(Algorithm(alg), a)
+  if isnothing(USV)
+    if any(isnan, a)
       println("SVD failed, the matrix you were trying to SVD contains NaNs.")
     else
       println(lapack_svd_error_message(alg))
     end
     return nothing
   end
-  MU, MS, MV = MUSV
-  conj!(MV)
 
-  P = MS .^ 2
-  if truncate
+  U, S, V = USV
+  conj!(V)
+
+  P = S .^ 2
+  if any(!isnothing, (maxdim, cutoff))
     P, truncerr, _ = truncate!!(
       P; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff
     )
@@ -79,25 +75,51 @@ function svd(
   end
   spec = Spectrum(P, truncerr)
   dS = length(P)
-  if dS < length(MS)
-    MU = MU[:, 1:dS]
+  if dS < length(S)
+    U = U[:, 1:dS]
     # Fails on some GPU backends like Metal.
     # resize!(MS, dS)
-    MS = MS[1:dS]
-    MV = MV[:, 1:dS]
+    S = S[1:dS]
+    V = V[:, 1:dS]
   end
+  return U, DiagonalMatrix(S), V, spec
+end
+
+# TODO: Rewrite this function to be more modern:
+# 1. Output `Spectrum` as a keyword argument that gets overwritten.
+# 2. Dispatch on `alg`.
+# 3. Make this into two layers, one that handles indices and one that works with `Matrix`.
+"""
+    svd(T::ArrayStorageTensor{<:Number,2}; kwargs...)
 
+svd of an order-2 DenseTensor
+"""
+function svd(
+  T::ArrayStorageTensor;
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  alg=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
+  # Only used by BlockSparse svd
+  min_blockdim=nothing,
+)
+  U, S, V, spec = tsvd(
+    storage(T); mindim, maxdim, cutoff, alg, use_absolute_cutoff, use_relative_cutoff
+  )
   # Make the new indices to go onto U and V
   # TODO: Put in a separate function, such as
   # `rewrap_inds` or something like that.
+  dS = length(S[DiagIndices()])
   indstype = typeof(inds(T))
   u = eltype(indstype)(dS)
   v = eltype(indstype)(dS)
   Uinds = indstype((ind(T, 1), u))
   Sinds = indstype((u, v))
   Vinds = indstype((ind(T, 2), v))
-  U = tensor(MU, Uinds)
-  S = tensor(Diag(MS), Sinds)
-  V = tensor(MV, Vinds)
-  return U, S, V, spec
+  TU = tensor(U, Uinds)
+  TS = tensor(S, Sinds)
+  TV = tensor(V, Vinds)
+  return TU, TS, TV, spec
 end
diff --git a/NDTensors/src/arraystorage/arraystorage/tensor/zeros.jl b/NDTensors/src/arraystorage/arraystorage/tensor/zeros.jl
index 66759892e4..c990e5db30 100644
--- a/NDTensors/src/arraystorage/arraystorage/tensor/zeros.jl
+++ b/NDTensors/src/arraystorage/arraystorage/tensor/zeros.jl
@@ -11,3 +11,8 @@ end
 function zeros(tensortype::Type{<:ArrayStorageTensor}, inds::Tuple{Vararg{Integer}})
   return _zeros(tensortype, inds)
 end
+
+# To resolve ambiguity error with `Base.zeros`.
+function zeros(tensortype::Type{<:ArrayStorageTensor}, inds::Tuple{})
+  return _zeros(tensortype, inds)
+end
diff --git a/NDTensors/src/arraystorage/diagonalarray/storage/contract.jl b/NDTensors/src/arraystorage/diagonalarray/storage/contract.jl
new file mode 100644
index 0000000000..b86753ac20
--- /dev/null
+++ b/NDTensors/src/arraystorage/diagonalarray/storage/contract.jl
@@ -0,0 +1,172 @@
+# TODO: Move to a different file.
+parenttype(::Type{<:DiagonalArray{<:Any,<:Any,P}}) where {P} = P
+
+# TODO: Move to a different file.
+function promote_rule(
+  storagetype1::Type{<:ArrayStorage}, storagetype2::Type{<:DiagonalArray}
+)
+  # TODO: Replace with `unwrap_type` once
+  # https://github.com/ITensor/ITensors.jl/pull/1220
+  # is merged.
+  return promote_type(storagetype1, leaf_parenttype(storagetype2))
+end
+
+# TODO: Move to a different file.
+function promote_rule(
+  storagetype1::Type{<:DiagonalArray}, storagetype2::Type{<:DiagonalArray}
+)
+  return error("Not implemented yet")
+end
+
+function contraction_output_type(
+  arraytype1::Type{<:DiagonalArray}, arraytype2::Type{<:DiagonalArray}, inds
+)
+  return error("Not implemented yet")
+end
+
+default_convert_to_dense() = true
+
+# TODO: Modernize this function, rewrite in terms of `Array` and `DiagonalArray`.
+# TODO: Move to `storage`.
+function contract!(
+  C::ArrayStorage,
+  Clabels,
+  A::DiagonalArray,
+  Alabels,
+  B::ArrayStorage,
+  Blabels,
+  α::Number=one(eltype(C)),
+  β::Number=zero(eltype(C));
+  convert_to_dense=default_convert_to_dense(),
+)
+  if convert_to_dense
+    contract_dense!(C, Clabels, A, Alabels, B, Blabels, α, β)
+    return C
+  end
+  if !isone(α) || !iszero(β)
+    error(
+      "`contract!(::ArrayStorageTensor, ::DiagTensor, ::ArrayStorageTensor, α, β; convert_to_dense = false)` with `α ≠ 1` or `β ≠ 0` is not currently supported. You can call it with `convert_to_dense = true` instead.",
+    )
+  end
+  astarts = zeros(Int, length(Alabels))
+  bstart = 0
+  cstart = 0
+  b_cstride = 0
+  nbu = 0
+  for ib in 1:length(Blabels)
+    ia = findfirst(==(Blabels[ib]), Alabels)
+    if !isnothing(ia)
+      b_cstride += stride(B, ib)
+      bstart += astarts[ia] * stride(B, ib)
+    else
+      nbu += 1
+    end
+  end
+  c_cstride = 0
+  for ic in 1:length(Clabels)
+    ia = findfirst(==(Clabels[ic]), Alabels)
+    if !isnothing(ia)
+      c_cstride += stride(C, ic)
+      cstart += astarts[ia] * stride(C, ic)
+    end
+  end
+  # strides of the uncontracted dimensions of
+  # B
+  bustride = zeros(Int, nbu)
+  custride = zeros(Int, nbu)
+  # size of the uncontracted dimensions of
+  # B, to be used in CartesianIndices
+  busize = zeros(Int, nbu)
+  n = 1
+  for ib in 1:length(Blabels)
+    if Blabels[ib] > 0
+      bustride[n] = stride(B, ib)
+      busize[n] = size(B, ib)
+      ic = findfirst(==(Blabels[ib]), Clabels)
+      custride[n] = stride(C, ic)
+      n += 1
+    end
+  end
+  boffset_orig = 1 - sum(strides(B))
+  coffset_orig = 1 - sum(strides(C))
+  cartesian_inds = CartesianIndices(Tuple(busize))
+  for inds in cartesian_inds
+    boffset = boffset_orig
+    coffset = coffset_orig
+    for i in 1:nbu
+      ii = inds[i]
+      boffset += ii * bustride[i]
+      coffset += ii * custride[i]
+    end
+    c = zero(eltype(C))
+    for j in 1:DiagonalArrays.diaglength(A)
+      # With α == 0 && β == 1
+      C[cstart + j * c_cstride + coffset] +=
+        A[DiagIndex(j)] * B[bstart + j * b_cstride + boffset]
+      # XXX: not sure if this is correct
+      #C[cstart+j*c_cstride+coffset] += α * A[DiagIndex(j)] * B[bstart+j*b_cstride+boffset] + β * C[cstart+j*c_cstride+coffset]
+    end
+  end
+end
+
+function contract!(
+  C::ArrayStorage{<:Any,0},
+  Clabels,
+  A::DiagonalArray,
+  Alabels,
+  B::ArrayStorage,
+  Blabels,
+  α::Number=one(eltype(C)),
+  β::Number=zero(eltype(C));
+  convert_to_dense=nothing,
+)
+  # If all of B is contracted
+  # TODO: can also check NC+NB==NA
+  min_dim = min(minimum(size(A)), minimum(size(B)))
+  if length(Clabels) == 0
+    # all indices are summed over, just add the product of the diagonal
+    # elements of A and B
+    # Assumes C starts set to 0
+    c₁ = zero(eltype(C))
+    for i in 1:min_dim
+      c₁ += A[DiagIndex(i)] * B[DiagIndex(i)]
+    end
+    C[DiagIndex(1)] = α * c₁ + β * C[DiagIndex(1)]
+  else
+    # not all indices are summed over, set the diagonals of the result
+    # to the product of the diagonals of A and B
+    # TODO: should we make this return a Diag storage?
+    for i in 1:min_dim
+      C[DiagIndex(i)] = α * A[DiagIndex(i)] * B[DiagIndex(i)] + β * C[DiagIndex(i)]
+    end
+  end
+  return C
+end
+
+function contract_dense!(
+  C::ArrayStorage,
+  Clabels,
+  A::DiagonalArray,
+  Alabels,
+  B::ArrayStorage,
+  Blabels,
+  α::Number=one(eltype(C)),
+  β::Number=zero(eltype(C)),
+)
+  return contract!(C, Clabels, densearray(A), Alabels, B, Blabels, α, β)
+end
+
+# Overspecifying types to fix ambiguity error.
+function contract!(
+  C::ArrayStorage,
+  Clabels,
+  A::ArrayStorage,
+  Alabels,
+  B::DiagonalArray,
+  Blabels,
+  α::Number=one(eltype(C)),
+  β::Number=zero(eltype(C));
+  convert_to_dense=default_convert_to_dense(),
+)
+  return contract!(C, Clabels, B, Blabels, A, Alabels, α, β; convert_to_dense)
+end
diff --git a/NDTensors/src/arraystorage/diagonalarray/tensor/contract.jl b/NDTensors/src/arraystorage/diagonalarray/tensor/contract.jl
index e7011896a1..18fe5a1cff 100644
--- a/NDTensors/src/arraystorage/diagonalarray/tensor/contract.jl
+++ b/NDTensors/src/arraystorage/diagonalarray/tensor/contract.jl
@@ -1,142 +1,25 @@
-function promote_rule(storagetype1::Type{<:ArrayStorage}, storagetype2::Type{<:Diag})
-  return promote_type(storagetype1, datatype(storagetype2))
+# The output must be initialized as zero since it is sparse, cannot be undefined
+# Overspecifying types to fix ambiguity error.
+# TODO: Rewrite in terms of `DiagonalArray` and `Array`, not `Tensor`.
+function contraction_output(
+  T1::Tensor{T,N,<:DiagonalArray{T,N,<:AbstractVector{T}}}, T2::ArrayStorageTensor, indsR
+) where {T,N}
+  return zero_contraction_output(T1, T2, indsR)
 end
 
-function contraction_output_type(
-  tensortype1::Type{<:DiagTensor}, tensortype2::Type{<:ArrayStorageTensor}, indsR
-)
-  return similartype(promote_type(tensortype1, tensortype2), indsR)
+# TODO: Rewrite in terms of `DiagonalArray` and `Array`, not `Tensor`.
+function contraction_output(
+  T1::ArrayStorageTensor, T2::Tensor{T,N,<:DiagonalArray{T,N,<:AbstractVector{T}}}, indsR
+) where {T,N}
+  return contraction_output(T2, T1, indsR)
 end
 
-function contraction_output_type(
-  tensortype1::Type{<:ArrayStorageTensor}, tensortype2::Type{<:DiagTensor}, indsR
-)
-  return contraction_output_type(tensortype2, tensortype1, indsR)
-end
-
-# TODO: Modernize this function, rewrite in terms of `Array` and `DiagonalArray`.
-function contract!(
-  C::ArrayStorageTensor{ElC,NC},
-  Clabels,
-  A::DiagTensor{ElA,NA},
-  Alabels,
-  B::ArrayStorageTensor{ElB,NB},
-  Blabels,
-  α::Number=one(ElC),
-  β::Number=zero(ElC);
-  convert_to_dense::Bool=true,
-) where {ElA,NA,ElB,NB,ElC,NC}
-  #@timeit_debug timer "diag-dense contract!" begin
-  if all(i -> i < 0, Blabels)
-    # If all of B is contracted
-    # TODO: can also check NC+NB==NA
-    min_dim = min(minimum(dims(A)), minimum(dims(B)))
-    if length(Clabels) == 0
-      # all indices are summed over, just add the product of the diagonal
-      # elements of A and B
-      # Assumes C starts set to 0
-      c₁ = zero(ElC)
-      for i in 1:min_dim
-        c₁ += getdiagindex(A, i) * getdiagindex(B, i)
-      end
-      setdiagindex!(C, α * c₁ + β * getdiagindex(C, 1), 1)
-    else
-      # not all indices are summed over, set the diagonals of the result
-      # to the product of the diagonals of A and B
-      # TODO: should we make this return a Diag storage?
-      for i in 1:min_dim
-        setdiagindex!(
-          C, α * getdiagindex(A, i) * getdiagindex(B, i) + β * getdiagindex(C, i), i
-        )
-      end
-    end
-  else
-    # Most general contraction
-    if convert_to_dense
-      # TODO: Define a conversion from `DiagonalArray` to `Array`.
-      contract!(C, Clabels, to_arraystorage(dense(A)), Alabels, B, Blabels, α, β)
-    else
-      if !isone(α) || !iszero(β)
-        error(
-          "`contract!(::ArrayStorageTensor, ::DiagTensor, ::ArrayStorageTensor, α, β; convert_to_dense = false)` with `α ≠ 1` or `β ≠ 0` is not currently supported. You can call it with `convert_to_dense = true` instead.",
-        )
-      end
-      astarts = zeros(Int, length(Alabels))
-      bstart = 0
-      cstart = 0
-      b_cstride = 0
-      nbu = 0
-      for ib in 1:length(Blabels)
-        ia = findfirst(==(Blabels[ib]), Alabels)
-        if !isnothing(ia)
-          b_cstride += stride(B, ib)
-          bstart += astarts[ia] * stride(B, ib)
-        else
-          nbu += 1
-        end
-      end
-
-      c_cstride = 0
-      for ic in 1:length(Clabels)
-        ia = findfirst(==(Clabels[ic]), Alabels)
-        if !isnothing(ia)
-          c_cstride += stride(C, ic)
-          cstart += astarts[ia] * stride(C, ic)
-        end
-      end
-
-      # strides of the uncontracted dimensions of
-      # B
-      bustride = zeros(Int, nbu)
-      custride = zeros(Int, nbu)
-      # size of the uncontracted dimensions of
-      # B, to be used in CartesianIndices
-      busize = zeros(Int, nbu)
-      n = 1
-      for ib in 1:length(Blabels)
-        if Blabels[ib] > 0
-          bustride[n] = stride(B, ib)
-          busize[n] = size(B, ib)
-          ic = findfirst(==(Blabels[ib]), Clabels)
-          custride[n] = stride(C, ic)
-          n += 1
-        end
-      end
-
-      boffset_orig = 1 - sum(strides(B))
-      coffset_orig = 1 - sum(strides(C))
-      cartesian_inds = CartesianIndices(Tuple(busize))
-      for inds in cartesian_inds
-        boffset = boffset_orig
-        coffset = coffset_orig
-        for i in 1:nbu
-          ii = inds[i]
-          boffset += ii * bustride[i]
-          coffset += ii * custride[i]
-        end
-        c = zero(ElC)
-        for j in 1:diaglength(A)
-          # With α == 0 && β == 1
-          C[cstart + j * c_cstride + coffset] +=
-            getdiagindex(A, j) * B[bstart + j * b_cstride + boffset]
-          # XXX: not sure if this is correct
-          #C[cstart+j*c_cstride+coffset] += α * getdiagindex(A, j)* B[bstart+j*b_cstride+boffset] + β * C[cstart+j*c_cstride+coffset]
-        end
-      end
-    end
-  end
-  #end # @timeit
-end
-
-function contract!(
-  C::ArrayStorageTensor,
-  Clabels,
-  A::ArrayStorageTensor,
-  Alabels,
-  B::DiagTensor,
-  Blabels,
-  α::Number=one(eltype(C)),
-  β::Number=zero(eltype(C)),
-)
-  return contract!(C, Clabels, B, Blabels, A, Alabels, α, β)
+# Overspecifying types to fix ambiguity error.
+# TODO: Rewrite in terms of `DiagonalArray`, not `Tensor`.
+function contraction_output(
+  tensor1::Tensor{T1,N1,<:DiagonalArray{T1,N1,<:AbstractVector{T1}}},
+  tensor2::Tensor{T2,N2,<:DiagonalArray{T2,N2,<:AbstractVector{T2}}},
+  indsR,
+) where {T1,N1,T2,N2}
+  return zero_contraction_output(tensor1, tensor2, indsR)
 end
diff --git a/NDTensors/src/blocksparse/linearalgebra.jl b/NDTensors/src/blocksparse/linearalgebra.jl
index b9fb8b6471..6d23b1c29d 100644
--- a/NDTensors/src/blocksparse/linearalgebra.jl
+++ b/NDTensors/src/blocksparse/linearalgebra.jl
@@ -4,8 +4,9 @@ const DiagBlockSparseMatrix{ElT,StoreT,IndsT} = DiagBlockSparseTensor{ElT,2,Stor
 const DiagMatrix{ElT,StoreT,IndsT} = DiagTensor{ElT,2,StoreT,IndsT}
 
 function _truncated_blockdim(
-  S::DiagMatrix, docut::Real; singular_values=false, truncate=true, min_blockdim=0
+  S::DiagMatrix, docut::Real; singular_values=false, truncate=true, min_blockdim=nothing
 )
+  min_blockdim = replace_nothing(min_blockdim, 0)
   # TODO: Replace `cpu` with `unwrap_type` dispatch.
   S = cpu(S)
   full_dim = diaglength(S)
@@ -34,13 +35,16 @@ per row/column, otherwise it fails.
 This assumption makes it so the result can be
 computed from the dense svds of seperate blocks.
 """
-function svd(T::Tensor{ElT,2,<:BlockSparse}; kwargs...) where {ElT}
-  alg::String = get(kwargs, :alg, "divide_and_conquer")
-  min_blockdim::Int = get(kwargs, :min_blockdim, 0)
-  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)
-
-  #@timeit_debug timer "block sparse svd" begin
-
+function svd(
+  T::Tensor{ElT,2,<:BlockSparse};
+  min_blockdim=nothing,
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  alg=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
+) where {ElT}
   Us = Vector{DenseTensor{ElT,2}}(undef, nnzblocks(T))
   Ss = Vector{DiagTensor{real(ElT),2}}(undef, nnzblocks(T))
   Vs = Vector{DenseTensor{ElT,2}}(undef, nnzblocks(T))
@@ -50,7 +54,7 @@ function svd(T::Tensor{ElT,2,<:BlockSparse}; kwargs...) where {ElT}
 
   for (n, b) in enumerate(eachnzblock(T))
     blockT = blockview(T, b)
-    USVb = svd(blockT; alg=alg)
+    USVb = svd(blockT; alg)
     if isnothing(USVb)
       return nothing
     end
@@ -77,11 +81,13 @@ function svd(T::Tensor{ElT,2,<:BlockSparse}; kwargs...) where {ElT}
   nzblocksT = nzblocks(T)
 
   dropblocks = Int[]
-  if truncate
-    d, truncerr, docut = truncate!!(d; kwargs...)
+  if any(!isnothing, (maxdim, cutoff))
+    d, truncerr, docut = truncate!!(
+      d; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff
+    )
     for n in 1:nnzblocks(T)
       blockdim = _truncated_blockdim(
-        Ss[n], docut; min_blockdim, singular_values=true, truncate
+        Ss[n], docut; min_blockdim, singular_values=true, truncate=true
       )
       if blockdim == 0
         push!(dropblocks, n)
@@ -197,7 +203,6 @@ function svd(T::Tensor{ElT,2,<:BlockSparse}; kwargs...) where {ElT}
     copyto!(blockview(V, blockV), Vb)
   end
   return U, S, V, Spectrum(d, truncerr)
-  #end # @timeit_debug
 end
 
 _eigen_eltypes(T::Hermitian{ElT,<:BlockSparseMatrix{ElT}}) where {ElT} = real(ElT), ElT
@@ -206,10 +211,13 @@ _eigen_eltypes(T::BlockSparseMatrix{ElT}) where {ElT} = complex(ElT), complex(El
 
 function eigen(
   T::Union{Hermitian{ElT,<:Tensor{ElT,2,<:BlockSparse}},Tensor{ElT,2,<:BlockSparse}};
-  kwargs...,
+  min_blockdim=nothing,
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
 ) where {ElT<:Union{Real,Complex}}
-  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)
-
   ElD, ElV = _eigen_eltypes(T)
 
   # Sorted eigenvalues
@@ -237,10 +245,14 @@ function eigen(
   dropblocks = Int[]
   sort!(d; rev=true, by=abs)
 
-  if truncate
-    d, truncerr, docut = truncate!!(d; kwargs...)
+  if any(!isnothing, (maxdim, cutoff))
+    d, truncerr, docut = truncate!!(
+      d; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff
+    )
     for n in 1:nnzblocks(T)
-      blockdim = _truncated_blockdim(Ds[n], docut)
+      blockdim = _truncated_blockdim(
+        Ds[n], docut; min_blockdim, singular_values=false, truncate=true
+      )
       if blockdim == 0
         push!(dropblocks, n)
       else
@@ -331,7 +343,7 @@ qr(T::BlockSparseTensor{<:Any,2}; kwargs...) = qx(qr, T; kwargs...)
 #  This code thanks to Niklas Tausendpfund 
 #  https://github.com/ntausend/variance_iTensor/blob/main/Hubig_variance_test.ipynb
 #
-function qx(qx::Function, T::BlockSparseTensor{<:Any,2}; kwargs...)
+function qx(qx::Function, T::BlockSparseTensor{<:Any,2}; positive=nothing)
   ElT = eltype(T)
   # getting total number of blocks
   nnzblocksT = nnzblocks(T)
@@ -342,7 +354,7 @@ function qx(qx::Function, T::BlockSparseTensor{<:Any,2}; kwargs...)
 
   for (jj, b) in enumerate(eachnzblock(T))
     blockT = blockview(T, b)
-    QXb = qx(blockT; kwargs...) #call dense qr at src/linearalgebra.jl 387
+    QXb = qx(blockT; positive)
 
     if (isnothing(QXb))
       return nothing
diff --git a/NDTensors/src/default_kwargs.jl b/NDTensors/src/default_kwargs.jl
new file mode 100644
index 0000000000..73d8756d84
--- /dev/null
+++ b/NDTensors/src/default_kwargs.jl
@@ -0,0 +1,9 @@
+replace_nothing(::Nothing, replacement) = replacement
+replace_nothing(value, replacement) = value
+
+default_maxdim(a) = minimum(size(a))
+default_mindim(a) = true
+default_cutoff(a) = zero(eltype(a))
+default_svd_alg(a) = "divide_and_conquer"
+default_use_absolute_cutoff(a) = false
+default_use_relative_cutoff(a) = true
diff --git a/NDTensors/src/linearalgebra/linearalgebra.jl b/NDTensors/src/linearalgebra/linearalgebra.jl
index 53f339e3d7..ce5e7e2cb0 100644
--- a/NDTensors/src/linearalgebra/linearalgebra.jl
+++ b/NDTensors/src/linearalgebra/linearalgebra.jl
@@ -94,38 +94,18 @@ end
 
 svd of an order-2 DenseTensor
 """
-function svd(T::DenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
-  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)
-
-  #
-  # Keyword argument deprecations
-  #
-  use_absolute_cutoff = false
-  if haskey(kwargs, :absoluteCutoff)
-    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
-    use_absolute_cutoff = get(kwargs, :absoluteCutoff, use_absolute_cutoff)
-  end
-
-  use_relative_cutoff = true
-  if haskey(kwargs, :doRelCutoff)
-    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
-    use_relative_cutoff = get(kwargs, :doRelCutoff, use_relative_cutoff)
-  end
-
-  if haskey(kwargs, :fastsvd) || haskey(kwargs, :fastSVD)
-    error(
-      "In svd, fastsvd/fastSVD keyword arguments are removed in favor of alg, see documentation for more details.",
-    )
-  end
-
-  maxdim::Int = get(kwargs, :maxdim, minimum(dims(T)))
-  mindim::Int = get(kwargs, :mindim, 1)
-  cutoff = get(kwargs, :cutoff, 0.0)
-  use_absolute_cutoff::Bool = get(kwargs, :use_absolute_cutoff, use_absolute_cutoff)
-  use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff)
-  alg::String = get(kwargs, :alg, "divide_and_conquer")
-
-  #@timeit_debug timer "dense svd" begin
+function svd(
+  T::DenseTensor{ElT,2,IndsT};
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
+  alg=nothing,
+  # Only used by BlockSparse svd
+  min_blockdim=nothing,
+) where {ElT,IndsT}
+  alg = replace_nothing(alg, default_svd_alg(T))
   if alg == "divide_and_conquer"
     MUSV = svd_catch_error(matrix(T); alg=LinearAlgebra.DivideAndConquer())
     if isnothing(MUSV)
@@ -148,7 +128,7 @@ function svd(T::DenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
   elseif alg == "recursive"
     MUSV = svd_recursive(matrix(T))
   elseif alg == "QRAlgorithm" || alg == "JacobiAlgorithm"
-    MUSV = svd_catch_error(matrix(T); alg=alg)
+    MUSV = svd_catch_error(matrix(T); alg)
   else
     error(
       "svd algorithm $alg is not currently supported. Please see the documentation for currently supported algorithms.",
@@ -167,9 +147,9 @@ function svd(T::DenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
   #end # @timeit_debug
 
   P = MS .^ 2
-  if truncate
+  if any(!isnothing, (maxdim, cutoff))
     P, truncerr, _ = truncate!!(
-      P; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff, kwargs...
+      P; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff
     )
   else
     truncerr = 0.0
@@ -197,27 +177,13 @@ function svd(T::DenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
 end
 
 function eigen(
-  T::Hermitian{ElT,<:DenseTensor{ElT,2,IndsT}}; kwargs...
+  T::Hermitian{ElT,<:DenseTensor{ElT,2,IndsT}};
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
 ) where {ElT<:Union{Real,Complex},IndsT}
-  # Keyword argument deprecations
-  use_absolute_cutoff = false
-  if haskey(kwargs, :absoluteCutoff)
-    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
-    use_absolute_cutoff = get(kwargs, :absoluteCutoff, use_absolute_cutoff)
-  end
-  use_relative_cutoff = true
-  if haskey(kwargs, :doRelCutoff)
-    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
-    use_relative_cutoff = get(kwargs, :doRelCutoff, use_relative_cutoff)
-  end
-
-  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)
-  maxdim::Int = get(kwargs, :maxdim, minimum(dims(T)))
-  mindim::Int = get(kwargs, :mindim, 1)
-  cutoff::Union{Nothing,Float64} = get(kwargs, :cutoff, 0.0)
-  use_absolute_cutoff::Bool = get(kwargs, :use_absolute_cutoff, use_absolute_cutoff)
-  use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff)
-
   matrixT = matrix(T)
   ## TODO Here I am calling parent to ensure that the correct `any` function
   ## is envoked for non-cpu matrices
@@ -237,9 +203,9 @@ function eigen(
   DM = DM[p]
   VM = VM[:, p]
 
-  if truncate
+  if any(!isnothing, (maxdim, cutoff))
     DM, truncerr, _ = truncate!!(
-      DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff, kwargs...
+      DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff
     )
     dD = length(DM)
     if dD < size(VM, 2)
@@ -317,27 +283,13 @@ random_orthog(::Type{ElT}, n::Int, m::Int) where {ElT<:Real} = random_unitary(El
 random_orthog(n::Int, m::Int) = random_orthog(Float64, n, m)
 
 function eigen(
-  T::DenseTensor{ElT,2,IndsT}; kwargs...
+  T::DenseTensor{ElT,2,IndsT};
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
 ) where {ElT<:Union{Real,Complex},IndsT}
-  # Keyword argument deprecations
-  use_absolute_cutoff = false
-  if haskey(kwargs, :absoluteCutoff)
-    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
-    use_absolute_cutoff = get(kwargs, :absoluteCutoff, use_absolute_cutoff)
-  end
-  use_relative_cutoff = true
-  if haskey(kwargs, :doRelCutoff)
-    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
-    use_relative_cutoff = get(kwargs, :doRelCutoff, use_relative_cutoff)
-  end
-
-  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)
-  maxdim::Int = get(kwargs, :maxdim, minimum(dims(T)))
-  mindim::Int = get(kwargs, :mindim, 1)
-  cutoff::Float64 = get(kwargs, :cutoff, 0.0)
-  use_absolute_cutoff::Bool = get(kwargs, :use_absolute_cutoff, use_absolute_cutoff)
-  use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff)
-
   matrixT = matrix(T)
   if any(!isfinite, matrixT)
     throw(
@@ -354,9 +306,9 @@ function eigen(
   #DM = DM[p]
   #VM = VM[:,p]
 
-  if truncate
+  if any(!isnothing, (maxdim, cutoff))
     DM, truncerr, _ = truncate!!(
-      DM; maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff, kwargs...
+      DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff
     )
     dD = length(DM)
     if dD < size(VM, 2)
@@ -380,23 +332,23 @@ function eigen(
   return D, V, spec
 end
 
-# NDTensors.Unwrap.qr
-function qr(T::DenseTensor{<:Any,2}; positive=false, kwargs...)
+# LinearAlgebra.qr
+function qr(T::DenseTensor{<:Any,2}; positive=false)
   qxf = positive ? qr_positive : qr
-  return qx(qxf, T; kwargs...)
+  return qx(qxf, T)
 end
 
 # NDTensors.Unwrap.ql
-function ql(T::DenseTensor{<:Any,2}; positive=false, kwargs...)
+function ql(T::DenseTensor{<:Any,2}; positive=false)
   qxf = positive ? ql_positive : ql
-  return qx(qxf, T; kwargs...)
+  return qx(qxf, T)
 end
 
 #
 #  Generic function for qr and ql decomposition of dense matrix.
 #  The X tensor = R or L.
 #
-function qx(qx::Function, T::DenseTensor{<:Any,2}; kwargs...)
+function qx(qx::Function, T::DenseTensor{<:Any,2})
   QM, XM = qx(expose(matrix(T)))
   # Be aware that if positive==false, then typeof(QM)=LinearAlgebra.QRCompactWYQ, not Matrix
   # It gets converted to matrix below.
@@ -485,7 +437,7 @@ end
 #  Lapack replaces A with Q & R carefully packed together.  So here we just copy a
 #  before letting lapack overwirte it. 
 #
-function ql(A::AbstractMatrix; kwargs...)
+function ql(A::AbstractMatrix)
   Base.require_one_based_indexing(A)
   T = eltype(A)
   AA = similar(A, LinearAlgebra._qreltype(T), size(A))
@@ -495,7 +447,7 @@ function ql(A::AbstractMatrix; kwargs...)
     cutype = unwrap_type(AA)
     AA = NDTensors.cpu(AA)
   end
-  Q, L = ql!(AA; kwargs...)
+  Q, L = ql!(AA)
   if iscuda
     Q = adapt(cutype, Q)
     L = adapt(cutype, L)
diff --git a/NDTensors/src/truncate.jl b/NDTensors/src/truncate.jl
index de691f786e..beba20d937 100644
--- a/NDTensors/src/truncate.jl
+++ b/NDTensors/src/truncate.jl
@@ -18,31 +18,22 @@ function truncate!!(::Type{<:AbstractArray}, P::AbstractArray; kwargs...)
 end
 
 # CPU implementation.
-function truncate!(P::AbstractVector{ElT}; cutoff=zero(eltype(P)), kwargs...) where {ElT}
-  if isnothing(cutoff)
-    cutoff = typemin(ElT)
-  end
-
-  # Keyword argument deprecations
-  use_absolute_cutoff = false
-  if haskey(kwargs, :absoluteCutoff)
-    @warn "In truncate!, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
-    use_absolute_cutoff = get(kwargs, :absoluteCutoff, use_absolute_cutoff)
-  end
-  use_relative_cutoff = true
-  if haskey(kwargs, :doRelCutoff)
-    @warn "In truncate!, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
-    use_relative_cutoff = get(kwargs, :doRelCutoff, use_relative_cutoff)
-  end
-
-  maxdim::Int = min(get(kwargs, :maxdim, length(P)), length(P))
-  mindim::Int = max(get(kwargs, :mindim, 1), 1)
-
-  use_absolute_cutoff::Bool = get(kwargs, :use_absolute_cutoff, use_absolute_cutoff)
-  use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff)
+function truncate!(
+  P::AbstractVector;
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
+)
+  mindim = replace_nothing(mindim, default_mindim(P))
+  maxdim = replace_nothing(maxdim, length(P))
+  cutoff = replace_nothing(cutoff, typemin(eltype(P)))
+  use_absolute_cutoff = replace_nothing(use_absolute_cutoff, default_use_absolute_cutoff(P))
+  use_relative_cutoff = replace_nothing(use_relative_cutoff, default_use_relative_cutoff(P))
 
   origm = length(P)
-  docut = zero(ElT)
+  docut = zero(eltype(P))
 
   #if P[1] <= 0.0
   #  P[1] = 0.0
@@ -52,7 +43,7 @@ function truncate!(P::AbstractVector{ElT}; cutoff=zero(eltype(P)), kwargs...) wh
 
   if origm == 1
     docut = abs(P[1]) / 2
-    return zero(ElT), docut
+    return zero(eltype(P)), docut
   end
 
   s = sign(P[1])
@@ -60,12 +51,12 @@ function truncate!(P::AbstractVector{ElT}; cutoff=zero(eltype(P)), kwargs...) wh
 
   #Zero out any negative weight
   for n in origm:-1:1
-    (P[n] >= zero(ElT)) && break
-    P[n] = zero(ElT)
+    (P[n] >= zero(eltype(P))) && break
+    P[n] = zero(eltype(P))
   end
 
   n = origm
-  truncerr = zero(ElT)
+  truncerr = zero(eltype(P))
   while n > maxdim
     truncerr += P[n]
     n -= 1
@@ -79,10 +70,10 @@ function truncate!(P::AbstractVector{ElT}; cutoff=zero(eltype(P)), kwargs...) wh
       n -= 1
     end
   else
-    scale = one(ElT)
+    scale = one(eltype(P))
     if use_relative_cutoff
       scale = sum(P)
-      (scale == zero(ElT)) && (scale = one(ElT))
+      (scale == zero(eltype(P))) && (scale = one(eltype(P)))
     end
 
     #Continue truncating until *sum* of discarded probability 
@@ -101,13 +92,12 @@ function truncate!(P::AbstractVector{ElT}; cutoff=zero(eltype(P)), kwargs...) wh
 
   if n < origm
     docut = (P[n] + P[n + 1]) / 2
-    if abs(P[n] - P[n + 1]) < ElT(1e-3) * P[n]
-      docut += ElT(1e-3) * P[n]
+    if abs(P[n] - P[n + 1]) < eltype(P)(1e-3) * P[n]
+      docut += eltype(P)(1e-3) * P[n]
     end
   end
 
   s < 0 && (P .*= s)
   resize!(P, n)
-
   return truncerr, docut
 end
diff --git a/NDTensors/test/arraytensor/diagonalarray.jl b/NDTensors/test/arraytensor/diagonalarray.jl
new file mode 100644
index 0000000000..70c0029dfc
--- /dev/null
+++ b/NDTensors/test/arraytensor/diagonalarray.jl
@@ -0,0 +1,26 @@
+using NDTensors
+using NDTensors.DiagonalArrays
+using LinearAlgebra
+using Test
+
+using NDTensors: storage, storagetype
+
+@testset "Tensor wrapping DiagonalArray" begin
+  D = DiagonalArray(randn(3), 3, 4, 5)
+  Dᵈ = densearray(D)
+  A = randn(3, 4, 5)
+
+  for convert_to_dense in (true, false)
+    @test contract(D, (-1, -2, -3), A, (-1, -2, -3); convert_to_dense) ≈
+      contract(Dᵈ, (-1, -2, -3), A, (-1, -2, -3))
+    @test contract(D, (-1, -2, 1), A, (-1, -2, 2); convert_to_dense) ≈
+      contract(Dᵈ, (-1, -2, 1), A, (-1, -2, 2))
+  end
+
+  # Tensor tests
+  Dᵗ = tensor(D, size(D))
+  Dᵈᵗ = tensor(Dᵈ, size(D))
+  Aᵗ = tensor(A, size(A))
+  @test contract(Dᵗ, (-1, -2, -3), Aᵗ, (-1, -2, -3)) ≈
+    contract(Dᵈᵗ, (-1, -2, -3), Aᵗ, (-1, -2, -3))
+end
diff --git a/src/mps/abstractmps.jl b/src/mps/abstractmps.jl
index a6309e7d2f..94722e301f 100644
--- a/src/mps/abstractmps.jl
+++ b/src/mps/abstractmps.jl
@@ -1592,7 +1592,8 @@ bond indices is performed. Afterward, tensors
 Either modify in-place with `orthogonalize!` or
 out-of-place with `orthogonalize`.
 """
-function orthogonalize!(M::AbstractMPS, j::Int; kwargs...)
+function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothing)
+  # TODO: Delete `maxdim` and `normalize` keyword arguments.
   @debug_check begin
     if !(1 <= j <= length(M))
       error("Input j=$j to `orthogonalize!` out of range (valid range = 1:$(length(M)))")
@@ -1608,7 +1609,7 @@ function orthogonalize!(M::AbstractMPS, j::Int; kwargs...)
     else
       ltags = TagSet("Link,l=$b")
     end
-    L, R = factorize(M[b], linds; tags=ltags, kwargs...)
+    L, R = factorize(M[b], linds; tags=ltags, maxdim)
     M[b] = L
     M[b + 1] *= R
     setleftlim!(M, b)
@@ -1629,7 +1630,7 @@ function orthogonalize!(M::AbstractMPS, j::Int; kwargs...)
     else
       ltags = TagSet("Link,l=$b")
     end
-    L, R = factorize(M[b + 1], rinds; tags=ltags, kwargs...)
+    L, R = factorize(M[b + 1], rinds; tags=ltags, maxdim)
     M[b + 1] = L
     M[b] *= R
 
diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl
index 5fecdd0abf..cd85d488f6 100644
--- a/src/mps/dmrg.jl
+++ b/src/mps/dmrg.jl
@@ -269,7 +269,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
             phi,
             1,
             eigsolve_which_eigenvalue;
-            ishermitian=ishermitian,
+            ishermitian,
             tol=eigsolve_tol,
             krylovdim=eigsolve_krylovdim,
             maxiter=eigsolve_maxiter,
@@ -311,10 +311,10 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
             mindim=mindim(sweeps, sw),
             cutoff=cutoff(sweeps, sw),
             eigen_perturbation=drho,
-            ortho=ortho,
+            ortho,
             normalize=true,
-            which_decomp=which_decomp,
-            svd_alg=svd_alg,
+            which_decomp,
+            svd_alg,
           )
         end
         maxtruncerr = max(maxtruncerr, spec.truncerr)
diff --git a/src/mps/mpo.jl b/src/mps/mpo.jl
index 84550e5558..c3b47fb5c0 100644
--- a/src/mps/mpo.jl
+++ b/src/mps/mpo.jl
@@ -597,13 +597,12 @@ end
 
 Apply(A::MPO, ψ::MPS; kwargs...) = Applied(apply, (A, ψ), NamedTuple(kwargs))
 
-function contract(A::MPO, ψ::MPS; alg="densitymatrix", kwargs...)
-  if haskey(kwargs, :method)
-    # Backwards compatibility, use `method`.
-    alg = get(kwargs, :method, "densitymatrix")
-  end
+function contract(A::MPO, ψ::MPS; alg=nothing, method=alg, kwargs...)
+  # TODO: Delete `method` since it is deprecated.
+  alg = NDTensors.replace_nothing(method, "densitymatrix")
 
   # Keyword argument deprecations
+  # TODO: Delete these.
   if alg == "DensityMatrix"
     @warn "In contract, method DensityMatrix is deprecated in favor of densitymatrix"
     alg = "densitymatrix"
diff --git a/src/mps/mps.jl b/src/mps/mps.jl
index 2a19b6c1e0..8d3f721687 100644
--- a/src/mps/mps.jl
+++ b/src/mps/mps.jl
@@ -527,19 +527,24 @@ Factorize the ITensor `phi` and replace the ITensors
 `b` and `b+1` of MPS `M` with the factors. Choose
 the orthogonality with `ortho="left"/"right"`.
 """
-function replacebond!(M::MPS, b::Int, phi::ITensor; kwargs...)
-  ortho::String = get(kwargs, :ortho, "left")
-  swapsites::Bool = get(kwargs, :swapsites, false)
-  which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
-  normalize::Bool = get(kwargs, :normalize, false)
-
-  # Deprecated keywords
-  if haskey(kwargs, :dir)
-    error(
-      """dir keyword in replacebond! has been replaced by ortho.
-      Note that the options are now the same as factorize, so use `left` instead of `fromleft` and `right` instead of `fromright`.""",
-    )
-  end
+function replacebond!(
+  M::MPS,
+  b::Int,
+  phi::ITensor;
+  normalize=nothing,
+  swapsites=nothing,
+  ortho=nothing,
+  # Decomposition kwargs
+  which_decomp=nothing,
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  eigen_perturbation=nothing,
+  svd_alg=nothing,
+)
+  normalize = NDTensors.replace_nothing(normalize, false)
+  swapsites = NDTensors.replace_nothing(swapsites, false)
+  ortho = NDTensors.replace_nothing(ortho, "left")
 
   indsMb = inds(M[b])
   if swapsites
@@ -548,7 +553,16 @@ function replacebond!(M::MPS, b::Int, phi::ITensor; kwargs...)
     indsMb = replaceind(indsMb, sb, sbp1)
   end
   L, R, spec = factorize(
-    phi, indsMb; which_decomp=which_decomp, tags=tags(linkind(M, b)), kwargs...
+    phi,
+    indsMb;
+    mindim,
+    maxdim,
+    cutoff,
+    ortho,
+    which_decomp,
+    eigen_perturbation,
+    svd_alg,
+    tags=tags(linkind(M, b)),
   )
   M[b] = L
   M[b + 1] = R
diff --git a/src/tensor_operations/matrix_decomposition.jl b/src/tensor_operations/matrix_decomposition.jl
index 9da68a2c00..378bdfa826 100644
--- a/src/tensor_operations/matrix_decomposition.jl
+++ b/src/tensor_operations/matrix_decomposition.jl
@@ -106,14 +106,30 @@ Utrunc2, Strunc2, Vtrunc2 = svd(A, i, k; cutoff=1e-10);
 
 See also: [`factorize`](@ref), [`eigen`](@ref)
 """
-function svd(A::ITensor, Linds...; leftdir=nothing, rightdir=nothing, kwargs...)
-  utags::TagSet = get(kwargs, :lefttags, get(kwargs, :utags, "Link,u"))
-  vtags::TagSet = get(kwargs, :righttags, get(kwargs, :vtags, "Link,v"))
+function svd(
+  A::ITensor,
+  Linds...;
+  leftdir=nothing,
+  rightdir=nothing,
+  lefttags=nothing,
+  righttags=nothing,
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  alg=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
+  min_blockdim=nothing,
+  # Deprecated
+  utags=lefttags,
+  vtags=righttags,
+)
+  lefttags = NDTensors.replace_nothing(lefttags, ts"Link,u")
+  righttags = NDTensors.replace_nothing(righttags, ts"Link,v")
 
-  # Keyword argument deprecations
-  #if haskey(kwargs, :utags) || haskey(kwargs, :vtags)
-  #  @warn "Keyword arguments `utags` and `vtags` are deprecated in favor of `leftags` and `righttags`."
-  #end
+  # Deprecated
+  utags = NDTensors.replace_nothing(utags, ts"Link,u")
+  vtags = NDTensors.replace_nothing(vtags, ts"Link,v")
 
   Lis = commoninds(A, indices(Linds...))
   Ris = uniqueinds(A, Lis)
@@ -142,7 +158,16 @@ function svd(A::ITensor, Linds...; leftdir=nothing, rightdir=nothing, kwargs...)
     AC = permute(AC, cL, cR)
   end
 
-  USVT = svd(tensor(AC); kwargs...)
+  USVT = svd(
+    tensor(AC);
+    mindim,
+    maxdim,
+    cutoff,
+    alg,
+    use_absolute_cutoff,
+    use_relative_cutoff,
+    min_blockdim,
+  )
   if isnothing(USVT)
     return nothing
   end
@@ -155,10 +180,10 @@ function svd(A::ITensor, Linds...; leftdir=nothing, rightdir=nothing, kwargs...)
   U = UC * dag(CL)
   V = VC * dag(CR)
 
-  settags!(U, utags, u)
-  settags!(S, utags, u)
-  settags!(S, vtags, v)
-  settags!(V, vtags, v)
+  U = settags(U, utags, u)
+  S = settags(S, utags, u)
+  S = settags(S, vtags, v)
+  V = settags(V, vtags, v)
 
   u = settags(u, utags)
   v = settags(v, vtags)
@@ -262,7 +287,31 @@ A * U ≈ Ul * D # true
 
 See also: [`svd`](@ref), [`factorize`](@ref)
 """
-function eigen(A::ITensor, Linds, Rinds; kwargs...)
+function eigen(
+  A::ITensor,
+  Linds,
+  Rinds;
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
+  ishermitian=nothing,
+  tags=nothing,
+  lefttags=nothing,
+  righttags=nothing,
+  plev=nothing,
+  leftplev=nothing,
+  rightplev=nothing,
+)
+  ishermitian = NDTensors.replace_nothing(ishermitian, false)
+  tags = NDTensors.replace_nothing(tags, ts"Link,eigen")
+  lefttags = NDTensors.replace_nothing(lefttags, tags)
+  righttags = NDTensors.replace_nothing(righttags, tags)
+  plev = NDTensors.replace_nothing(plev, 0)
+  leftplev = NDTensors.replace_nothing(leftplev, plev)
+  rightplev = NDTensors.replace_nothing(rightplev, plev)
+
   @debug_check begin
     if hasqns(A)
       @assert flux(A) == QN()
@@ -276,16 +325,6 @@ function eigen(A::ITensor, Linds, Rinds; kwargs...)
   N != NL + NR &&
     error("Number of left and right indices must add up to total number of indices")
 
-  ishermitian::Bool = get(kwargs, :ishermitian, false)
-
-  tags::TagSet = get(kwargs, :tags, "Link,eigen")
-  lefttags::TagSet = get(kwargs, :lefttags, tags)
-  righttags::TagSet = get(kwargs, :righttags, tags)
-
-  plev::Int = get(kwargs, :plev, 0)
-  leftplev::Int = get(kwargs, :leftplev, plev)
-  rightplev::Int = get(kwargs, :rightplev, plev)
-
   if lefttags == righttags && leftplev == rightplev
     leftplev = rightplev + 1
   end
@@ -335,7 +374,7 @@ function eigen(A::ITensor, Linds, Rinds; kwargs...)
 
   AT = ishermitian ? Hermitian(tensor(AC)) : tensor(AC)
 
-  DT, VT, spec = eigen(AT; kwargs...)
+  DT, VT, spec = eigen(AT; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff)
   D, VC = itensor(DT), itensor(VT)
 
   V = VC * dag(CR)
@@ -445,7 +484,9 @@ end
 #
 #  Generic function implementing both qr and ql decomposition. The X tensor = R or L. 
 #
-function qx(qx::Function, A::ITensor, Linds::Indices, Rinds::Indices; tags, kwargs...)
+function qx(
+  qx::Function, A::ITensor, Linds::Indices, Rinds::Indices; tags=ts"Link,qx", positive=false
+)
   # Strip out any extra indices that are not in A.
   # Unit test test/base/test_itensor.jl line 1469 will fail without this.
   Linds = commoninds(A, Linds)
@@ -468,7 +509,7 @@ function qx(qx::Function, A::ITensor, Linds::Indices, Rinds::Indices; tags, kwar
   #
   AC = permute(AC, cL, cR; allow_alias=true)
 
-  QT, XT = qx(tensor(AC); kwargs...) #pass order(AC)==2 matrix down to the NDTensors level where qr/ql are implemented.
+  QT, XT = qx(tensor(AC); positive) #pass order(AC)==2 matrix down to the NDTensors level where qr/ql are implemented.
   #
   #  Undo the combine oepration, to recover all tensor indices.
   #
@@ -491,8 +532,10 @@ end
 #  Generic function implementing both rq and lq decomposition. Implemented using qr/ql 
 #  with swapping the left and right indices.  The X tensor = R or L. 
 #
-function xq(qx::Function, A::ITensor, Linds::Indices, Rinds::Indices; tags, kwargs...)
-  Q, X, q = qx(A, Rinds, Linds; kwargs...)
+function xq(
+  qx::Function, A::ITensor, Linds::Indices, Rinds::Indices; tags=ts"Link,xq", positive=false
+)
+  Q, X, q = qx(A, Rinds, Linds; positive)
   #
   # fix up the tag name for the index between Q and L.
   #  
@@ -507,8 +550,8 @@ polar(A::ITensor; kwargs...) = error(noinds_error_message("polar"))
 
 # TODO: allow custom tags in internal indices?
 # TODO: return the new common indices?
-function polar(A::ITensor, Linds...; kwargs...)
-  U, S, V = svd(A, Linds...; kwargs...)
+function polar(A::ITensor, Linds...)
+  U, S, V = svd(A, Linds...)
   u = commoninds(S, U)
   v = commoninds(S, V)
   δᵤᵥ′ = δ(u..., v'...)
@@ -517,13 +560,12 @@ function polar(A::ITensor, Linds...; kwargs...)
   return Q, P, commoninds(Q, P)
 end
 
-function factorize_qr(A::ITensor, Linds...; kwargs...)
-  ortho::String = get(kwargs, :ortho, "left")
+function factorize_qr(A::ITensor, Linds...; ortho="left", tags=nothing, positive=false)
   if ortho == "left"
-    L, R, q = qr(A, Linds...; kwargs...)
+    L, R, q = qr(A, Linds...; tags, positive)
   elseif ortho == "right"
     Lis = uniqueinds(A, indices(Linds...))
-    R, L, q = qr(A, Lis...; kwargs...)
+    R, L, q = qr(A, Lis...; tags, positive)
   else
     error("In factorize using qr decomposition, ortho keyword
     $ortho not supported. Supported options are left or right.")
@@ -551,12 +593,27 @@ function factorize_svd(
   Linds...;
   (singular_values!)=nothing,
   ortho="left",
-  svd_alg="divide_and_conquer",
-  dir=ITensors.In,
-  kwargs...,
+  alg=nothing,
+  dir=nothing,
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  tags=nothing,
 )
+  dir = NDTensors.replace_nothing(dir, ITensors.In)
   leftdir, rightdir = -dir, -dir
-  USV = svd(A, Linds...; leftdir, rightdir, alg=svd_alg, kwargs...)
+  USV = svd(
+    A,
+    Linds...;
+    leftdir,
+    rightdir,
+    alg,
+    mindim,
+    maxdim,
+    cutoff,
+    lefttags=tags,
+    righttags=tags,
+  )
   if isnothing(USV)
     return nothing
   end
@@ -579,9 +636,17 @@ function factorize_svd(
   return L, R, spec
 end
 
-function factorize_eigen(A::ITensor, Linds...; kwargs...)
-  ortho::String = get(kwargs, :ortho, "left")
-  delta_A2 = get(kwargs, :eigen_perturbation, nothing)
+function factorize_eigen(
+  A::ITensor,
+  Linds...;
+  ortho="left",
+  eigen_perturbation=nothing,
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  tags=nothing,
+)
+  delta_A2 = eigen_perturbation
   if ortho == "left"
     Lis = commoninds(A, indices(Linds...))
   elseif ortho == "right"
@@ -599,7 +664,7 @@ function factorize_eigen(A::ITensor, Linds...; kwargs...)
     noprime!(delta_A2)
     A2 += delta_A2
   end
-  F = eigen(A2, Lis, simLis; ishermitian=true, kwargs...)
+  F = eigen(A2, Lis, simLis; ishermitian=true, mindim, maxdim, cutoff, tags)
   D, _, spec = F
   L = F.Vt
   R = dag(L) * A
@@ -654,13 +719,26 @@ Perform a factorization of `A` into ITensors `L` and `R` such that `A ≈ L * R`
 
 For truncation arguments, see: [`svd`](@ref)
 """
-function factorize(A::ITensor, Linds...; maxdim=nothing, kwargs...)
-  ortho::String = get(kwargs, :ortho, "left")
-  tags::TagSet = get(kwargs, :tags, "Link,fact")
-  plev::Int = get(kwargs, :plev, 0)
-  which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
-  cutoff = get(kwargs, :cutoff, nothing)
-  eigen_perturbation = get(kwargs, :eigen_perturbation, nothing)
+function factorize(
+  A::ITensor,
+  Linds...;
+  mindim=nothing,
+  maxdim=nothing,
+  cutoff=nothing,
+  ortho=nothing,
+  tags=nothing,
+  plev=nothing,
+  which_decomp=nothing,
+  # eigen
+  eigen_perturbation=nothing,
+  # svd
+  svd_alg=nothing,
+  use_absolute_cutoff=nothing,
+  use_relative_cutoff=nothing,
+  min_blockdim=nothing,
+  (singular_values!)=nothing,
+  dir=nothing,
+)
   if !isnothing(eigen_perturbation)
     if !(isnothing(which_decomp) || which_decomp == "eigen")
       error("""when passing a non-trivial eigen_perturbation to `factorize`,
@@ -669,11 +747,9 @@ function factorize(A::ITensor, Linds...; maxdim=nothing, kwargs...)
     end
     which_decomp = "eigen"
   end
-
-  if haskey(kwargs, :which_factorization)
-    error("""which_factorization keyword in factorize has
-    been replace by which_decomp.""")
-  end
+  ortho = NDTensors.replace_nothing(ortho, "left")
+  tags = NDTensors.replace_nothing(tags, ts"Link,fact")
+  plev = NDTensors.replace_nothing(plev, 0)
 
   # Determines when to use eigen vs. svd (eigen is less precise,
   # so eigen should only be used if a larger cutoff is requested)
@@ -697,17 +773,18 @@ function factorize(A::ITensor, Linds...; maxdim=nothing, kwargs...)
       which_decomp = "eigen"
     end
   end
-
   if which_decomp == "svd"
-    LR = factorize_svd(A, Linds...; kwargs..., maxdim=maxdim)
+    LR = factorize_svd(
+      A, Linds...; mindim, maxdim, cutoff, tags, ortho, alg=svd_alg, dir, singular_values!
+    )
     if isnothing(LR)
       return nothing
     end
     L, R, spec = LR
   elseif which_decomp == "eigen"
-    L, R, spec = factorize_eigen(A, Linds...; kwargs..., maxdim=maxdim)
+    L, R, spec = factorize_eigen(A, Linds...; mindim, maxdim, cutoff, tags, ortho)
   elseif which_decomp == "qr"
-    L, R = factorize_qr(A, Linds...; kwargs...)
+    L, R = factorize_qr(A, Linds...; ortho, tags)
     spec = Spectrum(nothing, 0.0)
   else
     throw(ArgumentError("""In factorize, factorization $which_decomp is not
@@ -717,8 +794,8 @@ function factorize(A::ITensor, Linds...; maxdim=nothing, kwargs...)
   # Set the tags and prime level
   l = commonind(L, R)
   l̃ = setprime(settags(l, tags), plev)
-  replaceind!(L, l, l̃)
-  replaceind!(R, l, l̃)
+  L = replaceind(L, l, l̃)
+  R = replaceind(R, l, l̃)
   l = l̃
 
   return L, R, spec, l