Simpler CTranspose + multiplication methods

src/arrays/arraymath.jl

@@ -0,0 +1,62 @@
+import Base: *
+
+##################
+# Multiplication #
+##################
+
+# The below needs to be defined more generically;
+# the below only defines multiplication on orthonormal, 
+# N<3 objects 
+
+typealias OrthonormalBasis{S<:Orthonormal} AbstractBasis{S}
+
+# bra * ket -> scalar
+function *{B<:OrthonormalBasis}(bra::DualVector{B}, ket::QuVector{B})
+    return dot(rawcoeffs(bra), rawcoeffs(ket))
+end
+
+# bra * operator -> bra
+function *{B<:OrthonormalBasis}(bra::DualVector{B}, op::QuMatrix{B}) 
+    return QuArray(Ac_mul_B(rawcoeffs(op), rawcoeffs(bra)), bases(op,2))'
+end
+
+function *{B<:OrthonormalBasis}(bra::DualVector{B}, op::DualMatrix{B}) 
+    return QuArray(rawcoeffs(op)*rawcoeffs(bra), bases(op,2))'
+end
+
+# operator * ket -> ket
+function *{B<:OrthonormalBasis}(op::QuMatrix{B}, ket::QuVector{B}) 
+    return QuArray(rawcoeffs(op)*rawcoeffs(ket), bases(op,1))
+end
+
+function *{B<:OrthonormalBasis}(op::DualMatrix{B}, ket::QuVector{B})
+    return QuArray(Ac_mul_B(rawcoeffs(op), rawcoeffs(ket)), bases(op,1))
+end
+
+# ket * bra -> operator
+function *{B<:OrthonormalBasis}(ket::QuVector{B}, bra::DualVector{B}) 
+    return QuArray(A_mul_Bc(rawcoeffs(ket), rawcoeffs(bra)), 
+                   bases(ket,1), 
+                   bases(bra,1))
+end
+
+# operator * operator -> operator
+function *{B<:OrthonormalBasis}(qm::QuMatrix{B}, dm::DualMatrix{B})
+    return QuArray(A_mul_Bc(rawcoeffs(qm), rawcoeffs(dm)), 
+                   bases(qm,1), 
+                   bases(dm,2))
+end
+
+function *{B<:OrthonormalBasis}(dm::DualMatrix{B}, qm::QuMatrix{B})
+    return QuArray(Ac_mul_B(rawcoeffs(dm), rawcoeffs(qm)), 
+                   bases(dm,1), 
+                   bases(qm,2))
+end
+
+function *{B<:OrthonormalBasis}(qm1::QuMatrix{B}, qm2::QuMatrix{B}) 
+    return QuArray(rawcoeffs(qm1)*rawcoeffs(qm2), 
+                   bases(qm1,1), 
+                   bases(qm2,2))
+end
+
+*(dm1::DualMatrix, dm2::DualMatrix) = (dm1.qarr*dm2.qarr)'

src/arrays/arrays.jl

@@ -2,7 +2,7 @@
 # Include Statements #
 ######################
     include("quarray.jl")
-    include("dual.jl")
+    include("arraymath.jl")
     include("constructors.jl")
     include("ladderops.jl")
 

src/arrays/quarray.jl

@@ -20,7 +20,7 @@
 
     QuArray{B<:AbstractBasis,T,N}(coeffs::AbstractArray{T,N}, bases::NTuple{N,B}) = QuArray{B,T,N,typeof(coeffs)}(coeffs, bases)    
     QuArray(coeffs, bases::AbstractBasis...) = QuArray(coeffs, bases)
-    QuArray(coeffs) = QuArray(coeffs, basesfordims(size(coeffs)))
+    QuArray(coeffs) = QuArray(coeffs, map(FiniteBasis, size(coeffs)))
 
     typealias QuVector{B<:AbstractBasis,T,A} QuArray{B,T,1,A}
     typealias QuMatrix{B<:AbstractBasis,T,A} QuArray{B,T,2,A}
@@ -46,24 +46,67 @@
 
     Base.in(c, qarr::QuArray) = in(c, rawcoeffs(qarr))
 
+##############
+# CTranspose #
+##############
+    immutable CTranspose{B,T,N,A} <: AbstractQuArray{B,T,N}
+        qarr::QuArray{B,T,N,A}
+        CTranspose(qarr::QuVector{B,T,A}) = new(qarr)
+        CTranspose(qarr::QuMatrix{B,T,A}) = new(qarr)
+        CTranspose(qarr::QuArray) = error("Conjugate transpose is unsupported for QuArrays of dimension $N")
+    end
+
+    CTranspose{B,T,N,A}(qa::QuArray{B,T,N,A}) = CTranspose{B,T,N,A}(qa)
+
+    typealias DualVector{B,T,A} CTranspose{B,T,1,A}
+    typealias DualMatrix{B,T,A} CTranspose{B,T,2,A}
+
     ######################
-    # Printing Functions #
+    # Accessor functions #
     ######################
-    Base.summary{B}(qarr::AbstractQuArray{B}) = "$(sizenotation(size(qarr))) $(typenotation(qarr)) in $B"              
+    rawcoeffs(ct::CTranspose) = rawcoeffs(ct.qarr)
+    coefftype{B,T,N,A}(::CTranspose{B,T,N,A}) = A
+
+    revind(len, i) = len - (i-1)
+    bases(ct::CTranspose, i) = bases(ct.qarr, revind(ndims(ct), i))
+
+    ########################
+    # Array-like functions #
+    ########################
+    Base.ndims(ct::CTranspose) = ndims(ct.qarr)
+    Base.length(ct::CTranspose) = length(ct.qarr)
+
+    Base.size(ct::CTranspose) = reverse(size(ct.qarr))
+    Base.size(ct::CTranspose, i) = size(ct, revind(ndims(ct), i))
+
+    Base.getindex(ct::CTranspose, i, j) = getindex(ct.qarr, j, i)'
+    Base.getindex(dv::DualVector, i) = getindex(dv.qarr, i)'
+
+    Base.setindex!(ct::CTranspose, x, i, j) = setindex!(ct.qarr, x', j, i)
+    Base.setindex!(dv::DualVector, x, i) = setindex!(dv.qarr, x', i)
+
+    Base.ctranspose(qarr::QuArray) = CTranspose(qarr)
+    Base.ctranspose(ct::CTranspose) = ct.qarr
+
+######################
+# Printing Functions #
+######################
+    Base.summary{B}(qarr::AbstractQuArray{B}) = "$(sizenotation(size(qarr))) $(typerepr(qarr)) in $B"              
 
     function Base.show(io::IO, qarr::AbstractQuArray)
         println(io, summary(qarr)*":")
-        println(io, "...transposed/conjugated?: $(istran(qarr))/$(isconj(qarr))")
         println(io, "...original coefficients: $(coefftype(qarr))")
         print(io, repr(rawcoeffs(qarr)))
     end
 
-    ####################
-    # Helper Functions #
-    ####################
-    typenotation(::AbstractQuVector) = "QuVector"
-    typenotation(::AbstractQuMatrix) = "QuMatrix"
-    typenotation(::AbstractQuArray) = "QuArray"
+####################
+# Helper Functions #
+####################
+    typerepr(::QuVector) = "QuVector"
+    typerepr(::QuMatrix) = "QuMatrix"
+    typerepr(::DualVector) = "DualVector"
+    typerepr(::DualMatrix) = "DualMatrix"
+    typerepr(::QuArray) = "QuArray"
 
     sizenotation(tup::(Int,)) = "$(first(tup))-element"
     sizenotation(tup::(Int...)) = reduce(*, map(s->"$(s)x", tup))[1:end-1] 
@@ -75,10 +118,4 @@
         return false
     end
 
-    # Assumes that every basis type passed in
-    # has a constructor B(::eltype(lens))
-    function basesfordims(lens::Tuple, B=ntuple(length(lens), x->FiniteBasis))
-        return ntuple(length(lens), n->B[n](lens[n]))
-    end
-
 export QuArray

src/bases/finitebasis.jl

@@ -24,7 +24,7 @@
         FiniteBasis(lens::Int...) = new(lens)
     end
 
-    FiniteBasis{N, S<:AbstractStructure}(lens::NTuple{N,Int}, ::Type{S}=AbstractStructure) = FiniteBasis{S,N}(lens)
+    FiniteBasis{N, S<:AbstractStructure}(lens::NTuple{N,Int}, ::Type{S}=Orthonormal) = FiniteBasis{S,N}(lens)
     FiniteBasis(lens::Int...) = FiniteBasis(lens)
 
     Base.convert{S,N}(::Type{FiniteBasis{S,N}}, f::FiniteBasis) = FiniteBasis(f.lens, S)

test/multest.jl

@@ -0,0 +1,23 @@
+m = [1+1im 2+2im 3+3im; 
+     4+4im 5+5im 6+6im; 
+     7+7im 8+8im 9+9im]
+
+v = [1.0,2.0,3.0]
+
+qm = QuArray(m)
+qv = QuArray(v)
+
+@assert rawcoeffs(qm*qm) == m*m
+@assert rawcoeffs(qm'*qm') == m*m
+@assert rawcoeffs(qm'*qm) == m'*m
+@assert rawcoeffs(qm*qm') == m*m'
+
+@assert rawcoeffs(qm*qv) == m*v
+@assert rawcoeffs(qv'*qm') == m*v
+@assert rawcoeffs(qm'*qv) == m'*v
+@assert rawcoeffs(qv'*qm) == m'*v
+
+@assert qv'*qv == dot(v,v)
+@assert rawcoeffs(qv*qv') == v*v'
+
+@assert qv'*qm*qv == first(v'*m*v)

test/runtests.jl

@@ -1,3 +1,5 @@
 using QuBase
+using QuBase.rawcoeffs
 using Base.Test
 
+include("multest.jl")