Skip to content

Commit

Permalink
[ITensorGPU] [ENHANCEMENT] Improve functionality for transferring dat…
Browse files Browse the repository at this point in the history
…a between CPU and GPU (#956)

* Fix `expect` and `correlation_matrix` for MPS on GPU.

* Make sure QR and SVD preserve element type (i.e. single precision) in more cases.
  • Loading branch information
mtfishman authored Aug 4, 2022
1 parent 0a8f086 commit d7d6ec8
Show file tree
Hide file tree
Showing 48 changed files with 803 additions and 533 deletions.
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,9 @@ ITensorMakie/Manifest.toml
ITensorMakie/test/Manifest.toml
ITensorGLMakie/Manifest.toml
ITensorGLMakie/test/Manifest.toml
ITensorGPU/test/Manifest.toml
ITensorVisualizationBase/Manifest.toml
ITensorVisualizationBase/test/Manifest.toml
NDTensors/Manifest.toml
/Manifest.toml
docs/Manifest.toml
3 changes: 3 additions & 0 deletions ITensorGPU/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -7,12 +7,14 @@ version = "0.0.5"
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
GPUCompiler = "61eb1bfa-7361-4325-ad38-22787b887f55"
HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
Expand All @@ -21,6 +23,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
Adapt = "3.3.1"
CUDA = "3.5.0"
Combinatorics = "1.0.2"
Functors = "0.3.0"
GPUArrays = "8.1.2"
GPUCompiler = "0.13.8"
HDF5 = "0.15.7"
Expand Down
12 changes: 10 additions & 2 deletions ITensorGPU/src/ITensorGPU.jl
Original file line number Diff line number Diff line change
@@ -1,19 +1,22 @@
module ITensorGPU

using CUDA
using CUDA.Adapt
using CUDA.CUTENSOR
using CUDA.CUBLAS
using CUDA.CUSOLVER
using Functors
using LinearAlgebra
using Random, Strided
using TimerOutputs
using SimpleTraits
using StaticArrays
using ITensors
using ITensors.NDTensors
using Strided
import CUDA: CuArray, CuMatrix, CuVector, cu
import CUDA.CUTENSOR: cutensorContractionPlan_t, cutensorAlgo_t

import CUDA.Adapt: adapt_structure
import CUDA.Mem: pin
#=
const devs = Ref{Vector{CUDAdrv.CuDevice}}()
Expand Down Expand Up @@ -48,6 +51,7 @@ import ITensors:
BroadcastStyle,
Indices
import ITensors.NDTensors:
can_contract,
similartype,
ContractionProperties,
contract!!,
Expand Down Expand Up @@ -82,8 +86,13 @@ import ITensors.NDTensors:
Ctrans,
_contract_scalar!,
_contract_scalar_noperm!

using ITensors.NDTensors: setdata, setstorage, cpu, IsWrappedArray, parenttype

import Base.*, Base.permutedims!
import Base: similar
include("traits.jl")
include("adapt.jl")
include("tensor/cudense.jl")
include("tensor/dense.jl")
include("tensor/culinearalgebra.jl")
Expand All @@ -92,7 +101,6 @@ include("tensor/cucombiner.jl")
include("tensor/cudiag.jl")
include("cuitensor.jl")
include("mps/cumps.jl")
include("mps/cumpo.jl")

#const ContractionPlans = Dict{String, Tuple{cutensorAlgo_t, cutensorContractionPlan_t}}()
const ContractionPlans = Dict{String,cutensorAlgo_t}()
Expand Down
25 changes: 25 additions & 0 deletions ITensorGPU/src/adapt.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
#
# Used to adapt `EmptyStorage` types
#

NDTensors.cu(eltype::Type{<:Number}, x) = fmap(x -> adapt(CuArray{eltype}, x), x)
NDTensors.cu(x) = fmap(x -> adapt(CuArray, x), x)

NDTensors.to_vector_type(arraytype::Type{CuArray}) = CuVector
NDTensors.to_vector_type(arraytype::Type{CuArray{T}}) where {T} = CuVector{T}

function NDTensors.set_eltype_if_unspecified(
arraytype::Type{CuVector{T}}, eltype::Type
) where {T}
return arraytype
end
function NDTensors.set_eltype_if_unspecified(arraytype::Type{CuVector}, eltype::Type)
return CuVector{eltype}
end

# Overload `CUDA.cu` for convenience
const ITensorType = Union{
TensorStorage,Tensor,ITensor,Array{ITensor},Array{<:Array{ITensor}},MPS,MPO
}
CUDA.cu(x::ITensorType) = NDTensors.cu(x)
CUDA.cu(eltype::Type{<:Number}, x::ITensorType) = NDTensors.cu(eltype, x)
47 changes: 17 additions & 30 deletions ITensorGPU/src/cuitensor.jl
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,20 @@ function ITensor(
data = CuArray{eltype}(as, A)
return itensor(Dense(data), inds)
end

# Fix ambiguity error
function ITensor(
as::NDTensors.AliasStyle, eltype::Type{<:Number}, A::CuArray{<:Number}, inds::Tuple{}
)
length(A) dim(inds) && throw(
DimensionMismatch(
"In ITensor(::CuArray, inds), length of AbstractArray ($(length(A))) must match total dimension of IndexSet ($(dim(inds)))",
),
)
data = CuArray{eltype}(as, A)
return itensor(Dense(data), inds)
end

# Helper functions for different view behaviors
CuArray{ElT,N}(::NeverAlias, A::AbstractArray) where {ElT,N} = CuArray{ElT,N}(A)
function CuArray{ElT,N}(::AllowAlias, A::AbstractArray) where {ElT,N}
Expand All @@ -50,38 +64,11 @@ function CuArray{<:Any,N}(as::AliasStyle, A::AbstractArray{ElTA,N}) where {N,ElT
return CuArray{ElTA,N}(as, A)
end

#TODO: check that the size of the Array matches the Index dimensions
function cuITensor(A::Array{S}, inds) where {S<:Number}
return ITensor(Dense(CuArray{S}(A)), inds)
end
function cuITensor(A::CuArray{S}, inds::IndexSet) where {S<:Number}
return ITensor(Dense(A), inds)
end
cuITensor(A::Array{S}, inds::Index...) where {S<:Number} = cuITensor(A, IndexSet(inds...))
cuITensor(A::CuArray{S}, inds::Index...) where {S<:Number} = cuITensor(A, IndexSet(inds...))

function cuITensor(A::ITensor)
return if storage(tensor(A)) isa ITensors.EmptyStorage
cuITensor(zero(eltype(storage(tensor(A)))), inds(A)...)
else
cuITensor(data(tensor(A)), inds(A)...)
end
end

cu(A::ITensor) = cuITensor(A)
cuITensor(data::Array, inds...) = cu(ITensor(data, inds...))

# Helpful for moving gate structures to GPU
cu(A::Array{ITensor}) = map(cu, A)
cu(A::Array{<:Array{ITensor}}) = map(cu, A)

function cpu(A::ITensor)
typeof(data(storage(A))) <: CuArray && return ITensor(cpu(storage(A)), inds(A))
return A
end
cuITensor(data::CuArray, inds...) = ITensor(data, inds...)

# Helpful for moving gate structures to CPU
cpu(A::Array{ITensor}) = map(cpu, A)
cpu(A::Array{<:Array{ITensor}}) = map(cpu, A)
cuITensor(A::ITensor) = cu(A)

function randomCuITensor(::Type{S}, inds::Indices) where {S<:Real}
T = cuITensor(S, inds)
Expand Down
31 changes: 0 additions & 31 deletions ITensorGPU/src/mps/cumpo.jl

This file was deleted.

143 changes: 10 additions & 133 deletions ITensorGPU/src/mps/cumps.jl
Original file line number Diff line number Diff line change
@@ -1,136 +1,13 @@
function cuMPS(psi::MPS)
phi = copy(psi)
for site in 1:length(psi)
phi.data[site] = cuITensor(psi.data[site])
end
return phi
end
# cu(ψ::Union{MPS,MPO}) = map(cu, ψ)
# cpu(ψ::Union{MPS,MPO}) = map(cpu, ψ)

cu::MPS) = cuMPS(ψ)
cuMPS::MPS) = cu(ψ)
cuMPS(args...; kwargs...) = cu(MPS(args...; kwargs...))
randomCuMPS(args...; kwargs...) = cu(randomMPS(args...; kwargs...))

cuMPS() = MPS()
function cuMPS(::Type{T}, sites::Vector{<:Index}; linkdims::Integer=1) where {T<:Number}
N = length(sites)
v = Vector{ITensor}(undef, N)
if N == 1
v[1] = emptyITensor(T, sites[1])
return MPS(v)
end
# For backwards compatibility
productCuMPS(args...; kwargs...) = cuMPS(args...; kwargs...)

space = if hasqns(sites)
[QN() => linkdims]
else
linkdims
end

l = [Index(1, "Link,l=$ii") for ii in 1:(N - 1)]
for ii in eachindex(sites)
s = sites[ii]
if ii == 1
v[ii] = cuITensor(l[ii], s)
elseif ii == N
v[ii] = cuITensor(l[ii - 1], s)
else
v[ii] = cuITensor(l[ii - 1], s, l[ii])
end
end
return MPS(v, 0, N + 1)
end
function cuMPS(sites::Vector{<:Index}, args...; kwargs...)
return cuMPS(Float64, sites, args...; kwargs...)
end
function cuMPS(N::Int; ortho_lims::UnitRange=1:N)
return MPS(Vector{ITensor}(undef, N); ortho_lims=ortho_lims)
end

function randomCuMPS(sites)
M = cuMPS(sites)
for i in eachindex(sites)
randn!(M[i])
normalize!(M[i])
end
M.llim = 1
M.rlim = length(M)
return M
end
function randomCuMPS(N::Int; ortho_lims::UnitRange=1:N)
return randomCuMPS(Vector{ITensor}(undef, N); ortho_lims=ortho_lims)
end

const productCuMPS = cuMPS

function cuMPS(::Type{T}, ivals::Vector{<:Pair{<:Index}}) where {T<:Number}
N = length(ivals)
As = Vector{ITensor}(undef, N)
links = Vector{Index}(undef, N)
for n in 1:N
s = ITensors.ind(ivals[n])
links[n] = Index(1, "Link,l=$n")
if n == 1
A = ITensor(T, s, links[n])
A[ivals[n], links[n](1)] = 1.0
elseif n == N
A = ITensor(T, links[n - 1], s)
A[links[n - 1](1), ivals[n]] = 1.0
else
A = ITensor(T, links[n - 1], s, links[n])
A[links[n - 1](1), ivals[n], links[n](1)] = 1.0
end
As[n] = cuITensor(A)
end
return MPS(As, 0, 2)
end
cuMPS(ivals::Vector{Pair{<:Index}}) = cuMPS(Float64, ivals)

function cuMPS(::Type{T}, sites::Vector{<:Index}, states_) where {T<:Number}
if length(sites) != length(states_)
throw(DimensionMismatch("Number of sites and and initial vals don't match"))
end
N = length(states_)
M = cuMPS(N)

if N == 1
M[1] = state(sites[1], states_[1])
return M
end

states = [state(sites[j], states_[j]) for j in 1:N]

if hasqns(states[1])
lflux = QN()
for j in 1:(N - 1)
lflux += flux(states[j])
end
links = Vector{QNIndex}(undef, N - 1)
for j in (N - 1):-1:1
links[j] = Index(lflux => 1; tags="Link,l=$j", dir=In)
lflux -= flux(states[j])
end
else
links = [Index(1; tags="Link,l=$n") for n in 1:N]
end

M[1] = cuITensor(T, sites[1], links[1])
M[1] += cuITensor(states[1] * state(links[1], 1))
for n in 2:(N - 1)
M[n] = cuITensor(T, dag(links[n - 1]), sites[n], links[n])
M[n] += cuITensor(state(dag(links[n - 1]), 1) * states[n] * state(links[n], 1))
end
M[N] = cuITensor(T, dag(links[N - 1]), sites[N])
M[N] += cuITensor(state(dag(links[N - 1]), 1) * states[N])

return M
end

function cuMPS(
::Type{T}, sites::Vector{<:Index}, state::Union{String,Integer}
) where {T<:Number}
return cuMPS(T, sites, fill(state, length(sites)))
end

function cuMPS(::Type{T}, sites::Vector{<:Index}, states::Function) where {T<:Number}
states_vec = [states(n) for n in 1:length(sites)]
return cuMPS(T, sites, states_vec)
end

cuMPS(sites::Vector{<:Index}, states) = cuMPS(Float64, sites, states)
cuMPO(M::MPO) = cu(M)
cuMPO(args...; kwargs...) = cu(MPO(args...; kwargs...))
randomCuMPO(args...; kwargs...) = cu(randomMPO(args...; kwargs...))
3 changes: 1 addition & 2 deletions ITensorGPU/src/tensor/cudense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,7 @@ function Dense{T,S}(x::T, size::Integer) where {T,S<:CuArray{<:T}}
fill!(arr, x)
return Dense{T,S}(arr)
end
cpu(x::CuDense{T}) where {T<:Number} = Dense(collect(x.data))
cpu(x::CuDenseTensor{T}) where {T<:Number} = Tensor(inds(x), cpu(store(x)))

function Base.complex(::Type{Dense{ElT,VT}}) where {ElT,VT<:CuArray}
return Dense{complex(ElT),CuVector{complex(ElT)}}
end
Expand Down
Loading

0 comments on commit d7d6ec8

Please sign in to comment.