From 49b67189bf610222893e5b317e604ed018fa3c24 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Thu, 22 Aug 2019 18:42:05 +0200
Subject: [PATCH 01/13] initial version o curand for randmat

---
 Project.toml                             |  9 +++++----
 REQUIRE                                  |  1 +
 curandommatrices/src/CuRandomMatrices.jl | 14 ++++++++++++++
 curandommatrices/src/ginibre.jl          | 10 ++++++++++
 curandommatrices/src/wigner.jl           |  4 ++++
 src/QuantumInformation.jl                |  9 ++++++++-
 6 files changed, 42 insertions(+), 5 deletions(-)
 create mode 100644 curandommatrices/src/CuRandomMatrices.jl
 create mode 100644 curandommatrices/src/ginibre.jl
 create mode 100644 curandommatrices/src/wigner.jl

diff --git a/Project.toml b/Project.toml
index e570e07..6c7a7a0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,20 +3,21 @@ uuid = "3c0b384b-479c-5684-b2ef-9d7a46dd931e"
 authors = ["Dariusz Kurzyk <dkurzyk@iitis.pl>", "Łukasz Pawela <lpawela@iitis.pl>", "Piotr Gawron <gawron@iitis.pl>", "Marcin Przewięźlikowski <m.przewie@gmail.com>"]
 version = "0.4.6"
 
-[compat]
-julia = "0.7, 1.0, 1.1, 1.2"
-
 [deps]
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
+[compat]
+julia = "0.7, 1.0, 1.1, 1.2"
+
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-        
+
 [targets]
 test = ["Test"]
diff --git a/REQUIRE b/REQUIRE
index a04bc71..c6bfad3 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -4,3 +4,4 @@ DocStringExtensions
 TensorOperations 
 Convex 
 SCS 
+Pkg
\ No newline at end of file
diff --git a/curandommatrices/src/CuRandomMatrices.jl b/curandommatrices/src/CuRandomMatrices.jl
new file mode 100644
index 0000000..dec6e6b
--- /dev/null
+++ b/curandommatrices/src/CuRandomMatrices.jl
@@ -0,0 +1,14 @@
+module CuRandomMatrices
+export curand
+
+using LinearAlgebra
+using CuArrays
+
+include("../../randommatrices/src/RandomMatrices.jl")
+using ..RandomMatrices
+
+include("ginibre.jl")
+# include("circular.jl")
+# include("wigner.jl")
+# include("wishart.jl")
+end
\ No newline at end of file
diff --git a/curandommatrices/src/ginibre.jl b/curandommatrices/src/ginibre.jl
new file mode 100644
index 0000000..052f217
--- /dev/null
+++ b/curandommatrices/src/ginibre.jl
@@ -0,0 +1,10 @@
+curand(g::GinibreEnsemble{1}) = CuArrays.randn(g.m, g.n)
+curand(g::GinibreEnsemble{2}) = CuArrays.randn( g.m, g.n)+1im*CuArrays.randn(g.m, g.n)
+
+function curand(g::GinibreEnsemble{4})
+    q0=CuArrays.randn(g.m, g.n)
+    q1=CuArrays.randn(g.m, g.n)
+    q2=CuArrays.randn(g.m, g.n)
+    q3=CuArrays.randn(g.m, g.n)
+    [q0+1im*q1 q2+1im*q3; -q2+1im*q3 q0-1im*q1]
+end
diff --git a/curandommatrices/src/wigner.jl b/curandommatrices/src/wigner.jl
new file mode 100644
index 0000000..515cd87
--- /dev/null
+++ b/curandommatrices/src/wigner.jl
@@ -0,0 +1,4 @@
+function curand(w::WignerEnsemble{β}) where β
+    z = curand(rng, GinibreEnsemble{β}(w.d))
+    (z + z') / 2sqrt(2β * w.d)
+end
\ No newline at end of file
diff --git a/src/QuantumInformation.jl b/src/QuantumInformation.jl
index 79e6c60..2b580df 100644
--- a/src/QuantumInformation.jl
+++ b/src/QuantumInformation.jl
@@ -8,7 +8,7 @@ using TensorOperations
 using Convex, SCS
 using Random: AbstractRNG, GLOBAL_RNG
 
-import Base: convert, size, length, kron, *, rand, show
+import Base: convert, size, length, kron, *, show, rand
 
 const ⊗ = kron
 
@@ -18,6 +18,13 @@ include("../randommatrices/src/RandomMatrices.jl")
 using .RandomMatrices
 eval(Expr(:export, names(RandomMatrices)...))
 
+using Pkg
+if "CuArrays" ∈ keys(Pkg.installed())
+    include("../curandommatrices/src/CuRandomMatrices.jl")
+    @eval using ..CuRandomMatrices
+    @eval export curand
+end
+
 include("base.jl")
 include("randomqobjects.jl")
 include("gates.jl")

From e52d5aae837345016d77f7e7a8cf3e58b6a2a0c9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Thu, 22 Aug 2019 18:52:27 +0200
Subject: [PATCH 02/13] fix wigner ensemble

---
 randommatrices/src/wigner.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/randommatrices/src/wigner.jl b/randommatrices/src/wigner.jl
index 0c7186c..5261619 100644
--- a/randommatrices/src/wigner.jl
+++ b/randommatrices/src/wigner.jl
@@ -2,16 +2,18 @@ export WignerEnsemble
 
 struct WignerEnsemble{β} <: QIContinuousMatrixDistribution
     d::Int
+    g::GinibreEnsemble{β}
 
     function WignerEnsemble{β}(d::Int) where β
         β == 4 && mod(d, 2) == 1 ? throw(ArgumentError("Dim must even")) : ()
-        new(d)
+        g = GinibreEnsemble{β}(d)
+        new(d, g)
     end
 end
 
 WignerEnsemble(d::Int) = WignerEnsemble{2}(d)
 
 function rand(rng::AbstractRNG, w::WignerEnsemble{β}) where β
-    z = rand(rng, GinibreEnsemble{β}(w.d))
+    z = rand(rng, w.g)
     (z + z') / 2sqrt(2β * w.d)
 end
\ No newline at end of file

From 08bb8fc76f7a8c2ffe8d3980082faacd09530f85 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Fri, 23 Aug 2019 02:10:25 +0200
Subject: [PATCH 03/13] almost working COE

---
 curandommatrices/src/CuRandomMatrices.jl |   8 +-
 curandommatrices/src/circular.jl         | 103 +++++++++++++++++++++++
 curandommatrices/src/wigner.jl           |   2 +-
 curandommatrices/src/wishart.jl          |   4 +
 4 files changed, 112 insertions(+), 5 deletions(-)
 create mode 100644 curandommatrices/src/circular.jl
 create mode 100644 curandommatrices/src/wishart.jl

diff --git a/curandommatrices/src/CuRandomMatrices.jl b/curandommatrices/src/CuRandomMatrices.jl
index dec6e6b..35dc464 100644
--- a/curandommatrices/src/CuRandomMatrices.jl
+++ b/curandommatrices/src/CuRandomMatrices.jl
@@ -2,13 +2,13 @@ module CuRandomMatrices
 export curand
 
 using LinearAlgebra
-using CuArrays
+using CuArrays, CUDAnative
 
 include("../../randommatrices/src/RandomMatrices.jl")
 using ..RandomMatrices
 
 include("ginibre.jl")
-# include("circular.jl")
-# include("wigner.jl")
-# include("wishart.jl")
+include("circular.jl")
+include("wigner.jl")
+include("wishart.jl")
 end
\ No newline at end of file
diff --git a/curandommatrices/src/circular.jl b/curandommatrices/src/circular.jl
new file mode 100644
index 0000000..9589663
--- /dev/null
+++ b/curandommatrices/src/circular.jl
@@ -0,0 +1,103 @@
+# FIXME: this can be accelerated by writing a custom kernel
+function angle!(a)
+    # @cuprintf("thread %ld", Int64(1))
+    ℜ = real.(a)
+    # @cuprintf("thread %ld", Int64(1))
+    ℑ = imag.(a)
+    # @cuprintf("thread %ld", Int64(1))
+    # @cuprintf("thread %ld", Int64(1))
+    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    # @cuprintf("thread %ld", Int64(1))
+    stride = blockDim().x * gridDim().x
+    # @cuprintf("thread %ld", Int64(1))
+    for i=index:stride:length(ℜ)
+        a[i] = (ℑ[i], ℜ[i])
+    end
+    return nothing
+end
+
+function _qr_fix!(z::CuMatrix)
+    q, r = CuArrays.qr!(z)
+    # @cuprintf("thread %ld", Int64(1))
+    ph = diag(r)
+    # @cuprintf("thread %ld", Int64(1))
+    @cuda threads=256 blocks=16 angle!(ph)
+    q = CuMatrix(q)
+    idim = size(r, 1)
+    for i=1:idim
+        q[:, i] .*= ph[i]
+    end
+    q[:, 1:idim]
+end
+
+function _qr_fix(z::CuMatrix)
+    a = copy(z)
+    _qr_fix!(a)
+end
+
+function curand(c::COE)
+    z = curand(c.g)
+    u = _qr_fix!(z)
+    transpose(u)*u
+end
+
+# function rand(rng::AbstractRNG, c::CUE)
+#     z = rand(rng, c.g)
+#     u = _qr_fix!(z)
+#     u
+# end
+
+# function rand(rng::AbstractRNG, c::CSE)
+#     z = rand(rng, c.g)
+#     u = _qr_fix!(z)
+#     ur = cat([[0 -1; 1 0] for _=1:c.d÷2]..., dims=[1,2])
+#     ur*u*ur'*transpose(u)
+# end
+
+# struct CircularRealEnsemble <: QIContinuousMatrixDistribution
+#     d::Int
+#     g::GinibreEnsemble{1}
+
+#     function CircularRealEnsemble(d::Int)
+#         g = GinibreEnsemble{1}(d)
+#         new(d, g)
+#     end
+# end
+
+# function rand(rng::AbstractRNG, c::CircularRealEnsemble)
+#     z = rand(rng, c.g)
+#     _qr_fix!(z)
+# end
+
+# struct CircularQuaternionEnsemble <: QIContinuousMatrixDistribution
+#     d::Int
+#     g::GinibreEnsemble{4}
+
+#     function CircularQuaternionEnsemble(d::Int)
+#         g = GinibreEnsemble{4}(d)
+#         new(d, g)
+#     end
+# end
+
+# function rand(rng::AbstractRNG, c::CircularQuaternionEnsemble)
+#     z = rand(rng, c.g)
+#     _qr_fix!(z)
+# end
+
+
+# struct HaarIsometry <: QIContinuousMatrixDistribution
+#     idim::Int
+#     odim::Int
+#     g::GinibreEnsemble{2}
+
+#     function HaarIsometry(idim::Int, odim::Int)
+#         idim <= odim || throw(ArgumentError("idim can't be greater than odim"))
+#         g = GinibreEnsemble{2}(odim, idim)
+#         new(idim, odim, g)
+#     end
+# end
+
+# function rand(rng::AbstractRNG, c::HaarIsometry)
+#     z = rand(rng, c.g)
+#     _qr_fix!(z)
+# end
\ No newline at end of file
diff --git a/curandommatrices/src/wigner.jl b/curandommatrices/src/wigner.jl
index 515cd87..b9a0aa2 100644
--- a/curandommatrices/src/wigner.jl
+++ b/curandommatrices/src/wigner.jl
@@ -1,4 +1,4 @@
 function curand(w::WignerEnsemble{β}) where β
-    z = curand(rng, GinibreEnsemble{β}(w.d))
+    z = curand(w.g)
     (z + z') / 2sqrt(2β * w.d)
 end
\ No newline at end of file
diff --git a/curandommatrices/src/wishart.jl b/curandommatrices/src/wishart.jl
new file mode 100644
index 0000000..202656d
--- /dev/null
+++ b/curandommatrices/src/wishart.jl
@@ -0,0 +1,4 @@
+function curand(w::WishartEnsemble{β, K}) where {β, K}
+    z = curand(w.g)/sqrt(2β * w.d)
+    z*z'
+end

From 6451b7f5760ae8d4d4bf9f06ee58aa422904c83a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lpawela@iitis.pl>
Date: Fri, 23 Aug 2019 17:45:05 +0200
Subject: [PATCH 04/13] add Requires to require

---
 Project.toml | 1 +
 REQUIRE      | 3 ++-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 6c7a7a0..e67a265 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
diff --git a/REQUIRE b/REQUIRE
index c6bfad3..021be06 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -4,4 +4,5 @@ DocStringExtensions
 TensorOperations 
 Convex 
 SCS 
-Pkg
\ No newline at end of file
+Pkg
+Requires
\ No newline at end of file

From 2e08c77238008d02a7d6d8d5de842661ec362a6f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lpawela@iitis.pl>
Date: Fri, 23 Aug 2019 17:45:29 +0200
Subject: [PATCH 05/13] Fix importing

---
 curandommatrices/src/CuRandomMatrices.jl |  3 ++-
 src/QuantumInformation.jl                | 10 ++++------
 2 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/curandommatrices/src/CuRandomMatrices.jl b/curandommatrices/src/CuRandomMatrices.jl
index 35dc464..c982c42 100644
--- a/curandommatrices/src/CuRandomMatrices.jl
+++ b/curandommatrices/src/CuRandomMatrices.jl
@@ -2,7 +2,8 @@ module CuRandomMatrices
 export curand
 
 using LinearAlgebra
-using CuArrays, CUDAnative
+using ..CuArrays
+using CUDAnative
 
 include("../../randommatrices/src/RandomMatrices.jl")
 using ..RandomMatrices
diff --git a/src/QuantumInformation.jl b/src/QuantumInformation.jl
index 2b580df..e5029e6 100644
--- a/src/QuantumInformation.jl
+++ b/src/QuantumInformation.jl
@@ -18,12 +18,10 @@ include("../randommatrices/src/RandomMatrices.jl")
 using .RandomMatrices
 eval(Expr(:export, names(RandomMatrices)...))
 
-using Pkg
-if "CuArrays" ∈ keys(Pkg.installed())
-    include("../curandommatrices/src/CuRandomMatrices.jl")
-    @eval using ..CuRandomMatrices
-    @eval export curand
-end
+using Requires
+@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" include("../curandommatrices/src/CuRandomMatrices.jl")
+@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" using ..CuRandomMatrices
+@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" export curand
 
 include("base.jl")
 include("randomqobjects.jl")

From 669b24f7612cec3999e19c9bc6b1954ea40cec2e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lpawela@iitis.pl>
Date: Fri, 23 Aug 2019 17:45:46 +0200
Subject: [PATCH 06/13] more erros in circular :(

---
 curandommatrices/src/circular.jl | 9 ++++++---
 curandommatrices/src/ginibre.jl  | 2 +-
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/curandommatrices/src/circular.jl b/curandommatrices/src/circular.jl
index 9589663..af34bb3 100644
--- a/curandommatrices/src/circular.jl
+++ b/curandommatrices/src/circular.jl
@@ -1,5 +1,6 @@
 # FIXME: this can be accelerated by writing a custom kernel
 function angle!(a)
+    a = CuArray(a)
     # @cuprintf("thread %ld", Int64(1))
     ℜ = real.(a)
     # @cuprintf("thread %ld", Int64(1))
@@ -10,9 +11,11 @@ function angle!(a)
     # @cuprintf("thread %ld", Int64(1))
     stride = blockDim().x * gridDim().x
     # @cuprintf("thread %ld", Int64(1))
-    for i=index:stride:length(ℜ)
-        a[i] = (ℑ[i], ℜ[i])
-    end
+    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
+    a[i] = ℜ[i] + ℑ[i]
+    # for i=index:stride:length(ℜ)
+    #     a[i] = +(ℑ[i], ℜ[i])
+    # end
     return nothing
 end
 
diff --git a/curandommatrices/src/ginibre.jl b/curandommatrices/src/ginibre.jl
index 052f217..3d8e138 100644
--- a/curandommatrices/src/ginibre.jl
+++ b/curandommatrices/src/ginibre.jl
@@ -1,5 +1,5 @@
 curand(g::GinibreEnsemble{1}) = CuArrays.randn(g.m, g.n)
-curand(g::GinibreEnsemble{2}) = CuArrays.randn( g.m, g.n)+1im*CuArrays.randn(g.m, g.n)
+curand(g::GinibreEnsemble{2}) = CuArrays.randn(g.m, g.n)+1im*CuArrays.randn(g.m, g.n)
 
 function curand(g::GinibreEnsemble{4})
     q0=CuArrays.randn(g.m, g.n)

From 2fd13e5bdaae78266eb6c9dd4b466d100d1d270e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Sat, 24 Aug 2019 02:05:20 +0200
Subject: [PATCH 07/13] working COE CUE CSE

---
 curandommatrices/src/circular.jl | 48 ++++++++++++--------------------
 1 file changed, 18 insertions(+), 30 deletions(-)

diff --git a/curandommatrices/src/circular.jl b/curandommatrices/src/circular.jl
index af34bb3..bad60f8 100644
--- a/curandommatrices/src/circular.jl
+++ b/curandommatrices/src/circular.jl
@@ -1,30 +1,18 @@
-# FIXME: this can be accelerated by writing a custom kernel
-function angle!(a)
-    a = CuArray(a)
-    # @cuprintf("thread %ld", Int64(1))
-    ℜ = real.(a)
-    # @cuprintf("thread %ld", Int64(1))
-    ℑ = imag.(a)
-    # @cuprintf("thread %ld", Int64(1))
-    # @cuprintf("thread %ld", Int64(1))
+# FIXME: this can be accelerated
+function cplx_phase!(a)
     index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
-    # @cuprintf("thread %ld", Int64(1))
-    stride = blockDim().x * gridDim().x
-    # @cuprintf("thread %ld", Int64(1))
     i = (blockIdx().x-1) * blockDim().x + threadIdx().x
-    a[i] = ℜ[i] + ℑ[i]
-    # for i=index:stride:length(ℜ)
-    #     a[i] = +(ℑ[i], ℜ[i])
-    # end
+    @inbounds a[i] = a[i] / sqrt(real(a[i])^2 + imag(a[i])^2)
     return nothing
 end
 
+# FIXME: use blocks to support larger than 1024
+
 function _qr_fix!(z::CuMatrix)
     q, r = CuArrays.qr!(z)
-    # @cuprintf("thread %ld", Int64(1))
     ph = diag(r)
-    # @cuprintf("thread %ld", Int64(1))
-    @cuda threads=256 blocks=16 angle!(ph)
+    len = length(ph)
+    @cuda threads=len cplx_phase!(ph)
     q = CuMatrix(q)
     idim = size(r, 1)
     for i=1:idim
@@ -44,18 +32,18 @@ function curand(c::COE)
     transpose(u)*u
 end
 
-# function rand(rng::AbstractRNG, c::CUE)
-#     z = rand(rng, c.g)
-#     u = _qr_fix!(z)
-#     u
-# end
+function curand(c::CUE)
+    z = curand(c.g)
+    u = _qr_fix!(z)
+    u
+end
 
-# function rand(rng::AbstractRNG, c::CSE)
-#     z = rand(rng, c.g)
-#     u = _qr_fix!(z)
-#     ur = cat([[0 -1; 1 0] for _=1:c.d÷2]..., dims=[1,2])
-#     ur*u*ur'*transpose(u)
-# end
+function curand(c::CSE)
+    z = curand(c.g)
+    u = _qr_fix!(z)
+    ur = cat([CuMatrix{Float32}([0 -1; 1 0]) for _=1:c.d÷2]..., dims=[1,2])
+    ur*u*ur'*transpose(u)
+end
 
 # struct CircularRealEnsemble <: QIContinuousMatrixDistribution
 #     d::Int

From f4f216de0dddca78b16436ed81a79e2d0037ffe4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Sat, 24 Aug 2019 22:27:08 +0200
Subject: [PATCH 08/13] finish circular ensemble

---
 curandommatrices/src/circular.jl | 60 +++++++-------------------------
 1 file changed, 13 insertions(+), 47 deletions(-)

diff --git a/curandommatrices/src/circular.jl b/curandommatrices/src/circular.jl
index bad60f8..5666c3e 100644
--- a/curandommatrices/src/circular.jl
+++ b/curandommatrices/src/circular.jl
@@ -34,8 +34,7 @@ end
 
 function curand(c::CUE)
     z = curand(c.g)
-    u = _qr_fix!(z)
-    u
+    _qr_fix!(z)
 end
 
 function curand(c::CSE)
@@ -45,50 +44,17 @@ function curand(c::CSE)
     ur*u*ur'*transpose(u)
 end
 
-# struct CircularRealEnsemble <: QIContinuousMatrixDistribution
-#     d::Int
-#     g::GinibreEnsemble{1}
-
-#     function CircularRealEnsemble(d::Int)
-#         g = GinibreEnsemble{1}(d)
-#         new(d, g)
-#     end
-# end
-
-# function rand(rng::AbstractRNG, c::CircularRealEnsemble)
-#     z = rand(rng, c.g)
-#     _qr_fix!(z)
-# end
-
-# struct CircularQuaternionEnsemble <: QIContinuousMatrixDistribution
-#     d::Int
-#     g::GinibreEnsemble{4}
-
-#     function CircularQuaternionEnsemble(d::Int)
-#         g = GinibreEnsemble{4}(d)
-#         new(d, g)
-#     end
-# end
-
-# function rand(rng::AbstractRNG, c::CircularQuaternionEnsemble)
-#     z = rand(rng, c.g)
-#     _qr_fix!(z)
-# end
-
-
-# struct HaarIsometry <: QIContinuousMatrixDistribution
-#     idim::Int
-#     odim::Int
-#     g::GinibreEnsemble{2}
+function curand(c::CircularRealEnsemble)
+    z = curand(c.g)
+    _qr_fix!(z)
+end
 
-#     function HaarIsometry(idim::Int, odim::Int)
-#         idim <= odim || throw(ArgumentError("idim can't be greater than odim"))
-#         g = GinibreEnsemble{2}(odim, idim)
-#         new(idim, odim, g)
-#     end
-# end
+function curand(c::CircularQuaternionEnsemble)
+    z = curand(c.g)
+    _qr_fix!(z)
+end
 
-# function rand(rng::AbstractRNG, c::HaarIsometry)
-#     z = rand(rng, c.g)
-#     _qr_fix!(z)
-# end
\ No newline at end of file
+function curand(c::HaarIsometry)
+    z = curand(c.g)
+    _qr_fix!(z)
+end
\ No newline at end of file

From b572cdb343147d985f9cad5b3de013101cb473d2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Sat, 24 Aug 2019 22:27:40 +0200
Subject: [PATCH 09/13] ?

---
 curandommatrices/test/circular.jl | 74 +++++++++++++++++++++++++++++++
 1 file changed, 74 insertions(+)
 create mode 100644 curandommatrices/test/circular.jl

diff --git a/curandommatrices/test/circular.jl b/curandommatrices/test/circular.jl
new file mode 100644
index 0000000..a7486da
--- /dev/null
+++ b/curandommatrices/test/circular.jl
@@ -0,0 +1,74 @@
+Random.seed!(42)
+
+@testset "CircularEnsemble" begin
+    @testset "CUE" begin
+        n=10
+        c = CUE(n)
+        u = CuRandomMatrices.curand(c)
+        @test norm(u*u' - I) ≈ 0 atol=1e-5
+
+        n = 100
+        c = CUE(n)
+        steps = 100
+        r = zeros(steps, n)
+
+        for i=1:steps
+            u = collect(CuRandomMatrices.curand(c))
+            r[i, :] = angle.(eigvals(u))
+        end
+        r = vec(r)
+        h = normalize(fit(Histogram, r, weights(ones(size(r))), -π:0.1π:π, closed=:left))
+        @test all(isapprox.(h.weights, 1/2π, atol=0.01))
+    end
+
+    @testset "COE" begin
+        n = 10
+        c = COE(n)
+        o = CuRandomMatrices.curand(c)
+        @test norm(o*o' - I) ≈ 0 atol=1e-5
+
+        n = 100
+        c = COE(n)
+        steps = 100
+        r = zeros(steps, n)
+        for i=1:steps
+            o = collect(CuRandomMatrices.curand(c))
+            r[i, :] = angle.(eigvals(o))
+        end
+        r = vec(r)
+        h = normalize(fit(Histogram, r, weights(ones(size(r))), -π:0.1π:π, closed=:left))
+        @test all(isapprox.(h.weights, 1/2π, atol=0.1))
+    end
+
+    @testset "CircularRealEnsemble" begin
+        c = CircularRealEnsemble(10)
+        o = CuRandomMatrices.curand(c)
+        @test size(o) == (10, 10)
+        @test eltype(o) <: Real
+    end
+end
+
+    @testset "HaarIsometry" begin
+        idim = 2
+        odim = 3
+        c = HaarIsometry(idim, odim)
+        u = CuRandomMatrices.curand(c)
+        @test size(u) == (odim, idim)
+        @test isapprox(norm(u'*u - I), 0, atol=1e-6)
+        @test_throws ArgumentError HaarIsometry(odim, idim)
+
+    @testset "CSE" begin
+        n = 10
+        c = CSE(n)
+        o = CuRandomMatrices.curand(c)
+        @test norm(o*o' - I) ≈ 0 atol=1e-5
+        @test size(o) == (n, n)
+    end
+
+    @testset "Circular quaternion ensemble" begin
+        c = CircularQuaternionEnsemble(10)
+        u = CuRandomMatrices.curand(c)
+        @test size(u) == (20, 20)
+        @test isapprox(norm(u'*u - I), 0, atol=1e-5)
+    end
+end
\ No newline at end of file

From 59098669f95ad068ceae0a43ee4d2ea245717c61 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Sat, 24 Aug 2019 22:28:03 +0200
Subject: [PATCH 10/13] add tests

---
 curandommatrices/test/ginibre.jl  | 23 +++++++++++++++++++++++
 curandommatrices/test/runtests.jl | 17 +++++++++++++++++
 curandommatrices/test/wigner.jl   | 11 +++++++++++
 curandommatrices/test/wishart.jl  | 15 +++++++++++++++
 4 files changed, 66 insertions(+)
 create mode 100644 curandommatrices/test/ginibre.jl
 create mode 100644 curandommatrices/test/runtests.jl
 create mode 100644 curandommatrices/test/wigner.jl
 create mode 100644 curandommatrices/test/wishart.jl

diff --git a/curandommatrices/test/ginibre.jl b/curandommatrices/test/ginibre.jl
new file mode 100644
index 0000000..2a43caf
--- /dev/null
+++ b/curandommatrices/test/ginibre.jl
@@ -0,0 +1,23 @@
+Random.seed!(42)
+
+@testset "GinibreEnsemble" begin
+    g = GinibreEnsemble{1}(10, 20)
+    z = rand(g)
+    @test eltype(z) <: Real
+    @test size(z) == (10, 20)
+    @test GinibreEnsemble(10, 20) == GinibreEnsemble{2}(10, 20)
+    @test GinibreEnsemble(10) == GinibreEnsemble{2}(10, 10)
+
+    @test_throws ArgumentError GinibreEnsemble{4}(11, 21)
+    g = GinibreEnsemble{4}(10, 20)
+    z = rand(g)
+    @test size(z) == (20, 40)
+    @test eltype(z) == ComplexF64
+
+@testset "_qr_fix" begin
+    a = rand(2, 2)
+    u1 = RandomMatrices._qr_fix(a)
+    u2 = RandomMatrices._qr_fix!(a)
+    @test norm(u1 - u2) ≈ 0
+end
+end
diff --git a/curandommatrices/test/runtests.jl b/curandommatrices/test/runtests.jl
new file mode 100644
index 0000000..487964c
--- /dev/null
+++ b/curandommatrices/test/runtests.jl
@@ -0,0 +1,17 @@
+using Test
+using Random
+using StatsBase
+using LinearAlgebra
+using CuArrays
+
+include("../../randommatrices/src/RandomMatrices.jl")
+using ..RandomMatrices
+
+include("../src/CuRandomMatrices.jl")
+using .CuRandomMatrices
+
+my_tests = ["circular.jl", "ginibre.jl", "wigner.jl", "wishart.jl"]
+
+for my_test in my_tests
+    include(my_test)
+end
diff --git a/curandommatrices/test/wigner.jl b/curandommatrices/test/wigner.jl
new file mode 100644
index 0000000..33ce00f
--- /dev/null
+++ b/curandommatrices/test/wigner.jl
@@ -0,0 +1,11 @@
+Random.seed!(42)
+
+@testset "WignerEnsemble" begin
+    w = WignerEnsemble{1}(10)
+    z = rand(w)
+
+    @test size(z) == (10, 10)
+    @test eltype(z) <: Real
+    @test ishermitian(z)
+    @test WignerEnsemble(10) == WignerEnsemble{2}(10)
+end
diff --git a/curandommatrices/test/wishart.jl b/curandommatrices/test/wishart.jl
new file mode 100644
index 0000000..c49e2ae
--- /dev/null
+++ b/curandommatrices/test/wishart.jl
@@ -0,0 +1,15 @@
+Random.seed!(42)
+
+@testset "WishartEnsemble" begin
+    w = WishartEnsemble{1, 0.1}(10)
+    z = rand(w)
+    ev = eigvals(z)
+
+    @test size(z) == (10, 10)
+    @test eltype(z) <: Real
+    @test length(ev[ev.>1e-5]) == 1
+    @test all(ev[ev.>1e-5] .> 0)
+
+    @test WishartEnsemble{1}(10) == WishartEnsemble{1, 1}(10)
+    @test WishartEnsemble(10) == WishartEnsemble{2, 1}(10)
+end

From 675145d9d388fba73fda30a3173f662af9946b7c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Sat, 24 Aug 2019 22:48:33 +0200
Subject: [PATCH 11/13] Allow arbitrary matrix dim

---
 curandommatrices/src/circular.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/curandommatrices/src/circular.jl b/curandommatrices/src/circular.jl
index 5666c3e..f578948 100644
--- a/curandommatrices/src/circular.jl
+++ b/curandommatrices/src/circular.jl
@@ -1,8 +1,9 @@
 # FIXME: this can be accelerated
 function cplx_phase!(a)
-    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     i = (blockIdx().x-1) * blockDim().x + threadIdx().x
-    @inbounds a[i] = a[i] / sqrt(real(a[i])^2 + imag(a[i])^2)
+    if i <= length(a)
+        @inbounds a[i] = a[i] / sqrt(real(a[i])^2 + imag(a[i])^2)
+    end
     return nothing
 end
 
@@ -11,8 +12,8 @@ end
 function _qr_fix!(z::CuMatrix)
     q, r = CuArrays.qr!(z)
     ph = diag(r)
-    len = length(ph)
-    @cuda threads=len cplx_phase!(ph)
+    len = min(length(ph), 1024)
+    @cuda threads=len blocks=16 cplx_phase!(ph)
     q = CuMatrix(q)
     idim = size(r, 1)
     for i=1:idim

From a95c835eb18862da23b8f3d4bbb3638222e4da01 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Sun, 25 Aug 2019 13:55:31 +0200
Subject: [PATCH 12/13] minor stylistic fixex

---
 curandommatrices/test/circular.jl | 2 +-
 randommatrices/src/ginibre.jl     | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/curandommatrices/test/circular.jl b/curandommatrices/test/circular.jl
index a7486da..d7fbcdb 100644
--- a/curandommatrices/test/circular.jl
+++ b/curandommatrices/test/circular.jl
@@ -2,7 +2,7 @@ Random.seed!(42)
 
 @testset "CircularEnsemble" begin
     @testset "CUE" begin
-        n=10
+        n = 10
         c = CUE(n)
         u = CuRandomMatrices.curand(c)
         @test norm(u*u' - I) ≈ 0 atol=1e-5
diff --git a/randommatrices/src/ginibre.jl b/randommatrices/src/ginibre.jl
index b790806..bbaea63 100644
--- a/randommatrices/src/ginibre.jl
+++ b/randommatrices/src/ginibre.jl
@@ -19,9 +19,9 @@ rand(rng::AbstractRNG, g::GinibreEnsemble{2}) = randn(rng, g.m, g.n)+1im*randn(r
 
 function rand(rng::AbstractRNG, g::GinibreEnsemble{4})
     # TODO: fix dimensions of blocks
-    q0=randn(rng, g.m, g.n)
-    q1=randn(rng, g.m, g.n)
-    q2=randn(rng, g.m, g.n)
-    q3=randn(rng, g.m, g.n)
+    q0 = randn(rng, g.m, g.n)
+    q1 = randn(rng, g.m, g.n)
+    q2 = randn(rng, g.m, g.n)
+    q3 = randn(rng, g.m, g.n)
     [q0+1im*q1 q2+1im*q3; -q2+1im*q3 q0-1im*q1]
 end

From 235400fe37ea674a76cedb16aeb14ccbcd600410 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <lukasz.pawela@gmail.com>
Date: Sun, 25 Aug 2019 13:55:57 +0200
Subject: [PATCH 13/13] loop over types

---
 curandommatrices/src/circular.jl | 29 +++++++++--------------------
 1 file changed, 9 insertions(+), 20 deletions(-)

diff --git a/curandommatrices/src/circular.jl b/curandommatrices/src/circular.jl
index f578948..96da3fe 100644
--- a/curandommatrices/src/circular.jl
+++ b/curandommatrices/src/circular.jl
@@ -12,11 +12,11 @@ end
 function _qr_fix!(z::CuMatrix)
     q, r = CuArrays.qr!(z)
     ph = diag(r)
-    len = min(length(ph), 1024)
+    len = min(length(ph), 1024) #hack, warpsize() segfaults
     @cuda threads=len blocks=16 cplx_phase!(ph)
     q = CuMatrix(q)
     idim = size(r, 1)
-    for i=1:idim
+    for i = 1:idim
         q[:, i] .*= ph[i]
     end
     q[:, 1:idim]
@@ -33,11 +33,6 @@ function curand(c::COE)
     transpose(u)*u
 end
 
-function curand(c::CUE)
-    z = curand(c.g)
-    _qr_fix!(z)
-end
-
 function curand(c::CSE)
     z = curand(c.g)
     u = _qr_fix!(z)
@@ -45,17 +40,11 @@ function curand(c::CSE)
     ur*u*ur'*transpose(u)
 end
 
-function curand(c::CircularRealEnsemble)
-    z = curand(c.g)
-    _qr_fix!(z)
-end
-
-function curand(c::CircularQuaternionEnsemble)
-    z = curand(c.g)
-    _qr_fix!(z)
-end
-
-function curand(c::HaarIsometry)
-    z = curand(c.g)
-    _qr_fix!(z)
+for T in (CUE, CircularRealEnsemble, CircularQuaternionEnsemble, HaarIsometry)
+    @eval begin
+        function curand(c::$T)
+            z = curand(c.g)
+            _qr_fix!(z)
+        end
+    end
 end
\ No newline at end of file