Skip to content

Commit

Permalink
Add Transpose immutable
Browse files Browse the repository at this point in the history
  • Loading branch information
@andreasnoack
andreasnoack committed May 14, 2014
1 parent f4d1e94 commit 37ebe9b
Show file tree
Hide file tree
Showing 5 changed files with 128 additions and 94 deletions.
52 changes: 0 additions & 52 deletions base/array.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1237,58 +1237,6 @@ function filter(f::Function, a::Vector)
return r
end

## Transpose ##

const sqrthalfcache = 1<<7
function transpose!{T<:Number}(B::Matrix{T}, A::Matrix{T})
m, n = size(A)
if size(B) != (n,m)
error("input and output must have same size")
end
elsz = isbits(T) ? sizeof(T) : sizeof(Ptr)
blocksize = ifloor(sqrthalfcache/elsz/1.4) # /1.4 to avoid complete fill of cache
if m*n <= 4*blocksize*blocksize
# For small sizes, use a simple linear-indexing algorithm
for i2 = 1:n
j = i2
offset = (j-1)*m
for i = offset+1:offset+m
B[j] = A[i]
j += n
end
end
return B
end
# For larger sizes, use a cache-friendly algorithm
for outer2 = 1:blocksize:size(A, 2)
for outer1 = 1:blocksize:size(A, 1)
for inner2 = outer2:min(n,outer2+blocksize)
i = (inner2-1)*m + outer1
j = inner2 + (outer1-1)*n
for inner1 = outer1:min(m,outer1+blocksize)
B[j] = A[i]
i += 1
j += n
end
end
end
end
B
end

function transpose{T<:Number}(A::Matrix{T})
B = similar(A, size(A, 2), size(A, 1))
transpose!(B, A)
end

ctranspose{T<:Real}(A::StridedVecOrMat{T}) = transpose(A)

transpose(x::StridedVector) = [ x[j] for i=1, j=1:size(x,1) ]
transpose(x::StridedMatrix) = [ x[j,i] for i=1:size(x,2), j=1:size(x,1) ]

ctranspose{T}(x::StridedVector{T}) = T[ conj(x[j]) for i=1, j=1:size(x,1) ]
ctranspose{T}(x::StridedMatrix{T}) = T[ conj(x[j,i]) for i=1:size(x,2), j=1:size(x,1) ]

# set-like operators for vectors
# These are moderately efficient, preserve order, and remove dupes.

Expand Down
4 changes: 3 additions & 1 deletion base/linalg.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
module LinAlg

importall Base
import Base: USE_BLAS64, size, copy, copy_transpose!, power_by_squaring, print_matrix, transpose!
import Base: USE_BLAS64, size, copy, copy_transpose!, power_by_squaring, print_matrix

export
# Modules
Expand Down Expand Up @@ -110,6 +110,7 @@ export
svdvals,
trace,
transpose,
ctranspose,
tril,
triu,
tril!,
Expand Down Expand Up @@ -189,6 +190,7 @@ Valid choices are 'U' (upper) or 'L' (lower).""")))
end

include("linalg/exceptions.jl")
include("linalg/transpose.jl")
include("linalg/generic.jl")

include("linalg/blas.jl")
Expand Down
16 changes: 10 additions & 6 deletions base/linalg/matmul.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ arithtype(T) = T
arithtype(::Type{Bool}) = Int

# multiply by diagonal matrix as vector
function scale!(C::Matrix, A::Matrix, b::Vector)
function scale!(C::AbstractMatrix, A::AbstractMatrix, b::AbstractVector)
m, n = size(A)
n==length(b) || throw(DimensionMismatch(""))
for j = 1:n
Expand All @@ -16,20 +16,21 @@ function scale!(C::Matrix, A::Matrix, b::Vector)
C
end

function scale!(C::Matrix, b::Vector, A::Matrix)
function scale!(C::AbstractMatrix, b::AbstractVector, A::AbstractMatrix)
m, n = size(A)
m==length(b) || throw(DimensionMismatch(""))
for j=1:n, i=1:m
C[i,j] = A[i,j]*b[i]
end
C
end
scale(A::Matrix, b::Vector) = scale!(similar(A, promote_type(eltype(A),eltype(b))), A, b)
scale(b::Vector, A::Matrix) = scale!(similar(b, promote_type(eltype(A),eltype(b)), size(A)), b, A)
scale(A::AbstractMatrix, b::AbstractVector) = scale!(similar(A, promote_type(eltype(A),eltype(b))), A, b)
scale(b::AbstractVector, A::AbstractMatrix) = scale!(similar(b, promote_type(eltype(A),eltype(b)), size(A)), b, A)

# Dot products

dot{T<:BlasReal}(x::Vector{T}, y::Vector{T}) = BLAS.dot(x, y)
# dot{T<:BlasReal}(x::Vector{T}, y::Vector{T}) = BLAS.dot(x, y)
*{T<:BlasReal,Conj}(x::VectorTranspose{T,Vector{T},Conj}, y::Vector{T}) = BLAS.dot(x.data, y)
dot{T<:BlasComplex}(x::Vector{T}, y::Vector{T}) = BLAS.dotc(x, y)
function dot{T<:BlasReal, TI<:Integer}(x::Vector{T}, rx::Union(UnitRange{TI},Range{TI}), y::Vector{T}, ry::Union(UnitRange{TI},Range{TI}))
length(rx)==length(ry) || throw(DimensionMismatch(""))
Expand Down Expand Up @@ -71,7 +72,10 @@ function (*){T,S}(A::AbstractMatrix{T}, x::AbstractVector{S})
TS = promote_type(arithtype(T),arithtype(S))
A_mul_B!(similar(x,TS,size(A,1)),A,x)
end
(*)(A::AbstractVector, B::AbstractMatrix) = reshape(A,length(A),1)*B

# (*)(A::AbstractVector, B::AbstractMatrix) = reshape(A,length(A),1)*B
*(x::Covector, A::AbstractMatrix) = (A'x')'
*(A::Adjoint, x::AbstractVector) = gemv!(similar(x), 'C', A.data, x)

A_mul_B!{T<:BlasFloat}(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}) = gemv!(y, 'N', A, x)
for elty in (Float32,Float64)
Expand Down
80 changes: 80 additions & 0 deletions base/linalg/transpose.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
immutable Transpose{T, N, S<:AbstractArray{T, N}, Conj} <: AbstractArray{T, N}
data::S
end
typealias MatrixTranspose{T, S, Conj} Transpose{T, 2, S, Conj}
typealias VectorTranspose{T, S, Conj} Transpose{T, 1, S, Conj}
typealias Adjoint{T, S} Transpose{T, 2, S, true}
typealias Covector{T, S} Transpose{T, 1, S, true}

ctranspose{T, N}(A::AbstractArray{T, N}) = N <= 2 ? Transpose{T, N, typeof(A), true}(A) : throw(ArgumentError("dimension cannot be larger than two"))
transpose{T<:Real, N}(A::AbstractArray{T, N}) = N <= 2 ? Transpose{T, N, typeof(A), true}(A) : throw(ArgumentError("dimension cannot be larger than two"))
transpose{T, N}(A::AbstractArray{T, N}) = N <= 2 ? Transpose{T, N, typeof(A), false}(A) : throw(ArgumentError("dimension cannot be larger than two"))

ctranspose{T,N,S}(A::Transpose{T,N,S,true}) = A.data
transpose{T,N,S}(A::Transpose{T,N,S,false}) = A.data

size(A::VectorTranspose, args...) = size(A.data, args...)
size(A::MatrixTranspose) = reverse(size(A.data))
size(A::MatrixTranspose, dim::Integer) = dim == 1 ? size(A.data, 2) : (dim == 2 ? size(A.data, 1) : size(A.data, dim))

getindex(A::VectorTranspose, i::Integer) = getindex(A.data, i)
getindex(A::MatrixTranspose, i::Integer, j::Integer) = getindex(A.data, j, i)

This comment has been minimized.

Sign in to view

Copy link
@stevengj

stevengj May 14, 2014

Member

Why no setindex!?

This comment has been minimized.

Sign in to view

Copy link
@StefanKarpinski

StefanKarpinski May 14, 2014

Member

I think that @andreasnoackjensen was worried about modifying the original array via a view, but I think that's where we're headed with making slices array views anyway, so we might as well start here.

This comment has been minimized.

Sign in to view

Copy link
@andreasnoack

andreasnoack May 14, 2014

Author Member

This is only a sketch so I have only a minimum of methods. If the idea is accepted I'll add setindex! and other relevant array methods.

## Transpose ##

const sqrthalfcache = 1<<7
function transpose!{T<:Number}(B::Matrix{T}, A::Matrix{T})

This comment has been minimized.

Sign in to view

Copy link
@stevengj

stevengj May 14, 2014

Member

Couldn't we use a cache-oblivious algorithm for this? See #6690.

This comment has been minimized.

Sign in to view

Copy link
@andreasnoack

andreasnoack May 14, 2014

Author Member

I just copied the version we already had into this file. When @Jutho's work is done it should go here I guess.
m, n = size(A)
if size(B) != (n,m)
error("input and output must have same size")
end
elsz = isbits(T) ? sizeof(T) : sizeof(Ptr)
blocksize = ifloor(sqrthalfcache/elsz/1.4) # /1.4 to avoid complete fill of cache
if m*n <= 4*blocksize*blocksize
# For small sizes, use a simple linear-indexing algorithm
for i2 = 1:n
j = i2
offset = (j-1)*m
for i = offset+1:offset+m
B[j] = A[i]
j += n
end
end
return B
end
# For larger sizes, use a cache-friendly algorithm
for outer2 = 1:blocksize:size(A, 2)
for outer1 = 1:blocksize:size(A, 1)
for inner2 = outer2:min(n,outer2+blocksize)
i = (inner2-1)*m + outer1
j = inner2 + (outer1-1)*n
for inner1 = outer1:min(m,outer1+blocksize)
B[j] = A[i]
i += 1
j += n
end
end
end
end
B
end

function full{T, S<:DenseMatrix}(A::MatrixTranspose{T, S, false})
B = similar(A, size(A, 2), size(A, 1))
transpose!(B, A)
end
function full{T, S<:DenseMatrix}(A::MatrixTranspose{T, S, true})
B = similar(A, size(A, 2), size(A, 1))
transpose!(B, A)
return conj!(B)
end

# full{T<:Real, S}(A::MatrixTranspose{T, S, true}) = transpose(A)

full(x::VectorTranspose) = x.data
full{T, S<:AbstractMatrix}(X::MatrixTranspose{T, S, false}) = [ X[i,j] for i=1:size(X,1), j=1:size(X,2) ]
full{T, S<:AbstractMatrix}(X::MatrixTranspose{T, S, true}) = [ conj(X[i,j]) for i=1:size(X,1), j=1:size(X,2) ]

# *(x::VectorTranspose, y::AbstractVector) = dot(x.data, y)

# *{T, S}(x::Covector{T, S}, A::AbstractMatrix{T}) = Transpose{T, 1, S, true}(Ac_mul_B(A, x.data))
70 changes: 35 additions & 35 deletions src/julia-syntax.scm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1442,29 +1442,29 @@
(call (top kwcall) ,f ,(length keys) ,@keyargs
,container ,@pa))))))))

(define (expand-transposed-op e ops)
(let ((a (caddr e))
(b (cadddr e)))
(cond ((ctrans? a)
(if (ctrans? b)
`(call ,(aref ops 0) #;Ac_mul_Bc ,(expand-forms (cadr a))
,(expand-forms (cadr b)))
`(call ,(aref ops 1) #;Ac_mul_B ,(expand-forms (cadr a))
,(expand-forms b))))
((trans? a)
(if (trans? b)
`(call ,(aref ops 2) #;At_mul_Bt ,(expand-forms (cadr a))
,(expand-forms (cadr b)))
`(call ,(aref ops 3) #;At_mul_B ,(expand-forms (cadr a))
,(expand-forms b))))
((ctrans? b)
`(call ,(aref ops 4) #;A_mul_Bc ,(expand-forms a)
,(expand-forms (cadr b))))
((trans? b)
`(call ,(aref ops 5) #;A_mul_Bt ,(expand-forms a)
,(expand-forms (cadr b))))
(else
`(call ,(cadr e) ,(expand-forms a) ,(expand-forms b))))))
; (define (expand-transposed-op e ops)
; (let ((a (caddr e))
; (b (cadddr e)))
; (cond ((ctrans? a)
; (if (ctrans? b)
; `(call ,(aref ops 0) #;Ac_mul_Bc ,(expand-forms (cadr a))
; ,(expand-forms (cadr b)))
; `(call ,(aref ops 1) #;Ac_mul_B ,(expand-forms (cadr a))
; ,(expand-forms b))))
; ((trans? a)
; (if (trans? b)
; `(call ,(aref ops 2) #;At_mul_Bt ,(expand-forms (cadr a))
; ,(expand-forms (cadr b)))
; `(call ,(aref ops 3) #;At_mul_B ,(expand-forms (cadr a))
; ,(expand-forms b))))
; ((ctrans? b)
; `(call ,(aref ops 4) #;A_mul_Bc ,(expand-forms a)
; ,(expand-forms (cadr b))))
; ((trans? b)
; `(call ,(aref ops 5) #;A_mul_Bt ,(expand-forms a)
; ,(expand-forms (cadr b))))
; (else
; `(call ,(cadr e) ,(expand-forms a) ,(expand-forms b))))))

(define (lower-update-op e)
(expand-forms
Expand Down Expand Up @@ -1665,18 +1665,18 @@
(expand-forms
`(call (top apply) ,f ,@(tuple-wrap argl '())))))

((and (eq? (cadr e) '*) (length= e 4))
(expand-transposed-op
e
#(Ac_mul_Bc Ac_mul_B At_mul_Bt At_mul_B A_mul_Bc A_mul_Bt)))
((and (eq? (cadr e) '/) (length= e 4))
(expand-transposed-op
e
#(Ac_rdiv_Bc Ac_rdiv_B At_rdiv_Bt At_rdiv_B A_rdiv_Bc A_rdiv_Bt)))
((and (eq? (cadr e) '\\) (length= e 4))
(expand-transposed-op
e
#(Ac_ldiv_Bc Ac_ldiv_B At_ldiv_Bt At_ldiv_B A_ldiv_Bc A_ldiv_Bt)))
; ((and (eq? (cadr e) '*) (length= e 4))
; (expand-transposed-op
; e
; #(Ac_mul_Bc Ac_mul_B At_mul_Bt At_mul_B A_mul_Bc A_mul_Bt)))
; ((and (eq? (cadr e) '/) (length= e 4))
; (expand-transposed-op
; e
; #(Ac_rdiv_Bc Ac_rdiv_B At_rdiv_Bt At_rdiv_B A_rdiv_Bc A_rdiv_Bt)))
; ((and (eq? (cadr e) '\\) (length= e 4))
; (expand-transposed-op
; e
; #(Ac_ldiv_Bc Ac_ldiv_B At_ldiv_Bt At_ldiv_B A_ldiv_Bc A_ldiv_Bt)))
(else
(map expand-forms e))))
(map expand-forms e)))
Expand Down

0 comments on commit 37ebe9b

Please sign in to comment.