Skip to content

Commit

Permalink
Merge pull request #7244 from JuliaLang/anj/block
Browse files Browse the repository at this point in the history 
Make (c)transpose work correctly for block matrices
  • Loading branch information
@andreasnoack
andreasnoack committed Jun 17, 2014
2 parents 80abdd3 + ab795c8 commit 4cfb9a1
Show file tree
Hide file tree
Showing 4 changed files with 48 additions and 43 deletions.
12 changes: 6 additions & 6 deletions base/array.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1235,7 +1235,7 @@ function transpose!(B::StridedMatrix,A::StridedMatrix)
@inbounds begin
for j = 1:n
for i = 1:m
B[j,i] = A[i,j]
B[j,i] = transpose(A[i,j])
end
end
end
Expand All @@ -1249,7 +1249,7 @@ function transposeblock!(B::StridedMatrix,A::StridedMatrix,m::Int,n::Int,offseti
@inbounds begin
for j = offsetj+(1:n)
for i = offseti+(1:m)
B[j,i] = A[i,j]
B[j,i] = transpose(A[i,j])
end
end
end
Expand All @@ -1272,7 +1272,7 @@ function ctranspose!(B::StridedMatrix,A::StridedMatrix)
@inbounds begin
for j = 1:n
for i = 1:m
B[j,i] = conj(A[i,j])
B[j,i] = ctranspose(A[i,j])
end
end
end
Expand All @@ -1286,7 +1286,7 @@ function ctransposeblock!(B::StridedMatrix,A::StridedMatrix,m::Int,n::Int,offset
@inbounds begin
for j = offsetj+(1:n)
for i = offseti+(1:m)
B[j,i] = conj(A[i,j])
B[j,i] = ctranspose(A[i,j])
end
end
end
Expand All @@ -1312,8 +1312,8 @@ function ctranspose(A::StridedMatrix)
end
ctranspose{T<:Real}(A::StridedVecOrMat{T}) = transpose(A)

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

# set-like operators for vectors
# These are moderately efficient, preserve order, and remove dupes.
Expand Down
4 changes: 2 additions & 2 deletions base/linalg/generic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,7 +245,7 @@ function issym(A::AbstractMatrix)
m, n = size(A)
m==n || return false
for i = 1:(n-1), j = (i+1):n
if A[i,j] != A[j,i]
if A[i,j] != transpose(A[j,i])
return false
end
end
Expand All @@ -258,7 +258,7 @@ function ishermitian(A::AbstractMatrix)
m, n = size(A)
m==n || return false
for i = 1:n, j = i:n
if A[i,j] != conj(A[j,i])
if A[i,j] != ctranspose(A[j,i])
return false
end
end
Expand Down
70 changes: 35 additions & 35 deletions base/linalg/matmul.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -325,7 +325,7 @@ function generic_matvecmul!{T,S,R}(C::AbstractVector{R}, tA, A::AbstractMatrix{T
aoffs = (k-1)*Astride
s = z
for i = 1:nA
s += A[aoffs+i] * B[i]
s += A[aoffs+i].'B[i]
end
C[k] = s
end
Expand All @@ -334,11 +334,11 @@ function generic_matvecmul!{T,S,R}(C::AbstractVector{R}, tA, A::AbstractMatrix{T
aoffs = (k-1)*Astride
s = z
for i = 1:nA
s += conj(A[aoffs+i]) * B[i]
s += A[aoffs+i]'B[i]
end
C[k] = s
end
else # tA == 'N'
else # tA == 'N'
fill!(C, z)
for k = 1:mB
aoffs = (k-1)*Astride
Expand Down Expand Up @@ -438,69 +438,69 @@ function generic_matmatmul!{T,S,R}(C::AbstractVecOrMat{R}, tA, tB, A::AbstractVe
end
elseif tB == 'T'
for i = 1:mA, j = 1:nB
Ctmp = A[i, 1]*B[j, 1]
Ctmp = A[i, 1]*B[j, 1].'
for k = 2:nA
Ctmp += A[i, k]*B[j, k]
Ctmp += A[i, k]*B[j, k].'
end
C[i,j] = Ctmp
end
else
for i = 1:mA, j = 1:nB
Ctmp = A[i, 1]*conj(B[j, 1])
Ctmp = A[i, 1]*B[j, 1]'
for k = 2:nA
Ctmp += A[i, k]*conj(B[j, k])
Ctmp += A[i, k]*B[j, k]'
end
C[i,j] = Ctmp
end
end
elseif tA == 'T'
if tB == 'N'
for i = 1:mA, j = 1:nB
Ctmp = A[1, i]*B[1, j]
Ctmp = A[1, i].'B[1, j]
for k = 2:nA
Ctmp += A[k, i]*B[k, j]
Ctmp += A[k, i].'B[k, j]
end
C[i,j] = Ctmp
end
elseif tB == 'T'
for i = 1:mA, j = 1:nB
Ctmp = A[1, i]*B[j, 1]
Ctmp = A[1, i].'B[j, 1].'
for k = 2:nA
Ctmp += A[k, i]*B[j, k]
Ctmp += A[k, i].'B[j, k].'
end
C[i,j] = Ctmp
end
else
for i = 1:mA, j = 1:nB
Ctmp = A[1, i]*conj(B[j, 1])
Ctmp = A[1, i].'B[j, 1]'
for k = 2:nA
Ctmp += A[k, i]*conj(B[j, k])
Ctmp += A[k, i].'B[j, k]'
end
C[i,j] = Ctmp
end
end
else
if tB == 'N'
for i = 1:mA, j = 1:nB
Ctmp = conj(A[1, i])*B[1, j]
Ctmp = A[1, i]'B[1, j]
for k = 2:nA
Ctmp += conj(A[k, i])*B[k, j]
Ctmp += A[k, i]'B[k, j]
end
C[i,j] = Ctmp
end
elseif tB == 'T'
for i = 1:mA, j = 1:nB
Ctmp = conj(A[1, i])*B[j, 1]
Ctmp = A[1, i]'B[j, 1].'
for k = 2:nA
Ctmp += conj(A[k, i])*B[j, k]
Ctmp += A[k, i]'B[j, k].'
end
C[i,j] = Ctmp
end
else
for i = 1:mA, j = 1:nB
Ctmp = conj(A[1, i])*conj(B[j, 1])
Ctmp = A[1, i]'B[j, 1]'
for k = 2:nA
Ctmp += conj(A[k, i])*conj(B[j, k])
Ctmp += A[k, i]'B[j, k]'
end
C[i,j] = Ctmp
end
Expand All @@ -519,16 +519,16 @@ end

function matmul2x2!{T,S,R}(C::AbstractMatrix{R}, tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
if tA == 'T'
A11 = A[1,1]; A12 = A[2,1]; A21 = A[1,2]; A22 = A[2,2]
A11 = transpose(A[1,1]); A12 = transpose(A[2,1]); A21 = transpose(A[1,2]); A22 = transpose(A[2,2])
elseif tA == 'C'
A11 = conj(A[1,1]); A12 = conj(A[2,1]); A21 = conj(A[1,2]); A22 = conj(A[2,2])
A11 = ctranspose(A[1,1]); A12 = ctranspose(A[2,1]); A21 = ctranspose(A[1,2]); A22 = ctranspose(A[2,2])
else
A11 = A[1,1]; A12 = A[1,2]; A21 = A[2,1]; A22 = A[2,2]
end
if tB == 'T'
B11 = B[1,1]; B12 = B[2,1]; B21 = B[1,2]; B22 = B[2,2]
B11 = transpose(B[1,1]); B12 = transpose(B[2,1]); B21 = transpose(B[1,2]); B22 = transpose(B[2,2])
elseif tB == 'C'
B11 = conj(B[1,1]); B12 = conj(B[2,1]); B21 = conj(B[1,2]); B22 = conj(B[2,2])
B11 = ctranspose(B[1,1]); B12 = ctranspose(B[2,1]); B21 = ctranspose(B[1,2]); B22 = ctranspose(B[2,2])
else
B11 = B[1,1]; B12 = B[1,2]; B21 = B[2,1]; B22 = B[2,2]
end
Expand All @@ -546,27 +546,27 @@ end

function matmul3x3!{T,S,R}(C::AbstractMatrix{R}, tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
if tA == 'T'
A11 = A[1,1]; A12 = A[2,1]; A13 = A[3,1];
A21 = A[1,2]; A22 = A[2,2]; A23 = A[3,2];
A31 = A[1,3]; A32 = A[2,3]; A33 = A[3,3];
A11 = transpose(A[1,1]); A12 = transpose(A[2,1]); A13 = transpose(A[3,1]);
A21 = transpose(A[1,2]); A22 = transpose(A[2,2]); A23 = transpose(A[3,2]);
A31 = transpose(A[1,3]); A32 = transpose(A[2,3]); A33 = transpose(A[3,3]);
elseif tA == 'C'
A11 = conj(A[1,1]); A12 = conj(A[2,1]); A13 = conj(A[3,1]);
A21 = conj(A[1,2]); A22 = conj(A[2,2]); A23 = conj(A[3,2]);
A31 = conj(A[1,3]); A32 = conj(A[2,3]); A33 = conj(A[3,3]);
A11 = ctranspose(A[1,1]); A12 = ctranspose(A[2,1]); A13 = ctranspose(A[3,1]);
A21 = ctranspose(A[1,2]); A22 = ctranspose(A[2,2]); A23 = ctranspose(A[3,2]);
A31 = ctranspose(A[1,3]); A32 = ctranspose(A[2,3]); A33 = ctranspose(A[3,3]);
else
A11 = A[1,1]; A12 = A[1,2]; A13 = A[1,3];
A21 = A[2,1]; A22 = A[2,2]; A23 = A[2,3];
A31 = A[3,1]; A32 = A[3,2]; A33 = A[3,3];
end

if tB == 'T'
B11 = B[1,1]; B12 = B[2,1]; B13 = B[3,1];
B21 = B[1,2]; B22 = B[2,2]; B23 = B[3,2];
B31 = B[1,3]; B32 = B[2,3]; B33 = B[3,3];
B11 = transpose(B[1,1]); B12 = transpose(B[2,1]); B13 = transpose(B[3,1]);
B21 = transpose(B[1,2]); B22 = transpose(B[2,2]); B23 = transpose(B[3,2]);
B31 = transpose(B[1,3]); B32 = transpose(B[2,3]); B33 = transpose(B[3,3]);
elseif tB == 'C'
B11 = conj(B[1,1]); B12 = conj(B[2,1]); B13 = conj(B[3,1]);
B21 = conj(B[1,2]); B22 = conj(B[2,2]); B23 = conj(B[3,2]);
B31 = conj(B[1,3]); B32 = conj(B[2,3]); B33 = conj(B[3,3]);
B11 = ctranspose(B[1,1]); B12 = ctranspose(B[2,1]); B13 = ctranspose(B[3,1]);
B21 = ctranspose(B[1,2]); B22 = ctranspose(B[2,2]); B23 = ctranspose(B[3,2]);
B31 = ctranspose(B[1,3]); B32 = ctranspose(B[2,3]); B33 = ctranspose(B[3,3]);
else
B11 = B[1,1]; B12 = B[1,2]; B13 = B[1,3];
B21 = B[2,1]; B22 = B[2,2]; B23 = B[2,3];
Expand Down
5 changes: 5 additions & 0 deletions base/operators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,11 @@ mod1{T<:Real}(x::T, y::T) = y-mod(y-x,y)
rem1{T<:Real}(x::T, y::T) = rem(x-1,y)+1
fld1{T<:Real}(x::T, y::T) = fld(x-1,y)+1

# transpose
transpose(x) = x
ctranspose(x) = conj(transpose(x))
conj(x) = x

# transposed multiply
Ac_mul_B (a,b) = ctranspose(a)*b
A_mul_Bc (a,b) = a*ctranspose(b)
Expand Down

0 comments on commit 4cfb9a1

Please sign in to comment.