Initial support for CUDA + factorizations (#336)

* Initial support for CUDA + factorizations

* Working adjoint

* Update test/cuda/factorizations.jl

Co-authored-by: Lukas Devos <[email protected]>

* Update test/cuda/factorizations.jl

Co-authored-by: Jutho <[email protected]>

* Fix posdef and copy

* Comments

---------

Co-authored-by: Lukas Devos <[email protected]>
Co-authored-by: Jutho <[email protected]>

Author: 3 people

ext/TensorKitCUDAExt/cutensormap.jl

@@ -119,7 +119,7 @@ function LinearAlgebra.isposdef(t::CuTensorMap)
     InnerProductStyle(spacetype(t)) === EuclideanInnerProduct() || return false
     for (c, b) in blocks(t)
         # do our own hermitian check
-        isherm = MatrixAlgebraKit.ishermitian(b; atol = eps(real(eltype(b))), rtol = eps(real(eltype(b))))
+        isherm = MatrixAlgebraKit.ishermitian(b)
         isherm || return false
         isposdef(Hermitian(b)) || return false
     end

src/factorizations/adjoint.jl

@@ -9,6 +9,9 @@ _adjoint(alg::MAK.LAPACK_HouseholderRQ) = MAK.LAPACK_HouseholderQL(; alg.kwargs.
 _adjoint(alg::MAK.PolarViaSVD) = MAK.PolarViaSVD(_adjoint(alg.svd_alg))
 _adjoint(alg::AbstractAlgorithm) = alg
 
+_adjoint(alg::MAK.CUSOLVER_HouseholderQR) = MAK.LQViaTransposedQR(alg)
+_adjoint(alg::MAK.LQViaTransposedQR) = alg.qr_alg
+
 for f in
     [
         :svd_compact, :svd_full, :svd_vals,
@@ -108,3 +111,7 @@ function MAK.svd_compact!(t::AdjointTensorMap, F, alg::DiagonalAlgorithm)
     F′ = svd_compact!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg))
     return reverse(adjoint.(F′))
 end
+
+function LinearAlgebra.isposdef(t::AdjointTensorMap)
+    return isposdef(adjoint(t))
+end

src/tensors/linalg.jl

@@ -196,14 +196,12 @@ LinearAlgebra.isdiag(t::AbstractTensorMap) = all(LinearAlgebra.isdiag ∘ last,
 
 # In-place methods
 #------------------
-# Wrapping the blocks in a StridedView enables multithreading if JULIA_NUM_THREADS > 1
-# TODO: reconsider this strategy, consider spawning different threads for different blocks
 
 # Copy, adjoint and fill:
 function Base.copy!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap)
     space(tdst) == space(tsrc) || throw(SpaceMismatch("$(space(tdst)) ≠ $(space(tsrc))"))
     for ((c, bdst), (_, bsrc)) in zip(blocks(tdst), blocks(tsrc))
-        copy!(StridedView(bdst), StridedView(bsrc))
+        copy!(bdst, bsrc)
     end
     return tdst
 end