From 9406a8fadedbb6f644cb85762cdce3cae55afe56 Mon Sep 17 00:00:00 2001
From: Will Tebbutt <will.tebbutt@invenialabs.co.uk>
Date: Thu, 12 Jan 2017 16:47:10 -0600
Subject: [PATCH 01/10] Basic DFT functionality added. Coverage of the expected
 operations is, however, somewhat patchy.

---
 src/SpecialMatrices.jl |  3 +-
 src/dft.jl             | 98 ++++++++++++++++++++++++++++++++++++++++++
 test/dft.jl            |  7 +++
 test/runtests.jl       |  1 +
 4 files changed, 108 insertions(+), 1 deletion(-)
 create mode 100644 src/dft.jl
 create mode 100644 test/dft.jl

diff --git a/src/SpecialMatrices.jl b/src/SpecialMatrices.jl
index 58ebfa1..642743e 100644
--- a/src/SpecialMatrices.jl
+++ b/src/SpecialMatrices.jl
@@ -1,9 +1,10 @@
 module SpecialMatrices
 
-import Base: A_mul_B!, full, getindex, inv, isassigned, size, *
+import Base: A_mul_B!, full, getindex, inv, isassigned, size, *, length
 
 include("cauchy.jl") #Cauchy matrix
 include("companion.jl") #Companion matrix
+include("dft.jl") #Discrete Fourier Transform matrices.
 include("frobenius.jl") #Frobenius matrix
 include("hankel.jl") #Hankel matrix
 include("hilbert.jl") #Hilbert matrix
diff --git a/src/dft.jl b/src/dft.jl
new file mode 100644
index 0000000..a6ec74b
--- /dev/null
+++ b/src/dft.jl
@@ -0,0 +1,98 @@
+export DFT, IDFT
+
+# Implementation of the DFT matrix. Not really an AbstractArray as it doesn't
+# make sense to write to this matrix. Should probably resolve this.
+
+# Non-unitary DFT matrices.
+abstract AbstractDFT
+immutable DFT <: AbstractDFT
+    N::Int
+end
+immutable IDFT <: AbstractDFT
+    N::Int
+end
+
+# Unitary DFT matrices.
+abstract AbstractUnitaryDFT <: AbstractDFT
+immutable UnitaryDFT <: AbstractUnitaryDFT
+    N::Int
+    rootN::Float64
+    UnitaryDFT(N) = new(N, sqrt(N))
+end
+immutable UnitaryIDFT <: AbstractUnitaryDFT
+    N::Int
+    rootN::Float64
+    UnitaryIDFT(N) = new(N, sqrt(N))
+end
+
+size(M::AbstractDFT) = (M.N, M.N)
+length(M::AbstractDFT) = M.N^2
+
+getindex(M::DFT, m::Int, n::Int) =
+    exp(-2π * im * (m-1) * (n-1) / M.N)
+getindex(N::IDFT, m::Int, n::Int) =
+    exp(2π * im * (m-1) * (n-1) / M.N) / M.N
+getindex(M::UnitaryDFT, m::Int, n::Int) =
+    exp(-2π * im * (m-1) * (n-1) / M.N) / sqrt(M.N)
+getindex(M::UnitaryIDFT, m::Int, n::Int) =
+    exp(2π * im * (m-1) * (n-1) / M.N) / sqrt(M.N)
+getindex(M::AbstractDFT, m::Int) = getindex(M, rem(m, M.N), div(m, M.N))
+
+isassigned(M::AbstractDFT, i::Int) = 1 <= i <= M.N ? true : false
+
+# Helper functions to check whether stuff is conformal or not.
+function conformal(U::AbstractDFT, x::Vector)
+    if U.N != length(x)
+        throw(ArgumentError("U and x are not conformal."))
+    end
+end
+function conformal(U::AbstractDFT, X::Matrix)
+    if U.N != size(X, 1)
+        throw(ArgumentError("U and X are not conformal."))
+    end
+end
+function conformal(X::Matrix, U::AbstractDFT)
+    if size(X, 2) != U.N
+        throw(ArgumentError("X and U are not conformal."))
+    end
+end
+
+# Multiplication from the left by the DFT matrices.
+*(U::DFT, x::Vector) = (conformal(U, x); fft(x))
+*(U::IDFT, x::Vector) = (conformal(U, x); ifft(x))
+*(U::UnitaryDFT, x::Vector) = (conformal(U, x); fft(x) / U.rootN)
+*(U::UnitaryIDFT, x::Vector) = (conformal(U, x); bfft(x) / U.rootN)
+
+*(U::DFT, X::Matrix) = (conformal(U, X); fft(X, 1))
+*(U::IDFT, X::Matrix) = (conformal(U, X); ifft(X, 1))
+*(U::UnitaryDFT, X::Matrix) = (conformal(U, X); fft(X, 1) / U.rootN)
+*(U::UnitaryIDFT, X::Matrix) = (conformal(U, X);  bfft(X, 1) / U.rootN)
+
+# Multiplication from the right by the DFT matrices.
+*(X::Matrix, U::DFT) = (conformal(X, U); fft(X, 2))
+*(X::Matrix, U::IDFT) = (conformal(X, U); ifft(X, 2))
+*(X::Matrix, U::UnitaryDFT) = (conformal(X, U); fft(X, 2) / U.rootN)
+*(X::Matrix, U::UnitaryIDFT) = (conformal(X, U); bfft(X, 2) / U.rootN)
+
+# In-place functionality.
+A_mul_B!(Y::Matrix, U::DFT, X::Matrix) =
+    (conformal(U, X); copy!(Y, X); fft!(Y, 1))
+A_mul_B!(Y::Matrix, U::IDFT, X::Matrix) =
+    (conformal(U, X); copy!(Y, X); ifft!(Y, 1))
+A_mul_B!(Y::Matrix, U::UnitaryDFT, X::Matrix) =
+    (conformal(U, X); copy!(Y, X); fft!(Y, 1); y ./= U.rootN)
+A_mul_B!(Y::Matrix, U::UnitaryIDFT, X::Matrix) =
+    (conformal(U, X); copy!(Y, X); bfft!(Y, 1); y ./= U.rootN)
+
+# Conjugate transpose of DFT times X.
+Ac_mul_B(U::DFT, X::Matrix) = bfft(X, 1)
+Ac_mul_B(U::IDFT, X::Matrix) = fft(X, 1) / U.N
+Ac_mul_B(U::UnitaryDFT, X::Matrix) = bfft(X, 1) / U.rootN
+Ac_mul_B(U::UnitaryIDFT, X::Matrix) = fft(X, 1) / U.rootN
+
+# All DFT matrices are symmetric.
+At_mul_B(U::AbstractDFT, X::Matrix) = U * X
+
+
+
+
diff --git a/test/dft.jl b/test/dft.jl
new file mode 100644
index 0000000..8148450
--- /dev/null
+++ b/test/dft.jl
@@ -0,0 +1,7 @@
+function check_vector()
+    N = 5
+    x = randn(N)
+    U = DFT(N)
+    all((U * x) .== (U * reshape(x, N, 1)))
+end
+@test check_vector()
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 10e8316..f59515c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,7 @@
 using SpecialMatrices
 using Base.Test  
 
+include("dft.jl") # Discrete Fourier Transform matrix.
 include("companion.jl") #Companion matrix
 include("frobenius.jl") #Frobenius matrix
 include("strang.jl") #Strang matrix

From ff5320a4f0e7cb52d04b11bad6aa132f15c504af Mon Sep 17 00:00:00 2001
From: Will Tebbutt <will.tebbutt@invenialabs.co.uk>
Date: Thu, 12 Jan 2017 17:54:39 -0600
Subject: [PATCH 02/10] Re-implmented existing Circulant multiplication
 functionality to correctly use dispatch to determine what to cast back to.
 Previous implementation used dynamic type checking and missed some edge
 cases. Not re-implemented the in-place multiplication functionality in this
 commit, nor have the tests been updated. Also added typealias for Circulant
 to reduce verbosity.

---
 src/dft.jl      | 57 ++++++++++++++++++++++++++++++++++-----------
 src/toeplitz.jl | 62 ++++++++++++++++++++++++++++++++++++++-----------
 2 files changed, 93 insertions(+), 26 deletions(-)

diff --git a/src/dft.jl b/src/dft.jl
index a6ec74b..3a983cf 100644
--- a/src/dft.jl
+++ b/src/dft.jl
@@ -1,4 +1,4 @@
-export DFT, IDFT
+export DFT, IDFT, UnitaryDFT, UnitaryIDFT
 
 # Implementation of the DFT matrix. Not really an AbstractArray as it doesn't
 # make sense to write to this matrix. Should probably resolve this.
@@ -16,12 +16,12 @@ end
 abstract AbstractUnitaryDFT <: AbstractDFT
 immutable UnitaryDFT <: AbstractUnitaryDFT
     N::Int
-    rootN::Float64
+    sqrtN::Float64
     UnitaryDFT(N) = new(N, sqrt(N))
 end
 immutable UnitaryIDFT <: AbstractUnitaryDFT
     N::Int
-    rootN::Float64
+    sqrtN::Float64
     UnitaryIDFT(N) = new(N, sqrt(N))
 end
 
@@ -60,35 +60,53 @@ end
 # Multiplication from the left by the DFT matrices.
 *(U::DFT, x::Vector) = (conformal(U, x); fft(x))
 *(U::IDFT, x::Vector) = (conformal(U, x); ifft(x))
-*(U::UnitaryDFT, x::Vector) = (conformal(U, x); fft(x) / U.rootN)
-*(U::UnitaryIDFT, x::Vector) = (conformal(U, x); bfft(x) / U.rootN)
+*(U::UnitaryDFT, x::Vector) = (conformal(U, x); fft(x) / U.sqrtN)
+*(U::UnitaryIDFT, x::Vector) = (conformal(U, x); bfft(x) / U.sqrtN)
 
 *(U::DFT, X::Matrix) = (conformal(U, X); fft(X, 1))
 *(U::IDFT, X::Matrix) = (conformal(U, X); ifft(X, 1))
-*(U::UnitaryDFT, X::Matrix) = (conformal(U, X); fft(X, 1) / U.rootN)
-*(U::UnitaryIDFT, X::Matrix) = (conformal(U, X);  bfft(X, 1) / U.rootN)
+*(U::UnitaryDFT, X::Matrix) = (conformal(U, X); fft(X, 1) / U.sqrtN)
+*(U::UnitaryIDFT, X::Matrix) = (conformal(U, X);  bfft(X, 1) / U.sqrtN)
 
 # Multiplication from the right by the DFT matrices.
 *(X::Matrix, U::DFT) = (conformal(X, U); fft(X, 2))
 *(X::Matrix, U::IDFT) = (conformal(X, U); ifft(X, 2))
-*(X::Matrix, U::UnitaryDFT) = (conformal(X, U); fft(X, 2) / U.rootN)
-*(X::Matrix, U::UnitaryIDFT) = (conformal(X, U); bfft(X, 2) / U.rootN)
+*(X::Matrix, U::UnitaryDFT) = (conformal(X, U); fft(X, 2) / U.sqrtN)
+*(X::Matrix, U::UnitaryIDFT) = (conformal(X, U); bfft(X, 2) / U.sqrtN)
 
 # In-place functionality.
+A_mul_B!(y::Vector, U::DFT, x::Vector) =
+    (conformal(U, x); copy!(y, x); fft!(y))
+A_mul_B!(y::Vector, U::IDFT, x::Vector) =
+    (conformal(U, x); copy!(y, x); ifft!(y))
+A_mul_B!(y::Vector, U::UnitaryDFT, x::Vector) =
+    (conformal(U, x); copy!(y, x); fft!(y); y ./= U.sqrtN)
+A_mul_B!(y::Vector, U::UnitaryIDFT, x::Vector) =
+    (conformal(U, x); copy!(y, x); bfft!(y); y ./= U.sqrtN)
+
 A_mul_B!(Y::Matrix, U::DFT, X::Matrix) =
     (conformal(U, X); copy!(Y, X); fft!(Y, 1))
 A_mul_B!(Y::Matrix, U::IDFT, X::Matrix) =
     (conformal(U, X); copy!(Y, X); ifft!(Y, 1))
 A_mul_B!(Y::Matrix, U::UnitaryDFT, X::Matrix) =
-    (conformal(U, X); copy!(Y, X); fft!(Y, 1); y ./= U.rootN)
+    (conformal(U, X); copy!(Y, X); fft!(Y, 1); y ./= U.sqrtN)
 A_mul_B!(Y::Matrix, U::UnitaryIDFT, X::Matrix) =
-    (conformal(U, X); copy!(Y, X); bfft!(Y, 1); y ./= U.rootN)
+    (conformal(U, X); copy!(Y, X); bfft!(Y, 1); y ./= U.sqrtN)
+
+A_mul_B!(Y::Matrix, X::Matrix, U::DFT) =
+    (conformal(X, U); copy!(Y, X); fft!(Y, 2))
+A_mul_B!(Y::Matrix, X::Matrix, U::IDFT) =
+    (conformal(X, U); copy!(Y, X); ifft!(Y, 2))
+A_mul_B!(Y::Matrix, X::Matrix, U::UnitaryDFT) =
+    (conformal(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
+A_mul_B!(Y::Matrix, X::Matrix, U::UnitaryIDFT) =
+    (conformal(X, U); copy!(Y, X); bfft!(Y, 2); y ./= U.sqrtN)
 
 # Conjugate transpose of DFT times X.
 Ac_mul_B(U::DFT, X::Matrix) = bfft(X, 1)
 Ac_mul_B(U::IDFT, X::Matrix) = fft(X, 1) / U.N
-Ac_mul_B(U::UnitaryDFT, X::Matrix) = bfft(X, 1) / U.rootN
-Ac_mul_B(U::UnitaryIDFT, X::Matrix) = fft(X, 1) / U.rootN
+Ac_mul_B(U::UnitaryDFT, X::Matrix) = bfft(X, 1) / U.sqrtN
+Ac_mul_B(U::UnitaryIDFT, X::Matrix) = fft(X, 1) / U.sqrtN
 
 # All DFT matrices are symmetric.
 At_mul_B(U::AbstractDFT, X::Matrix) = U * X
@@ -96,3 +114,16 @@ At_mul_B(U::AbstractDFT, X::Matrix) = U * X
 
 
 
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/toeplitz.jl b/src/toeplitz.jl
index 711974b..e229150 100644
--- a/src/toeplitz.jl
+++ b/src/toeplitz.jl
@@ -41,25 +41,61 @@ end
 
 
 
-immutable Circulant{T} <: AbstractArray{T, 2}
-	c :: Vector{T}
+immutable Circulant{T} <: AbstractMatrix{T}
+	c::Vector{T}
 end
 
-getindex(C::Circulant, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
-isassigned(C::Circulant, i::Int, j::Int) = isassigned(C.c,mod(i-j,length(C.c))+1)
-size(C::Circulant, r::Int) = (r==1 || r==2) ? length(C.c) :
+typealias Circ Circulant
+
+getindex(C::Circ, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
+isassigned(C::Circ, i::Int, j::Int) = isassigned(C.c,mod(i-j,length(C.c))+1)
+size(C::Circ, r::Int) = (r==1 || r==2) ? length(C.c) :
     throw(ArgumentError("Invalid dimension $r"))
-size(C::Circulant) = size(C,1), size(C,2)
+size(C::Circ) = size(C,1), size(C,2)
 
-# Fast matrix x vector via fft()
-# see Golub, van Loan, Matrix Computations, John Hopkins, Baltimore, 1996, p. 202 
-function *{T}(C::Circulant{T},x::Vector{T})
-    xt=fft(x)
-    vt=fft(C.c)
-    yt=vt.*xt
-    typeof(x[1])==Int ? map(Int,round(real(ifft(yt)))): ( (T <: Real) ? map(T,real(ifft(yt))) : ifft(yt))
+# Generic code for multiplication by a Circulant matrix. Dispatch is used later
+# on to ensure that the types returned are the correct ones. The idea is that
+# if both inputs are Integer, then the output should also be integer.
+# Similarly, if the inputs are both real, then the output should also be real.
+function circ_dot(C::Circ, X::VecOrMat)
+    if size(X, 1) != size(C.c, 1)
+        throw(ArgumentError("C and X are not conformal."))
+    end
+    U = DFT(size(X, 1))
+    return U * Diagonal(U * C.c) * U'X
+end
+function circ_dot(X::VecOrMat, C::Circ)
+    if size(X, 2) != size(C.c, 1)
+        throw(ArgumentError("C and X are not conformal."))
+    end
+    U = DFT(size(X, 1))
+    return X * U * Diagonal(U * C.c) * U'
 end
 
+# In general, the return type of the multiplication operation should be a
+# complex float.
+*(C::Circ, X::VecOrMat) = circ_dot(C, X)
+*(X::VecOrMat, C::Circ) = circ_dot(X, C)
+
+# If both arguments contain reals, then the result should really contain reals.
+*{T<:Real, V<:Real}(C::Circ{T}, X::VecOrMat{V}) = real(circ_dot(C, X))
+*{T<:Real, V<:Real}(X::VecOrMat{V}, C::Circ{T}) = real(circ_dot(X, C))
+
+# If both arguments contain (real) integers, then the result should be both
+# Integer-valued and real. Promote to common type for the purposes of
+# converting back to integer later.
+*{T<:Integer, V<:Integer}(C::Circ{T}, X::VecOrMat{V}) = *(promote(C, X)...)
+*{T<:Integer, V<:Integer}(X::VecOrMat{V}, C::Circ{T}) = *(promote(C, X)...)
+*{T<:Integer}(C::Circ{T}, X::VecOrMat{T}) =
+    convert(T, round.(real.(circ_dot(C, X))))
+*{T<:Integer}(X::VecOrMat{T}, C::Circ{T}) =
+    convert(T, round.(real.(circ_dot(C, X))))
+
+# TODO: Create functionality to allow complex integers to be appropriately handled.
+
+
+
+
 function A_mul_B!{T}(y::StridedVector{T},C::Circulant{T},x::StridedVector{T})
     xt=fft(x)
     vt=fft(C.c)

From 20ed4b4c8f1725aa19851f3dc094fce36318d051 Mon Sep 17 00:00:00 2001
From: Will Tebbutt <will.tebbutt@invenialabs.co.uk>
Date: Sat, 14 Jan 2017 00:59:32 -0600
Subject: [PATCH 03/10] Completed + tested basic functionality for circulant
 matrices, including multiplication in terms of DFT matrices, (conjugate)
 transposition, addition, subtraction, eigenvalues + eigenvectors. dft.jl has
 been further extended to cover more functionality and a few bugs have been
 fixed.

---
 src/SpecialMatrices.jl |  30 ++++----
 src/circulant.jl       |  94 ++++++++++++++++++++++++
 src/dft.jl             | 157 ++++++++++++++---------------------------
 src/toeplitz.jl        |  83 +---------------------
 test/circulant.jl      |  81 +++++++++++++++++++++
 test/dft.jl            |  24 +++++--
 test/runtests.jl       |   7 +-
 7 files changed, 270 insertions(+), 206 deletions(-)
 create mode 100644 src/circulant.jl
 create mode 100644 test/circulant.jl

diff --git a/src/SpecialMatrices.jl b/src/SpecialMatrices.jl
index 642743e..8683de5 100644
--- a/src/SpecialMatrices.jl
+++ b/src/SpecialMatrices.jl
@@ -1,17 +1,23 @@
 module SpecialMatrices
 
-import Base: A_mul_B!, full, getindex, inv, isassigned, size, *, length
+import Base: A_mul_B!, full, getindex, inv, isassigned, size, *, length, +, -,
+    Ac_mul_B, A_mul_Bc, At_mul_B, A_mul_Bt, eig
 
-include("cauchy.jl") #Cauchy matrix
-include("companion.jl") #Companion matrix
-include("dft.jl") #Discrete Fourier Transform matrices.
-include("frobenius.jl") #Frobenius matrix
-include("hankel.jl") #Hankel matrix
-include("hilbert.jl") #Hilbert matrix
-include("kahan.jl") #Kahan matrix
-include("riemann.jl") #Riemann matrix
-include("strang.jl") #Strang matrix
-include("toeplitz.jl") #Toeplitz matrix, Circulant matrix
-include("vandermonde.jl") #Vandermonde matrix
+typealias SV StridedVector
+typealias SM StridedMatrix
+typealias SVM StridedVecOrMat
+
+include("cauchy.jl") # Cauchy matrix
+include("circulant.jl") # Circulant matrix.
+include("companion.jl") # Companion matrix
+include("dft.jl") # Discrete Fourier Transform matrices.
+include("frobenius.jl") # Frobenius matrix
+include("hankel.jl") # Hankel matrix
+include("hilbert.jl") # Hilbert matrix
+include("kahan.jl") # Kahan matrix
+include("riemann.jl") # Riemann matrix
+include("strang.jl") # Strang matrix
+include("toeplitz.jl") # Toeplitz matrix
+include("vandermonde.jl") # Vandermonde matrix
 
 end # module
diff --git a/src/circulant.jl b/src/circulant.jl
new file mode 100644
index 0000000..f3d0e53
--- /dev/null
+++ b/src/circulant.jl
@@ -0,0 +1,94 @@
+export Circulant
+
+immutable Circulant{T} <: AbstractMatrix{T}
+    c::Vector{T}
+end
+typealias Circ Circulant
+
+function full{T}(C::Circ{T})
+    n = size(C, 1)
+    M = Array(T, n, n)
+    for i=1:n
+        M[i:n,i] = C.c[1:n-i+1]
+        M[1:i-1,i] = C.c[n-i+2:n]
+    end
+    M
+end
+
+getindex(C::Circ, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
+isassigned(C::Circ, i::Int, j::Int) = isassigned(C.c,mod(i-j,length(C.c))+1)
+size(C::Circ, r::Int) = (r==1 || r==2) ? length(C.c) :
+    throw(ArgumentError("Invalid dimension $r"))
+size(C::Circ) = size(C,1), size(C,2)
+
+Base.copy(C::Circ) = Circ(C.c)
+Base.conj(C::Circ) = Circ(conj(C.c))
+Base.conj!(C::Circ) = (conj!(C.c); C)
+
+Base.transpose(C::Circ) = Circ(circshift(reverse(C.c), 1))
+Base.ctranspose(C::Circ) = conj!(transpose(C))
+
++(C::Circ, a::Number) = Circ(C.c + a)
++(a::Number, C::Circ) = Circ(a + C.c)
++(C::Circ, D::Circ) = Circ(C.c + D.c)
+-(C::Circ, a::Number) = Circ(C.c - a)
+-(a::Number, C::Circ) = Circ(a - C.c)
+-(C::Circ, D::Circ) = Circ(C.c - D.c)
+-(C::Circ) = Circ(-C.c)
+
+# # Return the eigen-factorisation of C. Constituents are efficient to compute
+# # and the eigenvectors are represented efficiently using a UnitaryDFT object.
+# function Base.eigfact{T}(C::Circ{T})
+#     n = length(C.c)
+#     Base.LinAlg.Eigen(DFT{T}(n) * C.c, UnitaryDFT{T}(n))
+# end
+
+eig{T}(C::Circ{T}) = (N = length(C.c); (DFT{T}(N) * C.c, UnitaryDFT{T}(N)))
+
+# Generic code for multiplication by a Circulant matrix. Dispatch is used later
+# on to ensure that the types returned are the correct ones. The idea is that
+# if both inputs are Integer, then the output should also be integer.
+# Similarly, if the inputs are both real, then the output should also be real.
+function circ_dot{T}(C::Circ{T}, X::StridedVecOrMat)
+    if size(X, 1) != size(C, 2)
+        throw(ArgumentError("C and X are not conformal."))
+    end
+    γ, U = eig(C)
+    return Ac_mul_B(U, (Diagonal(γ) * (U * X)))
+end
+function circ_dot(X::StridedVecOrMat, C::Circ)
+    if size(X, 2) != size(C, 1)
+        println(size(X, 2))
+        println(size(C, 1))
+        throw(ArgumentError("C and X are not conformal."))
+    end
+    γ, U = eig(C)
+    return (A_mul_Bc(X, U) * Diagonal(γ)) * U
+end
+
+# In general, the return type of the * operation should be a complex float.
+*(C::Circ, X::SVM) = circ_dot(C, X)
+*(X::SVM, C::Circ) = circ_dot(X, C)
+
+# If both arguments contain reals, then the result should really contain reals.
+*{T<:Real, V<:Real}(C::Circ{T}, X::SVM{V}) = real!.(circ_dot(C, X))
+*{T<:Real, V<:Real}(X::SVM{V}, C::Circ{T}) = real!.(circ_dot(X, C))
+
+# TODO: How attached are people to this? This should probably change. Also, this doesn't
+# really get you anything in terms of performance. There's not reduction in memory use as far as I can see
+# as the in-place fourier transform is not used.
+function A_mul_B!{T}(y::StridedVector{T},C::Circulant{T},x::StridedVector{T})
+    xt=fft(x)
+    vt=fft(C.c)
+    yt=ifft(vt.*xt)
+    if T<: Int
+        map!(round,y,yt) 
+    elseif T<: Real
+        map!(real,y,yt)
+    else
+        copy!(y,yt)
+    end
+    return y
+end
+
+# TODO: Implement all of the other functionality that you would expect to find with a matrix.
diff --git a/src/dft.jl b/src/dft.jl
index 3a983cf..945ecec 100644
--- a/src/dft.jl
+++ b/src/dft.jl
@@ -1,129 +1,80 @@
-export DFT, IDFT, UnitaryDFT, UnitaryIDFT
+export DFT, UnitaryDFT
 
-# Implementation of the DFT matrix. Not really an AbstractArray as it doesn't
-# make sense to write to this matrix. Should probably resolve this.
-
-# Non-unitary DFT matrices.
-abstract AbstractDFT
-immutable DFT <: AbstractDFT
-    N::Int
-end
-immutable IDFT <: AbstractDFT
+# Non-unitary DFT matrix.
+abstract AbstractDFT{T} <: AbstractMatrix{T}
+immutable DFT{T} <: AbstractDFT{T}
     N::Int
 end
 
-# Unitary DFT matrices.
-abstract AbstractUnitaryDFT <: AbstractDFT
-immutable UnitaryDFT <: AbstractUnitaryDFT
+# Unitary DFT Matrix.
+immutable UnitaryDFT{T} <: AbstractDFT{T}
     N::Int
     sqrtN::Float64
     UnitaryDFT(N) = new(N, sqrt(N))
 end
-immutable UnitaryIDFT <: AbstractUnitaryDFT
-    N::Int
-    sqrtN::Float64
-    UnitaryIDFT(N) = new(N, sqrt(N))
-end
 
-size(M::AbstractDFT) = (M.N, M.N)
+Base.show(io::IO, U::DFT) = print(io, "DFT matrix.")
+Base.show(io::IO, U::UnitaryDFT) = print(io, "Unitary DFT matrix.")
+
+size(M::AbstractDFT) = M.N, M.N
 length(M::AbstractDFT) = M.N^2
 
-getindex(M::DFT, m::Int, n::Int) =
-    exp(-2π * im * (m-1) * (n-1) / M.N)
-getindex(N::IDFT, m::Int, n::Int) =
-    exp(2π * im * (m-1) * (n-1) / M.N) / M.N
-getindex(M::UnitaryDFT, m::Int, n::Int) =
-    exp(-2π * im * (m-1) * (n-1) / M.N) / sqrt(M.N)
-getindex(M::UnitaryIDFT, m::Int, n::Int) =
-    exp(2π * im * (m-1) * (n-1) / M.N) / sqrt(M.N)
-getindex(M::AbstractDFT, m::Int) = getindex(M, rem(m, M.N), div(m, M.N))
+# getindex(M::DFT, m::Int, n::Int) =
+#     exp(-2π * im * (m-1) * (n-1) / M.N)
+# getindex(M::UnitaryDFT, m::Int, n::Int) =
+#     exp(-2π * im * (m-1) * (n-1) / M.N) / M.sqrtN
+# getindex(M::AbstractDFT, m::Int) = getindex(M, rem(m, M.N), div(m, M.N))
 
-isassigned(M::AbstractDFT, i::Int) = 1 <= i <= M.N ? true : false
+# isassigned(M::AbstractDFT, i::Int) = 1 <= i <= M.N ? true : false
 
 # Helper functions to check whether stuff is conformal or not.
-function conformal(U::AbstractDFT, x::Vector)
+function conf(U::AbstractDFT, x::SV)
     if U.N != length(x)
-        throw(ArgumentError("U and x are not conformal."))
+        throw(ArgumentError("U and x are not conformal (in this order)."))
     end
 end
-function conformal(U::AbstractDFT, X::Matrix)
+# Always throw an error as it does not make sense to multiply a column Vector
+# by a matrix in this order.
+function conf(x::SV, U::AbstractDFT)
+    throw(ArgumentError("x and U are not conformal (in this order)."))
+end
+function conf(U::AbstractDFT, X::SM)
     if U.N != size(X, 1)
-        throw(ArgumentError("U and X are not conformal."))
+        throw(ArgumentError("U and X are not conformal (in this order)."))
     end
 end
-function conformal(X::Matrix, U::AbstractDFT)
+function conf(X::SM, U::AbstractDFT)
     if size(X, 2) != U.N
-        throw(ArgumentError("X and U are not conformal."))
+        throw(ArgumentError("X and U are not conformal (in this order)."))
     end
 end
 
-# Multiplication from the left by the DFT matrices.
-*(U::DFT, x::Vector) = (conformal(U, x); fft(x))
-*(U::IDFT, x::Vector) = (conformal(U, x); ifft(x))
-*(U::UnitaryDFT, x::Vector) = (conformal(U, x); fft(x) / U.sqrtN)
-*(U::UnitaryIDFT, x::Vector) = (conformal(U, x); bfft(x) / U.sqrtN)
-
-*(U::DFT, X::Matrix) = (conformal(U, X); fft(X, 1))
-*(U::IDFT, X::Matrix) = (conformal(U, X); ifft(X, 1))
-*(U::UnitaryDFT, X::Matrix) = (conformal(U, X); fft(X, 1) / U.sqrtN)
-*(U::UnitaryIDFT, X::Matrix) = (conformal(U, X);  bfft(X, 1) / U.sqrtN)
-
-# Multiplication from the right by the DFT matrices.
-*(X::Matrix, U::DFT) = (conformal(X, U); fft(X, 2))
-*(X::Matrix, U::IDFT) = (conformal(X, U); ifft(X, 2))
-*(X::Matrix, U::UnitaryDFT) = (conformal(X, U); fft(X, 2) / U.sqrtN)
-*(X::Matrix, U::UnitaryIDFT) = (conformal(X, U); bfft(X, 2) / U.sqrtN)
-
-# In-place functionality.
-A_mul_B!(y::Vector, U::DFT, x::Vector) =
-    (conformal(U, x); copy!(y, x); fft!(y))
-A_mul_B!(y::Vector, U::IDFT, x::Vector) =
-    (conformal(U, x); copy!(y, x); ifft!(y))
-A_mul_B!(y::Vector, U::UnitaryDFT, x::Vector) =
-    (conformal(U, x); copy!(y, x); fft!(y); y ./= U.sqrtN)
-A_mul_B!(y::Vector, U::UnitaryIDFT, x::Vector) =
-    (conformal(U, x); copy!(y, x); bfft!(y); y ./= U.sqrtN)
-
-A_mul_B!(Y::Matrix, U::DFT, X::Matrix) =
-    (conformal(U, X); copy!(Y, X); fft!(Y, 1))
-A_mul_B!(Y::Matrix, U::IDFT, X::Matrix) =
-    (conformal(U, X); copy!(Y, X); ifft!(Y, 1))
-A_mul_B!(Y::Matrix, U::UnitaryDFT, X::Matrix) =
-    (conformal(U, X); copy!(Y, X); fft!(Y, 1); y ./= U.sqrtN)
-A_mul_B!(Y::Matrix, U::UnitaryIDFT, X::Matrix) =
-    (conformal(U, X); copy!(Y, X); bfft!(Y, 1); y ./= U.sqrtN)
-
-A_mul_B!(Y::Matrix, X::Matrix, U::DFT) =
-    (conformal(X, U); copy!(Y, X); fft!(Y, 2))
-A_mul_B!(Y::Matrix, X::Matrix, U::IDFT) =
-    (conformal(X, U); copy!(Y, X); ifft!(Y, 2))
-A_mul_B!(Y::Matrix, X::Matrix, U::UnitaryDFT) =
-    (conformal(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
-A_mul_B!(Y::Matrix, X::Matrix, U::UnitaryIDFT) =
-    (conformal(X, U); copy!(Y, X); bfft!(Y, 2); y ./= U.sqrtN)
-
-# Conjugate transpose of DFT times X.
-Ac_mul_B(U::DFT, X::Matrix) = bfft(X, 1)
-Ac_mul_B(U::IDFT, X::Matrix) = fft(X, 1) / U.N
-Ac_mul_B(U::UnitaryDFT, X::Matrix) = bfft(X, 1) / U.sqrtN
-Ac_mul_B(U::UnitaryIDFT, X::Matrix) = fft(X, 1) / U.sqrtN
+# Multiplication from the left and right by the DFT matrices.
+*(U::DFT, X::SVM) = (conf(U, X); fft(X, 1))
+*(X::SVM, U::DFT) = (conf(X, U); fft(X, 2))
+*(U::UnitaryDFT, X::SVM) = (conf(U, X); fft(X, 1) / U.sqrtN)
+*(X::SVM, U::UnitaryDFT) = (conf(X, U); fft(X, 2) / U.sqrtN)
+
+# In-place functionality. Note that diagonal dispatch is used as both arguments should be of the same type.
+A_mul_B!(Y::SV, U::DFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y))
+A_mul_B!(Y::SM, U::DFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1))
+A_mul_B!(Y::SV, X::SV, U::DFT) = (conf(X, U); copy!(Y, X); fft!(Y))
+A_mul_B!(Y::SM, X::SM, U::DFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2))
+A_mul_B!(Y::SV, U::UnitaryDFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y); y ./= U.sqrtN)
+A_mul_B!(Y::SM, U::UnitaryDFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1); y ./= U.sqrtN)
+A_mul_B!(Y::SV, X::SV, U::UnitaryDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
+A_mul_B!(Y::SM, X::SM, U::UnitaryDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
+
+# Conjugate transpose of U'X and X*U'.
+Ac_mul_B(U::DFT, X::SVM) = (conf(U, X); bfft(X, 1))
+A_mul_Bc(X::SVM, U::DFT) = (conf(X, U); bfft(X, 2))
+Ac_mul_B(U::UnitaryDFT, X::SVM) = (conf(U, X); bfft(X, 1) / U.sqrtN)
+A_mul_Bc(X::SVM, U::UnitaryDFT) = (conf(X, U); bfft(X, 2) / U.sqrtN)
+
+# Explicit implementation of X'U and U*X' to avoid fallback in Base being called.
+Ac_mul_B(X::SVM, U::AbstractDFT) = (Xc = ctranspose(X); conf(Xc, U); A_mul_B!(Xc, Xc, U))
+A_mul_Bc(U::AbstractDFT, X::SVM) = (Xc = ctranspose(X); conf(U, Xc); A_mul_B!(Xc, U, Xc))
 
 # All DFT matrices are symmetric.
-At_mul_B(U::AbstractDFT, X::Matrix) = U * X
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+At_mul_B(U::AbstractDFT, X::SVM) = U * X
+A_mul_Bt(X::SVM, U::AbstractDFT) = X * U
diff --git a/src/toeplitz.jl b/src/toeplitz.jl
index e229150..6f23b2d 100644
--- a/src/toeplitz.jl
+++ b/src/toeplitz.jl
@@ -1,4 +1,4 @@
-export Circulant, Toeplitz
+export Toeplitz
 
 immutable Toeplitz{T} <: AbstractArray{T, 2}
 	c :: Vector{T}
@@ -38,84 +38,3 @@ function full{T}(To::Toeplitz{T})
 	end
 	M
 end
-
-
-
-immutable Circulant{T} <: AbstractMatrix{T}
-	c::Vector{T}
-end
-
-typealias Circ Circulant
-
-getindex(C::Circ, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
-isassigned(C::Circ, i::Int, j::Int) = isassigned(C.c,mod(i-j,length(C.c))+1)
-size(C::Circ, r::Int) = (r==1 || r==2) ? length(C.c) :
-    throw(ArgumentError("Invalid dimension $r"))
-size(C::Circ) = size(C,1), size(C,2)
-
-# Generic code for multiplication by a Circulant matrix. Dispatch is used later
-# on to ensure that the types returned are the correct ones. The idea is that
-# if both inputs are Integer, then the output should also be integer.
-# Similarly, if the inputs are both real, then the output should also be real.
-function circ_dot(C::Circ, X::VecOrMat)
-    if size(X, 1) != size(C.c, 1)
-        throw(ArgumentError("C and X are not conformal."))
-    end
-    U = DFT(size(X, 1))
-    return U * Diagonal(U * C.c) * U'X
-end
-function circ_dot(X::VecOrMat, C::Circ)
-    if size(X, 2) != size(C.c, 1)
-        throw(ArgumentError("C and X are not conformal."))
-    end
-    U = DFT(size(X, 1))
-    return X * U * Diagonal(U * C.c) * U'
-end
-
-# In general, the return type of the multiplication operation should be a
-# complex float.
-*(C::Circ, X::VecOrMat) = circ_dot(C, X)
-*(X::VecOrMat, C::Circ) = circ_dot(X, C)
-
-# If both arguments contain reals, then the result should really contain reals.
-*{T<:Real, V<:Real}(C::Circ{T}, X::VecOrMat{V}) = real(circ_dot(C, X))
-*{T<:Real, V<:Real}(X::VecOrMat{V}, C::Circ{T}) = real(circ_dot(X, C))
-
-# If both arguments contain (real) integers, then the result should be both
-# Integer-valued and real. Promote to common type for the purposes of
-# converting back to integer later.
-*{T<:Integer, V<:Integer}(C::Circ{T}, X::VecOrMat{V}) = *(promote(C, X)...)
-*{T<:Integer, V<:Integer}(X::VecOrMat{V}, C::Circ{T}) = *(promote(C, X)...)
-*{T<:Integer}(C::Circ{T}, X::VecOrMat{T}) =
-    convert(T, round.(real.(circ_dot(C, X))))
-*{T<:Integer}(X::VecOrMat{T}, C::Circ{T}) =
-    convert(T, round.(real.(circ_dot(C, X))))
-
-# TODO: Create functionality to allow complex integers to be appropriately handled.
-
-
-
-
-function A_mul_B!{T}(y::StridedVector{T},C::Circulant{T},x::StridedVector{T})
-    xt=fft(x)
-    vt=fft(C.c)
-    yt=ifft(vt.*xt)
-    if T<: Int
-        map!(round,y,yt) 
-    elseif T<: Real
-        map!(real,y,yt)
-    else
-        copy!(y,yt)
-    end
-    return y
-end
-
-function full{T}(C::Circulant{T})
-	n=size(C, 1)
-	M=Array(T, n, n)
-	for i=1:n
-		M[i:n,i] = C.c[1:n-i+1]
-		M[1:i-1,i] = C.c[n-i+2:n]
-	end
-	M
-end
diff --git a/test/circulant.jl b/test/circulant.jl
new file mode 100644
index 0000000..b994dfd
--- /dev/null
+++ b/test/circulant.jl
@@ -0,0 +1,81 @@
+c = [1, 2 + 3.1im, 3, 4 - 2.5im]
+C = Circulant(c)
+x = convert(Vector{Complex{Float64}}, randn(size(C, 1)))
+
+# Important to check this as everything else depends upon it working correctly.
+function check_full(C=C)
+    Ĉ = [C.c circshift(C.c, 1) circshift(C.c, 2) circshift(C.c, 3)]
+    all(C == Ĉ)
+end
+@test check_full()
+
+# Check that, upon copying, the correct semantics are respected.
+function check_copy(C=C)
+    D = copy(C)
+    D.c[1] = 5
+    C.c[1] == 5
+end
+@test check_copy()
+
+# Memory should not be shared between the result and the original.
+function check_conj(C=C)
+    D = conj(C)
+    all(conj(full(C)) == full(conj(C))) && !all(full(D) == full(C))
+end
+@test check_conj()
+
+# Memory should be shared between the new and the original.
+function check_conj!(C=C)
+    copy_C = deepcopy(C)
+    D = conj!(copy_C)
+    all(conj(full(C)) == full(D)) && all(full(D) == full(copy_C))
+end
+@test check_conj!()
+
+function check_transpose(C=C)
+    all(transpose(full(C)) == full(transpose(C)))
+end
+@test check_transpose()
+
+function check_ctranspose()
+    all(ctranspose(full(C)) == full(ctranspose(C)))
+end
+@test check_ctranspose()
+
+function check_add(C=C)
+    all((full(C) + full(C)) == full(C + C)) &&
+        all((full(C) + 5) == full(C + 5)) &&
+        all((5 + full(C)) == full(5 + C))
+end
+@test check_add()
+
+function check_minus(C=C)
+    all(full(-C) == -full(C)) &&
+        all((full(C) - full(C)) == full(C - C)) &&
+        all((full(C) - 5) == full(C - 5))
+end
+@test check_minus()
+
+function check_eigen(C=C, x=x)
+    γ, U = eig(C)
+    abs(sum(Ac_mul_B(U, (Diagonal(γ) * (U * x))) - full(C) * x)) < 1e-9
+end
+@test check_eigen()
+
+check_left_dot(C=C, x=x) = abs(sum(C * x - full(C) * x)) < 1e-9
+@test check_left_dot()
+
+function check_right_dot(C=C, x=x)
+    xt = transpose(x)
+    abs(sum(xt * C - xt * full(C))) < 1e-9
+end
+@test check_right_dot()
+
+
+
+
+
+
+
+
+
diff --git a/test/dft.jl b/test/dft.jl
index 8148450..312a146 100644
--- a/test/dft.jl
+++ b/test/dft.jl
@@ -1,7 +1,19 @@
-function check_vector()
-    N = 5
-    x = randn(N)
-    U = DFT(N)
-    all((U * x) .== (U * reshape(x, N, 1)))
+N = 5
+x = randn(N)
+X = reshape(x, N, 1)
+U = DFT{Float64}(N)
+
+function check_dot(N=N, x=x, X=X, U=U)
+    all((U * x) .== (U * X))
 end
-@test check_vector()
\ No newline at end of file
+@test check_dot()
+
+function check_inplace_dot(N=N, x=x, X=X, U=U)
+    x = convert(Vector{Complex{Float64}}, x)
+    X = convert(Matrix{Complex{Float64}}, X)
+    y = zeros(Complex{Float64}, N)
+    Y = zeros(Complex{Float64}, N, 1)
+    all(A_mul_B!(y, U, x) == U * x) && all(A_mul_B!(Y, U, X) == U * X) &&
+        all(A_mul_B!(Y', X', U) == X'U)
+end
+@test check_inplace_dot()
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index f59515c..5003007 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,11 +1,12 @@
 using SpecialMatrices
 using Base.Test  
 
-include("dft.jl") # Discrete Fourier Transform matrix.
+include("circulant.jl")
 include("companion.jl") #Companion matrix
+include("dft.jl") # Discrete Fourier Transform matrix.
 include("frobenius.jl") #Frobenius matrix
-include("strang.jl") #Strang matrix
 include("hankel.jl") #Hankel matrix
 include("hilbert.jl") #Hilbert matrix
+include("strang.jl") #Strang matrix
 include("toeplitz.jl") #Toeplitz matrix
-include("vandermonde.jl") #Vandermonde matrix
+include("vandermonde.jl") #Vandermonde matrix
\ No newline at end of file

From f449a93f3079ccf1e1ab4025b7463e2fdab7b5d0 Mon Sep 17 00:00:00 2001
From: Will Tebbutt <will.tebbutt@invenialabs.co.uk>
Date: Mon, 16 Jan 2017 18:04:00 +0000
Subject: [PATCH 04/10] Commit to bring remote up to speed.

---
 src/SpecialMatrices.jl |  5 ++-
 src/circulant.jl       | 98 +++++++++++++++++++++++++++++-------------
 src/dft.jl             | 29 +++++++------
 src/hankel.jl          |  3 +-
 test/circulant.jl      | 33 ++++++++------
 test/dft.jl            |  2 +-
 6 files changed, 108 insertions(+), 62 deletions(-)

diff --git a/src/SpecialMatrices.jl b/src/SpecialMatrices.jl
index 8683de5..68f7c08 100644
--- a/src/SpecialMatrices.jl
+++ b/src/SpecialMatrices.jl
@@ -1,16 +1,17 @@
 module SpecialMatrices
 
 import Base: A_mul_B!, full, getindex, inv, isassigned, size, *, length, +, -,
-    Ac_mul_B, A_mul_Bc, At_mul_B, A_mul_Bt, eig
+    Ac_mul_B, A_mul_Bc, At_mul_B, A_mul_Bt, eig, eigfact, Eigen, \, /
 
 typealias SV StridedVector
 typealias SM StridedMatrix
 typealias SVM StridedVecOrMat
 
+include("dft.jl") # Discrete Fourier Transform matrices.
+
 include("cauchy.jl") # Cauchy matrix
 include("circulant.jl") # Circulant matrix.
 include("companion.jl") # Companion matrix
-include("dft.jl") # Discrete Fourier Transform matrices.
 include("frobenius.jl") # Frobenius matrix
 include("hankel.jl") # Hankel matrix
 include("hilbert.jl") # Hilbert matrix
diff --git a/src/circulant.jl b/src/circulant.jl
index f3d0e53..94e7725 100644
--- a/src/circulant.jl
+++ b/src/circulant.jl
@@ -1,4 +1,4 @@
-export Circulant
+export Circulant, CircEig
 
 immutable Circulant{T} <: AbstractMatrix{T}
     c::Vector{T}
@@ -7,7 +7,7 @@ typealias Circ Circulant
 
 function full{T}(C::Circ{T})
     n = size(C, 1)
-    M = Array(T, n, n)
+    M = Array(T, n, n) 
     for i=1:n
         M[i:n,i] = C.c[1:n-i+1]
         M[1:i-1,i] = C.c[n-i+2:n]
@@ -36,43 +36,79 @@ Base.ctranspose(C::Circ) = conj!(transpose(C))
 -(C::Circ, D::Circ) = Circ(C.c - D.c)
 -(C::Circ) = Circ(-C.c)
 
-# # Return the eigen-factorisation of C. Constituents are efficient to compute
-# # and the eigenvectors are represented efficiently using a UnitaryDFT object.
-# function Base.eigfact{T}(C::Circ{T})
-#     n = length(C.c)
-#     Base.LinAlg.Eigen(DFT{T}(n) * C.c, UnitaryDFT{T}(n))
-# end
+# Return the eigen-factorisation of C. Constituents are efficient to compute
+# and the eigenvectors are represented efficiently using a UnitaryDFT object.
+eig{T<:Real}(C::Circ{T}) = (N = length(C.c); (DFT{Complex{T}}(N) * C.c, UDFT{Complex{T}}(N)))
+eig{T<:Complex}(C::Circ{T}) = (N = length(C.c); (DFT{T}(N) * C.c, UDFT{T}(N)))
+eigfact(C::Circ) = Eigen(eig(C)...)
+typealias CircEig{T<:Complex} Eigen{T, T, UDFT{T}, Vector{T}}
 
-eig{T}(C::Circ{T}) = (N = length(C.c); (DFT{T}(N) * C.c, UnitaryDFT{T}(N)))
-
-# Generic code for multiplication by a Circulant matrix. Dispatch is used later
-# on to ensure that the types returned are the correct ones. The idea is that
-# if both inputs are Integer, then the output should also be integer.
-# Similarly, if the inputs are both real, then the output should also be real.
-function circ_dot{T}(C::Circ{T}, X::StridedVecOrMat)
-    if size(X, 1) != size(C, 2)
+# Check for conformal arguments.
+function conf{T}(C::CircEig{T}, X::SVM)
+    if size(X, 1) != size(C.vectors, 2)
         throw(ArgumentError("C and X are not conformal."))
     end
-    γ, U = eig(C)
-    return Ac_mul_B(U, (Diagonal(γ) * (U * X)))
 end
-function circ_dot(X::StridedVecOrMat, C::Circ)
-    if size(X, 2) != size(C, 1)
-        println(size(X, 2))
-        println(size(C, 1))
-        throw(ArgumentError("C and X are not conformal."))
+function conf{T}(X::SVM, C::CircEig{T})
+    if size(X, 2) != size(C.vectors, 1)
+        throw(ArgumentError("X and C are not conformal."))
     end
-    γ, U = eig(C)
-    return (A_mul_Bc(X, U) * Diagonal(γ)) * U
 end
 
-# In general, the return type of the * operation should be a complex float.
-*(C::Circ, X::SVM) = circ_dot(C, X)
-*(X::SVM, C::Circ) = circ_dot(X, C)
+# TODO: How do I correctly do the in-place bits of this? Do some profiling
+# before moving forwards with more functionality.
+
+# Perform the multiplication A * B in place, overwriting B.
+# TODO: TEST THIS!
+function A_mul_B!(A::Diagonal, B::SVM)
+    M, N = size(B)
+    if size(A, 2) != M
+        throw(ArgumentError("A and B are not conformal."))
+    end
+    d = D.diag
+    @show M, N
+    for j in 1:N
+        for i in 1:M
+            B[i, j] *= d[i]
+        end
+    end
+end
 
-# If both arguments contain reals, then the result should really contain reals.
-*{T<:Real, V<:Real}(C::Circ{T}, X::SVM{V}) = real!.(circ_dot(C, X))
-*{T<:Real, V<:Real}(X::SVM{V}, C::Circ{T}) = real!.(circ_dot(X, C))
+# Always compute the eigen-factorisation to compute the multiplication. If
+# this operation is to be performed multiple times, then we should just cache
+# the factorisation to save re-computing it for each multiplication operation.
+function *{T}(C::CircEig{T}, X::SVM)
+    conf(C, X)
+    Γ, U = Diagonal(C.values), C.vectors
+    Ac_mul_B(U, (Γ * (U * X)))
+    tmp = U * X
+    bfft!(A_mul_B!(tmp, Γ, tmp), 1) / U.sqrtN
+    @show Ac_mul_B(U, (Γ * (U * X)))
+    tmp = U * X
+    @show bfft!(A_mul_B!(Γ, tmp), 1) / U.sqrtN
+end
+function *{T}(X::SVM, C::CircEig{T})
+    conf(X, C)
+    Γ, U = Diagonal(C.values), C.vectors
+    (A_mul_Bc(X, U) * Γ) * U
+end
+*(C::Circ, X::SVM) = eigfact(C) * X
+*(X::SVM, C::Circ) = X * eigfact(C)
+
+# Use the eigen-factorisation to compute the inv(C) * X. The inverse
+# factorisation should be cached if the operation is required multiple times.
+function \{T}(C::CircEig{T}, X::SVM)
+    conf(C, X)
+    Γ, U = Diagonal(C.values), C.vectors
+    Ac_mul_B(U, (Γ \ (U * X)))
+end
+function /{T}(X::SVM, C::CircEig{T})
+    conf(X, C)
+    Γ, U = Diagonal(C.values), C.vectors
+    (A_mul_Bc(X, U) / Γ) * U
+end
+\(C::Circ, X::SVM) = eigfact(C) \ X
+/(X::SVM, C::Circ) = X / eigfact(C)
 
 # TODO: How attached are people to this? This should probably change. Also, this doesn't
 # really get you anything in terms of performance. There's not reduction in memory use as far as I can see
diff --git a/src/dft.jl b/src/dft.jl
index 945ecec..4fb3248 100644
--- a/src/dft.jl
+++ b/src/dft.jl
@@ -1,27 +1,28 @@
-export DFT, UnitaryDFT
+export DFT, UnitaryDFT, UDFT
 
 # Non-unitary DFT matrix.
 abstract AbstractDFT{T} <: AbstractMatrix{T}
-immutable DFT{T} <: AbstractDFT{T}
+immutable DFT{T<:Complex} <: AbstractDFT{T}
     N::Int
 end
 
 # Unitary DFT Matrix.
-immutable UnitaryDFT{T} <: AbstractDFT{T}
+immutable UnitaryDFT{T<:Complex} <: AbstractDFT{T}
     N::Int
     sqrtN::Float64
     UnitaryDFT(N) = new(N, sqrt(N))
 end
+typealias UDFT UnitaryDFT
 
-Base.show(io::IO, U::DFT) = print(io, "DFT matrix.")
-Base.show(io::IO, U::UnitaryDFT) = print(io, "Unitary DFT matrix.")
+Base.show(io::IO, U::DFT) = print(io, "DFT matrix of size $U.N.")
+Base.show(io::IO, U::UDFT) = print(io, "Unitary DFT matrix of size $U.N.")
 
 size(M::AbstractDFT) = M.N, M.N
 length(M::AbstractDFT) = M.N^2
 
 # getindex(M::DFT, m::Int, n::Int) =
 #     exp(-2π * im * (m-1) * (n-1) / M.N)
-# getindex(M::UnitaryDFT, m::Int, n::Int) =
+# getindex(M::UnitaryDFT, m::Int, n::Int) =1
 #     exp(-2π * im * (m-1) * (n-1) / M.N) / M.sqrtN
 # getindex(M::AbstractDFT, m::Int) = getindex(M, rem(m, M.N), div(m, M.N))
 
@@ -52,24 +53,24 @@ end
 # Multiplication from the left and right by the DFT matrices.
 *(U::DFT, X::SVM) = (conf(U, X); fft(X, 1))
 *(X::SVM, U::DFT) = (conf(X, U); fft(X, 2))
-*(U::UnitaryDFT, X::SVM) = (conf(U, X); fft(X, 1) / U.sqrtN)
-*(X::SVM, U::UnitaryDFT) = (conf(X, U); fft(X, 2) / U.sqrtN)
+*(U::UDFT, X::SVM) = (conf(U, X); fft(X, 1) / U.sqrtN)
+*(X::SVM, U::UDFT) = (conf(X, U); fft(X, 2) / U.sqrtN)
 
 # In-place functionality. Note that diagonal dispatch is used as both arguments should be of the same type.
 A_mul_B!(Y::SV, U::DFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y))
 A_mul_B!(Y::SM, U::DFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1))
 A_mul_B!(Y::SV, X::SV, U::DFT) = (conf(X, U); copy!(Y, X); fft!(Y))
 A_mul_B!(Y::SM, X::SM, U::DFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2))
-A_mul_B!(Y::SV, U::UnitaryDFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y); y ./= U.sqrtN)
-A_mul_B!(Y::SM, U::UnitaryDFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1); y ./= U.sqrtN)
-A_mul_B!(Y::SV, X::SV, U::UnitaryDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
-A_mul_B!(Y::SM, X::SM, U::UnitaryDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
+A_mul_B!(Y::SV, U::UDFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y); y ./= U.sqrtN)
+A_mul_B!(Y::SM, U::UDFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1); y ./= U.sqrtN)
+A_mul_B!(Y::SV, X::SV, U::UDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
+A_mul_B!(Y::SM, X::SM, U::UDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
 
 # Conjugate transpose of U'X and X*U'.
 Ac_mul_B(U::DFT, X::SVM) = (conf(U, X); bfft(X, 1))
 A_mul_Bc(X::SVM, U::DFT) = (conf(X, U); bfft(X, 2))
-Ac_mul_B(U::UnitaryDFT, X::SVM) = (conf(U, X); bfft(X, 1) / U.sqrtN)
-A_mul_Bc(X::SVM, U::UnitaryDFT) = (conf(X, U); bfft(X, 2) / U.sqrtN)
+Ac_mul_B(U::UDFT, X::SVM) = (conf(U, X); bfft(X, 1) / U.sqrtN)
+A_mul_Bc(X::SVM, U::UDFT) = (conf(X, U); bfft(X, 2) / U.sqrtN)
 
 # Explicit implementation of X'U and U*X' to avoid fallback in Base being called.
 Ac_mul_B(X::SVM, U::AbstractDFT) = (Xc = ctranspose(X); conf(Xc, U); A_mul_B!(Xc, Xc, U))
diff --git a/src/hankel.jl b/src/hankel.jl
index 8fd2630..f198dd3 100644
--- a/src/hankel.jl
+++ b/src/hankel.jl
@@ -3,7 +3,8 @@ export Hankel
 immutable Hankel{T} <: AbstractArray{T, 2}
     c :: Vector{T}
 end
-Hankel{T}(c::Vector{T}) = length(c) % 2 == 1 ? Hankel{T}(c) : throw(ArgumentError(""))
+(::Hankel{T}){T}(c::Vector{T}) = length(c) % 2 == 1 ? Hankel{T}(c) : throw(ArgumentError(""))
+
 
 getindex(H::Hankel, i::Int, j::Int) = H.c[i+j-1]
 isassigned(H::Hankel, i::Int, j::Int) = isassigned(H.c, i+j-1)
diff --git a/test/circulant.jl b/test/circulant.jl
index b994dfd..9a79a82 100644
--- a/test/circulant.jl
+++ b/test/circulant.jl
@@ -1,6 +1,8 @@
 c = [1, 2 + 3.1im, 3, 4 - 2.5im]
 C = Circulant(c)
 x = convert(Vector{Complex{Float64}}, randn(size(C, 1)))
+xt = transpose(x)
+Cr = Circulant(randn(4))
 
 # Important to check this as everything else depends upon it working correctly.
 function check_full(C=C)
@@ -57,25 +59,30 @@ end
 @test check_minus()
 
 function check_eigen(C=C, x=x)
-    γ, U = eig(C)
-    abs(sum(Ac_mul_B(U, (Diagonal(γ) * (U * x))) - full(C) * x)) < 1e-9
+    eigen = eigfact(C)
+    Γ = Diagonal(eigen.values)
+    U = eigen.vectors
+    abs(sum(Ac_mul_B(U, (Γ * (U * x))) - full(C) * x)) < 1e-9
 end
-@test check_eigen()
+@test check_eigen(Cr)
 
 check_left_dot(C=C, x=x) = abs(sum(C * x - full(C) * x)) < 1e-9
 @test check_left_dot()
 
-function check_right_dot(C=C, x=x)
-    xt = transpose(x)
-    abs(sum(xt * C - xt * full(C))) < 1e-9
-end
+check_right_dot(C=C, xt=xt) = abs(sum(xt * C - xt * full(C))) < 1e-9
 @test check_right_dot()
 
+check_backslash(C=C, x=x) = abs(sum(C * (C \ x) - x)) < 1e-9
+@test check_backslash()
 
+check_forwardslash(C=C, xt=xt) = abs(sum((xt / C) * C - xt)) < 1e-9
+@test check_forwardslash()
 
-
-
-
-
-
-
+function check_inplace_diag_mult()
+    x = randn(100)
+    d = randn(100)
+    D = Diagonal(d)
+    tmp = D * x
+    all(tmp == A_mul_B!(D, x))
+end
+check_inplace_diag_mult()
\ No newline at end of file
diff --git a/test/dft.jl b/test/dft.jl
index 312a146..75d91d8 100644
--- a/test/dft.jl
+++ b/test/dft.jl
@@ -1,7 +1,7 @@
 N = 5
 x = randn(N)
 X = reshape(x, N, 1)
-U = DFT{Float64}(N)
+U = DFT{Complex128}(N)
 
 function check_dot(N=N, x=x, X=X, U=U)
     all((U * x) .== (U * X))

From da86f31e871101e9e9795ea1d05a47d3fb39af44 Mon Sep 17 00:00:00 2001
From: Will Tebbutt <will.tebbutt@invenialabs.co.uk>
Date: Thu, 19 Jan 2017 17:07:45 +0000
Subject: [PATCH 05/10] Circulant functionality sufficiently complete for a
 first vesion.

---
 src/SpecialMatrices.jl |  4 +-
 src/circulant.jl       | 90 ++++++++++++++++++++----------------------
 src/dft.jl             | 63 +++++++++++++++--------------
 test/circulant.jl      | 69 +++++++++++++++++++++++++-------
 4 files changed, 133 insertions(+), 93 deletions(-)

diff --git a/src/SpecialMatrices.jl b/src/SpecialMatrices.jl
index 68f7c08..4025cb8 100644
--- a/src/SpecialMatrices.jl
+++ b/src/SpecialMatrices.jl
@@ -1,7 +1,9 @@
 module SpecialMatrices
 
 import Base: A_mul_B!, full, getindex, inv, isassigned, size, *, length, +, -,
-    Ac_mul_B, A_mul_Bc, At_mul_B, A_mul_Bt, eig, eigfact, Eigen, \, /
+    Ac_mul_B, A_mul_Bc, At_mul_B, A_mul_Bt, eig, eigfact, Eigen, \, /,
+    transpose, ctranspose, copy, conj, conj!, inv, det, logdet, real, convert,
+    round
 
 typealias SV StridedVector
 typealias SM StridedMatrix
diff --git a/src/circulant.jl b/src/circulant.jl
index 94e7725..bc09f20 100644
--- a/src/circulant.jl
+++ b/src/circulant.jl
@@ -1,4 +1,4 @@
-export Circulant, CircEig
+export Circulant, CircEig, tocirc
 
 immutable Circulant{T} <: AbstractMatrix{T}
     c::Vector{T}
@@ -15,18 +15,18 @@ function full{T}(C::Circ{T})
     M
 end
 
-getindex(C::Circ, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
+# getindex(C::Circ, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
 isassigned(C::Circ, i::Int, j::Int) = isassigned(C.c,mod(i-j,length(C.c))+1)
 size(C::Circ, r::Int) = (r==1 || r==2) ? length(C.c) :
     throw(ArgumentError("Invalid dimension $r"))
 size(C::Circ) = size(C,1), size(C,2)
 
-Base.copy(C::Circ) = Circ(C.c)
-Base.conj(C::Circ) = Circ(conj(C.c))
-Base.conj!(C::Circ) = (conj!(C.c); C)
+copy(C::Circ) = Circ(C.c)
+conj(C::Circ) = Circ(conj(C.c))
+conj!(C::Circ) = (conj!(C.c); C)
 
-Base.transpose(C::Circ) = Circ(circshift(reverse(C.c), 1))
-Base.ctranspose(C::Circ) = conj!(transpose(C))
+transpose(C::Circ) = Circ(circshift(reverse(C.c), 1))
+ctranspose(C::Circ) = conj!(transpose(C))
 
 +(C::Circ, a::Number) = Circ(C.c + a)
 +(a::Number, C::Circ) = Circ(a + C.c)
@@ -43,6 +43,9 @@ eig{T<:Complex}(C::Circ{T}) = (N = length(C.c); (DFT{T}(N) * C.c, UDFT{T}(N)))
 eigfact(C::Circ) = Eigen(eig(C)...)
 typealias CircEig{T<:Complex} Eigen{T, T, UDFT{T}, Vector{T}}
 
+# Utility function to convert back to usual circulant formulation.
+tocirc{T}(C::CircEig{T}) = (N = length(C.values); Circ(DFT{T}(N)'C.values / N))
+
 # Check for conformal arguments.
 function conf{T}(C::CircEig{T}, X::SVM)
     if size(X, 1) != size(C.vectors, 2)
@@ -55,25 +58,6 @@ function conf{T}(X::SVM, C::CircEig{T})
     end
 end
 
-# TODO: How do I correctly do the in-place bits of this? Do some profiling
-# before moving forwards with more functionality.
-
-# Perform the multiplication A * B in place, overwriting B.
-# TODO: TEST THIS!
-function A_mul_B!(A::Diagonal, B::SVM)
-    M, N = size(B)
-    if size(A, 2) != M
-        throw(ArgumentError("A and B are not conformal."))
-    end
-    d = D.diag
-    @show M, N
-    for j in 1:N
-        for i in 1:M
-            B[i, j] *= d[i]
-        end
-    end
-end
-
 # Always compute the eigen-factorisation to compute the multiplication. If
 # this operation is to be performed multiple times, then we should just cache
 # the factorisation to save re-computing it for each multiplication operation.
@@ -81,20 +65,42 @@ function *{T}(C::CircEig{T}, X::SVM)
     conf(C, X)
     Γ, U = Diagonal(C.values), C.vectors
     Ac_mul_B(U, (Γ * (U * X)))
-    tmp = U * X
-    bfft!(A_mul_B!(tmp, Γ, tmp), 1) / U.sqrtN
-    @show Ac_mul_B(U, (Γ * (U * X)))
-    tmp = U * X
-    @show bfft!(A_mul_B!(Γ, tmp), 1) / U.sqrtN
+    Y = U * X
+    Ac_mul_B!(U, Γ * (U * X))
 end
 function *{T}(X::SVM, C::CircEig{T})
     conf(X, C)
     Γ, U = Diagonal(C.values), C.vectors
     (A_mul_Bc(X, U) * Γ) * U
+    Y = A_mul_Bc(X, U)
+    A_mul_B!(A_mul_B!(Y, Y, Γ), U)
 end
 *(C::Circ, X::SVM) = eigfact(C) * X
 *(X::SVM, C::Circ) = X * eigfact(C)
 
+# "In-place" multiplication.
+function A_mul_B!{T<:Number, V<:Complex}(Y::SVM{V}, C::CircEig{T}, X::SVM{V})
+    conf(C, X) 
+    Γ, U = Diagonal(C.values), C.vectors
+    Ac_mul_B!(U, A_mul_B!(Y, Γ, A_mul_B!(U, X)))
+end
+function A_mul_B!{T<:Number, V<:Complex}(Y::SVM{V}, X::SVM{V}, C::CircEig{T})
+    conf(X, C)
+    Γ, U = Diagonal(C.values), C.vectors
+    A_mul_B!(A_mul_B!(Y, A_mul_Bc!(X, U), Γ), U)
+end
+A_mul_B!(Y::SVM, C::Circ, X::SVM) = A_mul_B!(Y, eigfact(C), X)
+A_mul_B!(Y::SVM, X::SVM, C::Circ)  = A_mul_B!(Y, X, eigfact(C))
+
+inv{T}(C::CircEig{T}) = Eigen(1 ./ C.values, UDFT{T}(length(C.values)))
+inv(C::Circ) = tocirc(inv(eigfact(C)))
+
+det{T}(C::CircEig{T}) = prod(C.values)
+det(C::Circ) = det(eigfact(C))
+
+logdet{T}(C::CircEig{T}) = sum(log(C.values))
+logdet(C::Circ) = logdet(eigfact(C))
+
 # Use the eigen-factorisation to compute the inv(C) * X. The inverse
 # factorisation should be cached if the operation is required multiple times.
 function \{T}(C::CircEig{T}, X::SVM)
@@ -110,21 +116,9 @@ end
 \(C::Circ, X::SVM) = eigfact(C) \ X
 /(X::SVM, C::Circ) = X / eigfact(C)
 
-# TODO: How attached are people to this? This should probably change. Also, this doesn't
-# really get you anything in terms of performance. There's not reduction in memory use as far as I can see
-# as the in-place fourier transform is not used.
-function A_mul_B!{T}(y::StridedVector{T},C::Circulant{T},x::StridedVector{T})
-    xt=fft(x)
-    vt=fft(C.c)
-    yt=ifft(vt.*xt)
-    if T<: Int
-        map!(round,y,yt) 
-    elseif T<: Real
-        map!(real,y,yt)
-    else
-        copy!(y,yt)
-    end
-    return y
-end
+# Helper functionality for casting.
+real(C::Circ) = Circ(real(C.c))
+convert{T}(::Type{Circ{T}}, C::Circ) = Circ(convert(Vector{T}, C.c))
 
-# TODO: Implement all of the other functionality that you would expect to find with a matrix.
+round(C::Circ) = Circ(round(C.c))
+round{T}(::Type{Circ{T}}, C::Circ) = Circ(round(Vector{T}, C.c))
diff --git a/src/dft.jl b/src/dft.jl
index 4fb3248..86005b1 100644
--- a/src/dft.jl
+++ b/src/dft.jl
@@ -20,6 +20,13 @@ Base.show(io::IO, U::UDFT) = print(io, "Unitary DFT matrix of size $U.N.")
 size(M::AbstractDFT) = M.N, M.N
 length(M::AbstractDFT) = M.N^2
 
+# Don't bother to implement an in-place transpose as U is sufficiently
+# light-weight not to care in the slightest.
+transpose(U::AbstractDFT) = copy(U)
+
+# Indexing is explicitely disabled as accidentally using the naive version of
+# the DFT is a really bad idea. Basically any functionality that would allow
+# you to do so has been disabled.
 # getindex(M::DFT, m::Int, n::Int) =
 #     exp(-2π * im * (m-1) * (n-1) / M.N)
 # getindex(M::UnitaryDFT, m::Int, n::Int) =1
@@ -29,26 +36,10 @@ length(M::AbstractDFT) = M.N^2
 # isassigned(M::AbstractDFT, i::Int) = 1 <= i <= M.N ? true : false
 
 # Helper functions to check whether stuff is conformal or not.
-function conf(U::AbstractDFT, x::SV)
-    if U.N != length(x)
-        throw(ArgumentError("U and x are not conformal (in this order)."))
-    end
-end
-# Always throw an error as it does not make sense to multiply a column Vector
-# by a matrix in this order.
-function conf(x::SV, U::AbstractDFT)
-    throw(ArgumentError("x and U are not conformal (in this order)."))
-end
-function conf(U::AbstractDFT, X::SM)
-    if U.N != size(X, 1)
-        throw(ArgumentError("U and X are not conformal (in this order)."))
-    end
-end
-function conf(X::SM, U::AbstractDFT)
-    if size(X, 2) != U.N
-        throw(ArgumentError("X and U are not conformal (in this order)."))
-    end
-end
+conf(U::AbstractDFT, x::SV) = assert(U.N == length(x))
+conf(x::SV, U::AbstractDFT) = assert(false)
+conf(U::AbstractDFT, X::SM) = assert(U.N == size(X, 1))
+conf(X::SM, U::AbstractDFT) = assert(size(X, 2) == U.N)
 
 # Multiplication from the left and right by the DFT matrices.
 *(U::DFT, X::SVM) = (conf(U, X); fft(X, 1))
@@ -56,25 +47,39 @@ end
 *(U::UDFT, X::SVM) = (conf(U, X); fft(X, 1) / U.sqrtN)
 *(X::SVM, U::UDFT) = (conf(X, U); fft(X, 2) / U.sqrtN)
 
-# In-place functionality. Note that diagonal dispatch is used as both arguments should be of the same type.
+# In-place functionality matrix multiplication.
 A_mul_B!(Y::SV, U::DFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y))
-A_mul_B!(Y::SM, U::DFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1))
 A_mul_B!(Y::SV, X::SV, U::DFT) = (conf(X, U); copy!(Y, X); fft!(Y))
+A_mul_B!(Y::SV, U::UDFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y); Y ./= U.sqrtN)
+A_mul_B!(Y::SV, X::SV, U::UDFT) = (conf(X, U); copy!(Y, X); fft!(Y); Y ./= U.sqrtN)
+A_mul_B!(Y::SM, U::DFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1))
 A_mul_B!(Y::SM, X::SM, U::DFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2))
-A_mul_B!(Y::SV, U::UDFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y); y ./= U.sqrtN)
-A_mul_B!(Y::SM, U::UDFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1); y ./= U.sqrtN)
-A_mul_B!(Y::SV, X::SV, U::UDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
-A_mul_B!(Y::SM, X::SM, U::UDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); y ./= U.sqrtN)
+A_mul_B!(Y::SM, U::UDFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1); Y ./= U.sqrtN)
+A_mul_B!(Y::SM, X::SM, U::UDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); Y ./= U.sqrtN)
+A_mul_B!(U::DFT, X::SVM) = (conf(U, X); fft!(X, 1))
+A_mul_B!(X::SVM, U::DFT) = (conf(X, U); fft!(X, 2))
+A_mul_B!(U::UDFT, X::SVM) = (conf(U, X); fft!(X, 1); X ./= U.sqrtN)
+A_mul_B!(X::SVM, U::UDFT) = (conf(X, U); fft!(X, 2); X ./= U.sqrtN)
 
-# Conjugate transpose of U'X and X*U'.
+# Conjugate transpose U'X and X*U'.
 Ac_mul_B(U::DFT, X::SVM) = (conf(U, X); bfft(X, 1))
 A_mul_Bc(X::SVM, U::DFT) = (conf(X, U); bfft(X, 2))
 Ac_mul_B(U::UDFT, X::SVM) = (conf(U, X); bfft(X, 1) / U.sqrtN)
 A_mul_Bc(X::SVM, U::UDFT) = (conf(X, U); bfft(X, 2) / U.sqrtN)
 
+# Conjugate transpose operations in-place.
+Ac_mul_B!{T<:SVM}(Y::T, U::DFT, X::T) = (conf(U, X); copy!(Y, X); bfft!(Y, 1))
+A_mul_Bc!{T<:SVM}(Y::T, X::T, U::DFT) = (conf(X, U); copy!(Y, X); bfft!(Y, 2))
+Ac_mul_B!{T<:SVM}(Y::T, U::UDFT, X::T) = (conf(U, X); copy!(Y, X); bfft!(Y, 1); Y ./= U.sqrtN)
+A_mul_Bc!{T<:SVM}(Y::T, X::T, U::UDFT) = (conf(X, U); copy!(Y, X); bfft!(Y, 2); Y ./= U.sqrtN)
+Ac_mul_B!(U::DFT, X::SVM) = (conf(U, X); bfft!(X, 1))
+A_mul_Bc!(X::SVM, U::DFT) = (conf(X, U); bfft!(X, 2))
+Ac_mul_B!(U::UDFT, X::SVM) = (conf(U, X); bfft!(X, 1); X ./= U.sqrtN)
+A_mul_Bc!(X::SVM, U::UDFT) = (conf(X, U); bfft!(X, 2); X ./= U.sqrtN)
+
 # Explicit implementation of X'U and U*X' to avoid fallback in Base being called.
-Ac_mul_B(X::SVM, U::AbstractDFT) = (Xc = ctranspose(X); conf(Xc, U); A_mul_B!(Xc, Xc, U))
-A_mul_Bc(U::AbstractDFT, X::SVM) = (Xc = ctranspose(X); conf(U, Xc); A_mul_B!(Xc, U, Xc))
+Ac_mul_B(X::SVM, U::AbstractDFT) = (Xc = ctranspose(X); conf(Xc, U); A_mul_B!(Xc, U))
+A_mul_Bc(U::AbstractDFT, X::SVM) = (Xc = ctranspose(X); conf(U, Xc); A_mul_B!(U, Xc))
 
 # All DFT matrices are symmetric.
 At_mul_B(U::AbstractDFT, X::SVM) = U * X
diff --git a/test/circulant.jl b/test/circulant.jl
index 9a79a82..f56fbf1 100644
--- a/test/circulant.jl
+++ b/test/circulant.jl
@@ -1,13 +1,15 @@
-c = [1, 2 + 3.1im, 3, 4 - 2.5im]
+c = [1, 2.5 + 3.1im, 3, 4 - 2.5im]
 C = Circulant(c)
 x = convert(Vector{Complex{Float64}}, randn(size(C, 1)))
 xt = transpose(x)
 Cr = Circulant(randn(4))
 
+ϵ = 1e-12
+
 # Important to check this as everything else depends upon it working correctly.
 function check_full(C=C)
     Ĉ = [C.c circshift(C.c, 1) circshift(C.c, 2) circshift(C.c, 3)]
-    all(C == Ĉ)
+    all(full(C) == Ĉ)
 end
 @test check_full()
 
@@ -47,7 +49,7 @@ end
 function check_add(C=C)
     all((full(C) + full(C)) == full(C + C)) &&
         all((full(C) + 5) == full(C + 5)) &&
-        all((5 + full(C)) == full(5 + C))
+        all((5 + full(C)) == full(5 + C))   
 end
 @test check_add()
 
@@ -62,27 +64,64 @@ function check_eigen(C=C, x=x)
     eigen = eigfact(C)
     Γ = Diagonal(eigen.values)
     U = eigen.vectors
-    abs(sum(Ac_mul_B(U, (Γ * (U * x))) - full(C) * x)) < 1e-9
+    abs(sum(Ac_mul_B(U, (Γ * (U * x))) - full(C) * x)) < ϵ
 end
 @test check_eigen(Cr)
 
-check_left_dot(C=C, x=x) = abs(sum(C * x - full(C) * x)) < 1e-9
+function check_to_circ(C=C)
+    out = tocirc(eigfact(C))
+    abs(sum(full(C) - full(out))) < ϵ
+end
+@test check_to_circ()
+
+check_left_dot(C=C, x=x) = abs(sum(C * x - full(C) * x)) < ϵ
 @test check_left_dot()
 
-check_right_dot(C=C, xt=xt) = abs(sum(xt * C - xt * full(C))) < 1e-9
+check_right_dot(C=C, xt=xt) = abs(sum(xt * C - xt * full(C))) < ϵ
 @test check_right_dot()
 
-check_backslash(C=C, x=x) = abs(sum(C * (C \ x) - x)) < 1e-9
+check_backslash(C=C, x=x) = abs(sum(C * (C \ x) - x)) < ϵ
 @test check_backslash()
 
-check_forwardslash(C=C, xt=xt) = abs(sum((xt / C) * C - xt)) < 1e-9
+check_forwardslash(C=C, xt=xt) = abs(sum((xt / C) * C - xt)) < ϵ
 @test check_forwardslash()
 
-function check_inplace_diag_mult()
-    x = randn(100)
-    d = randn(100)
-    D = Diagonal(d)
-    tmp = D * x
-    all(tmp == A_mul_B!(D, x))
+function check_left_dot!(C=C, x=x)
+    x_loc = deepcopy(x)
+    x_out = deepcopy(x_loc)
+    tmp = C * x_loc
+    A_mul_B!(x_out, C, x_loc)
+    abs(sum(x_out - tmp)) < ϵ
+end
+@test check_left_dot!()
+
+function check_right_dot!(C=C, xt=xt)
+    xt_loc = deepcopy(xt)
+    tmp = xt_loc * C
+    out = A_mul_B!(xt_loc, xt_loc, C)
+    abs(sum(x - tmp)) < ϵ
+end
+
+check_inv(C=C) = abs(sum(inv(full(C)) - full(inv(C)))) < ϵ
+@test check_inv()
+
+check_det(C=C) = abs(det(full(C)) - det(C)) < ϵ
+@test check_det()
+
+check_logdet(C=C) = abs(logdet(full(C)) - logdet(C)) < ϵ
+@test check_det()
+
+check_real(C=C) = abs(sum(real(full(C)) - full(real(C)))) < ϵ
+@test check_real()
+
+function check_round(C=C)
+    abs(sum(round(full(real(C))) - full(round(real(C))))) < ϵ
+end
+@test check_round()
+
+function check_convert(C=C)
+    from_full = convert(Matrix{Int64}, full(round(real(C))))
+    from_circ = full(convert(Circulant{Int64}, real(round(C))))
+    abs(sum(from_full - from_circ)) < ϵ
 end
-check_inplace_diag_mult()
\ No newline at end of file
+@test check_convert()

From 93b767350bad5f85bf4882c74a3fca0c98893d8b Mon Sep 17 00:00:00 2001
From: Will Tebbutt <wt0881@my.bristol.ac.uk>
Date: Fri, 20 Jan 2017 11:40:29 +0000
Subject: [PATCH 06/10] Stupid mistake in circulant

Removed two lines of redundant test code from circulant.jl.
---
 src/circulant.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/circulant.jl b/src/circulant.jl
index bc09f20..0dfcbb0 100644
--- a/src/circulant.jl
+++ b/src/circulant.jl
@@ -64,14 +64,12 @@ end
 function *{T}(C::CircEig{T}, X::SVM)
     conf(C, X)
     Γ, U = Diagonal(C.values), C.vectors
-    Ac_mul_B(U, (Γ * (U * X)))
     Y = U * X
-    Ac_mul_B!(U, Γ * (U * X))
+    Ac_mul_B!(U, Γ * Y)
 end
 function *{T}(X::SVM, C::CircEig{T})
     conf(X, C)
     Γ, U = Diagonal(C.values), C.vectors
-    (A_mul_Bc(X, U) * Γ) * U
     Y = A_mul_Bc(X, U)
     A_mul_B!(A_mul_B!(Y, Y, Γ), U)
 end

From 7eb980f9c785bde0e26cb42a434f83f7c38c4179 Mon Sep 17 00:00:00 2001
From: Will Tebbutt <wt0881@my.bristol.ac.uk>
Date: Fri, 27 Jan 2017 22:10:13 +0000
Subject: [PATCH 07/10] Change over syntax to preferred.

---
 src/circulant.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/circulant.jl b/src/circulant.jl
index 0dfcbb0..0b88ca2 100644
--- a/src/circulant.jl
+++ b/src/circulant.jl
@@ -7,7 +7,7 @@ typealias Circ Circulant
 
 function full{T}(C::Circ{T})
     n = size(C, 1)
-    M = Array(T, n, n) 
+    M = Matrix{T}(n, n) 
     for i=1:n
         M[i:n,i] = C.c[1:n-i+1]
         M[1:i-1,i] = C.c[n-i+2:n]

From 746c4acf48e79613bcf853c4ed17fb302e9c6de4 Mon Sep 17 00:00:00 2001
From: Will Tebbutt <wt0881@my.bristol.ac.uk>
Date: Fri, 27 Jan 2017 22:13:47 +0000
Subject: [PATCH 08/10] Added in-place transpose method.

---
 src/dft.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/dft.jl b/src/dft.jl
index 86005b1..86f87d7 100644
--- a/src/dft.jl
+++ b/src/dft.jl
@@ -23,6 +23,7 @@ length(M::AbstractDFT) = M.N^2
 # Don't bother to implement an in-place transpose as U is sufficiently
 # light-weight not to care in the slightest.
 transpose(U::AbstractDFT) = copy(U)
+transpose!(U::AbstractDFT) = copy(U)
 
 # Indexing is explicitely disabled as accidentally using the naive version of
 # the DFT is a really bad idea. Basically any functionality that would allow

From b092d59821af4e4cc66b725b5ad5fac117c6a696 Mon Sep 17 00:00:00 2001
From: Will Tebbutt <wt0881@my.bristol.ac.uk>
Date: Fri, 27 Jan 2017 22:14:47 +0000
Subject: [PATCH 09/10] Added new-line at end

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 5003007..ecad7ae 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,4 +9,4 @@ include("hankel.jl") #Hankel matrix
 include("hilbert.jl") #Hilbert matrix
 include("strang.jl") #Strang matrix
 include("toeplitz.jl") #Toeplitz matrix
-include("vandermonde.jl") #Vandermonde matrix
\ No newline at end of file
+include("vandermonde.jl") #Vandermonde matrix

From 15beb2a292120fccbb6be2041348d74829a7c7a3 Mon Sep 17 00:00:00 2001
From: Will Tebbutt <wt0881@my.bristol.ac.uk>
Date: Fri, 27 Jan 2017 22:15:26 +0000
Subject: [PATCH 10/10] Added new-line at end of file.

---
 test/dft.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/dft.jl b/test/dft.jl
index 75d91d8..54feca4 100644
--- a/test/dft.jl
+++ b/test/dft.jl
@@ -16,4 +16,4 @@ function check_inplace_dot(N=N, x=x, X=X, U=U)
     all(A_mul_B!(y, U, x) == U * x) && all(A_mul_B!(Y, U, X) == U * X) &&
         all(A_mul_B!(Y', X', U) == X'U)
 end
-@test check_inplace_dot()
\ No newline at end of file
+@test check_inplace_dot()