Skip to content

Commit

Permalink
Support for complex matrices
Browse files Browse the repository at this point in the history 
Enable travis testing with BigFloat
  • Loading branch information
@andreasnoack
andreasnoack committed Jul 8, 2016
1 parent 8c6acef commit 758640a
Show file tree
Hide file tree
Showing 4 changed files with 204 additions and 129 deletions.
149 changes: 95 additions & 54 deletions src/ToeplitzMatrices.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,10 +54,10 @@ function A_mul_B!{T}(α::T, A::AbstractToeplitz{T}, x::StridedVector{T}, β::T,
return y
end
for i = 1:n
A.tmp[i] = complex(x[i])
A.tmp[i] = x[i]
end
for i = n+1:N
A.tmp[i] = zero(Complex{T})
A.tmp[i] = 0
end
A_mul_B!(A.tmp, A.dft, A.tmp)
for i = 1:N
Expand All @@ -66,7 +66,7 @@ function A_mul_B!{T}(α::T, A::AbstractToeplitz{T}, x::StridedVector{T}, β::T,
A.dft \ A.tmp
for i = 1:m
y[i] *= β
y[i] += α * real(A.tmp[i])
y[i] += α * (T <: Real ? real(A.tmp[i]) : A.tmp[i])
end
return y
end
Expand All @@ -84,9 +84,15 @@ function A_mul_B!{T}(α::T, A::AbstractToeplitz{T}, B::StridedMatrix{T}, β::T,
return C
end

# Translate three to five argument A_mul_B!
A_mul_B!(y, A::AbstractToeplitz, x) = A_mul_B!(one(eltype(A)), A, x, zero(eltype(A)), x)

# * operator
(*){T}(A::AbstractToeplitz{T}, B::StridedVecOrMat{T}) =
A_mul_B!(one(T), A, B, zero(T), size(B,2) == 1 ? zeros(T, size(A, 1)) : zeros(T, size(A, 1), size(B, 2)))
function (*)(A::AbstractToeplitz, B::StridedVecOrMat)
T = promote_type(eltype(A), eltype(B))
A_mul_B!(one(T), convert(AbstractMatrix{T}, A), convert(AbstractArray{T}, B), zero(T),
size(B,2) == 1 ? zeros(T, size(A, 1)) : zeros(T, size(A, 1), size(B, 2)))
end

# Left division of a general matrix B by a general Toeplitz matrix A, i.e. the solution x of Ax=B.
function A_ldiv_B!(A::AbstractToeplitz, B::StridedMatrix)
Expand All @@ -99,37 +105,42 @@ function A_ldiv_B!(A::AbstractToeplitz, B::StridedMatrix)
return B
end

function (\)(A::AbstractToeplitz, b::AbstractVector)
T = promote_type(eltype(A), eltype(b))
if T != eltype(A)
throw(ArgumentError("promotion of Toeplitz matrices not handled yet"))
end
bb = similar(b, T)
copy!(bb, b)
A_ldiv_B!(A, bb)
end

# General Toeplitz matrix
type Toeplitz{T<:Number} <: AbstractToeplitz{T}
type Toeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
vc::Vector{T}
vr::Vector{T}
vcvr_dft::Vector{Complex{T}}
tmp::Vector{Complex{T}}
dft::Base.DFT.Plan{Complex{T}}
vcvr_dft::Vector{S}
tmp::Vector{S}
dft::Base.DFT.Plan{S}
end

# Ctor
function Toeplitz{T<:Number}(vc::Vector{T}, vr::Vector{T})
m = length(vc)
function Toeplitz(vc::Vector, vr::Vector)
m, n = length(vc), length(vr)
if vc[1] != vr[1]
error("First element of the vectors must be the same")
end

tmp = complex([vc; reverse(vr[2:end])])
dft = plan_fft!(tmp)
return Toeplitz(vc, vr, dft*tmp, similar(tmp), dft)
end

# Conversion to Float for integer inputs
Toeplitz{T<:Integer}(vc::Vector{T}, vr::Vector{T}) =
Toeplitz(Vector{Float64}(vc),Vector{Float64}(vr))
Toeplitz{T<:Integer}(vc::Vector{Complex{T}}, vr::Vector{Complex{T}}) =
Toeplitz(Vector{Complex128}(vc),Vector{Complex128}(vr))
T = promote_type(eltype(vc), eltype(vr), Float32)
vcp, vrp = Vector{T}(vc), Vector{T}(vr)

# Input promotion
function Toeplitz{T1<:Number,T2<:Number}(vc::Vector{T1},vr::Vector{T2})
T=promote_type(T1,T2)
Toeplitz(Vector{T}(vc),Vector{T}(vr))
tmp = Array(promote_type(T, Complex{Float32}), m + n - 1)
copy!(tmp, vcp)
for i = 1:n - 1
tmp[i + m] = vrp[n - i + 1]
end
dft = plan_fft!(tmp)
return Toeplitz(vcp, vrp, dft*tmp, similar(tmp), dft)
end

# Size of a general Toeplitz matrix
Expand Down Expand Up @@ -161,35 +172,35 @@ function getindex(A::Toeplitz, i::Integer, j::Integer)
end

# Form a lower triangular Toeplitz matrix by annihilating all entries above the k-th diaganal
function tril{T}(A::Toeplitz{T}, k = 0)
function tril(A::Toeplitz, k = 0)
if k > 0
error("Second argument cannot be positive")
end
Al = TriangularToeplitz(copy(A.vc), 'L', length(A.vr))
if k < 0
for i in -1:-1:k
Al.ve[-i] = zero(T)
Al.ve[-i] = zero(eltype(A))
end
end
return Al
end

# Form a lower triangular Toeplitz matrix by annihilating all entries below the k-th diaganal
function triu{T}(A::Toeplitz{T}, k = 0)
function triu(A::Toeplitz, k = 0)
if k < 0
error("Second argument cannot be negative")
end
Al = TriangularToeplitz(copy(A.vr), 'U', length(A.vc))
if k > 0
for i in 1:k
Al.ve[i] = zero(T)
Al.ve[i] = zero(eltype(A))
end
end
return Al
end

A_ldiv_B!(A::Toeplitz, b::StridedVector) =
copy!(b, IterativeLinearSolvers.cgs(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1])
copy!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, strang(A), 1000, 100eps())[1])

# Symmetric
type SymmetricToeplitz{T<:BlasReal} <: AbstractToeplitz{T}
Expand Down Expand Up @@ -218,15 +229,21 @@ A_ldiv_B!(A::SymmetricToeplitz, b::StridedVector) =
copy!(b, IterativeLinearSolvers.cg(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1])

# Circulant
type Circulant{T<:BlasReal} <: AbstractToeplitz{T}
type Circulant{T<:Number,S<:Number} <: AbstractToeplitz{T}
vc::Vector{T}
vcvr_dft::Vector{Complex{T}}
tmp::Vector{Complex{T}}
vcvr_dft::Vector{S}
tmp::Vector{S}
dft::Base.DFT.Plan
end

function Circulant{T<:BlasReal}(vc::Vector{T})
tmp = zeros(Complex{T}, length(vc))
function Circulant{T<:Real}(vc::Vector{T})
TT = promote_type(eltype(vc), Float32)
tmp = zeros(promote_type(TT, Complex{Float32}), length(vc))
return Circulant(real(vc), fft(vc), tmp, plan_fft!(tmp))
end
function Circulant(vc::Vector)
T = promote_type(eltype(vc), Float32)
tmp = zeros(promote_type(T, Complex{Float32}), length(vc))
return Circulant(vc, fft(vc), tmp, plan_fft!(tmp))
end

Expand All @@ -246,12 +263,21 @@ function getindex(C::Circulant, i::Integer, j::Integer)
return C.vc[mod(i - j, length(C.vc)) + 1]
end

function Ac_mul_B(A::Circulant,B::Circulant)
function Ac_mul_B{T<:Real}(A::Circulant{T}, B::Circulant{T})
tmp = similar(A.vcvr_dft)
for i = 1:length(tmp)
tmp[i] = conj(A.vcvr_dft[i]) * B.vcvr_dft[i]
end
return Circulant(real(A.dft \ tmp), tmp, A.tmp, A.dft)
return Circulant(real(A.vc)(A.dft \ tmp), tmp, A.tmp, A.dft)
end
function Ac_mul_B(A::Circulant, B::Circulant)
T = promote_type(eltype(A), eltype(B))
tmp = similar(A.vcvr_dft, T)
for i = 1:length(tmp)
tmp[i] = conj(A.vcvr_dft[i]) * B.vcvr_dft[i]
end
tmp2 = A.dft \ tmp
return Circulant(Vector{T}(tmp2), tmp, Vector{T}(A.tmp), eltype(A) == T ? A.dft : plan_fft!(tmp2))
end

function A_ldiv_B!{T}(C::Circulant{T}, b::AbstractVector{T})
Expand All @@ -266,16 +292,19 @@ function A_ldiv_B!{T}(C::Circulant{T}, b::AbstractVector{T})
end
C.dft \ C.tmp
for i = 1:n
b[i] = real(C.tmp[i])
b[i] = (T <: Real ? real(C.tmp[i]) : C.tmp[i])
end
return b
end
(\)(C::Circulant, b::AbstractVector) = A_ldiv_B!(C, copy(b))

function inv{T<:BlasReal}(C::Circulant{T})
function inv{T<:Real}(C::Circulant{T})
vdft = 1 ./ C.vcvr_dft
return Circulant(real(C.dft \ vdft), copy(vdft), similar(vdft), C.dft)
end
function inv(C::Circulant)
vdft = 1 ./ C.vcvr_dft
return Circulant(C.dft \ vdft, copy(vdft), similar(vdft), C.dft)
end

function strang{T}(A::AbstractMatrix{T})
n = size(A, 1)
Expand All @@ -300,28 +329,40 @@ function chan{T}(A::AbstractMatrix{T})
end

# Triangular
type TriangularToeplitz{T<:Number} <: AbstractToeplitz{T}
type TriangularToeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
ve::Vector{T}
uplo::Char
vcvr_dft::Vector{Complex{T}}
tmp::Vector{Complex{T}}
vcvr_dft::Vector{S}
tmp::Vector{S}
dft::Base.DFT.Plan
end

function Toeplitz(A::TriangularToeplitz)
if A.uplo == 'L'
Toeplitz(A.ve,[A.ve[1];zeros(length(A.ve)-1)])
function TriangularToeplitz(ve::Vector, uplo::Symbol)
n = length(ve)

T = promote_type(eltype(ve), Float32)
vep = Vector{T}(ve)

tmp = zeros(promote_type(T, Complex{Float32}), 2n - 1)
if uplo == :L
copy!(tmp, vep)
else
@assert A.uplo == 'U'
Toeplitz([A.ve[1];zeros(length(A.ve)-1)],A.ve)
tmp[1] = vep[1]
for i = 1:n - 1
tmp[n + i] = ve[n - i + 1]
end
end
dft = plan_fft!(tmp)
return TriangularToeplitz(vep, string(uplo)[1], dft * tmp, similar(tmp), dft)
end

function TriangularToeplitz{T<:Number}(ve::Vector{T}, uplo::Symbol)
n = length(ve)
tmp = uplo == :L ? complex([ve; zeros(n)]) : complex([ve[1]; zeros(T, n); reverse(ve[2:end])])
dft = plan_fft!(tmp)
return TriangularToeplitz(ve, string(uplo)[1], dft * tmp, similar(tmp), dft)
function convert(::Type{Toeplitz}, A::TriangularToeplitz)
if A.uplo == 'L'
Toeplitz(A.ve, [A.ve[1]; zeros(length(A.ve) - 1)])
else
@assert A.uplo == 'U'
Toeplitz([A.ve[1]; zeros(length(A.ve) - 1)], A.ve)
end
end

function size(A::TriangularToeplitz, dim::Int)
Expand Down Expand Up @@ -389,7 +430,7 @@ end

# A_ldiv_B!(A::TriangularToeplitz,b::StridedVector) = inv(A)*b
A_ldiv_B!(A::TriangularToeplitz, b::StridedVector) =
copy!(b, IterativeLinearSolvers.cgs(A, zeros(length(b)), b, chan(A), 1000, 100eps())[1])
copy!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, chan(A), 1000, 100eps())[1])

# extend levinson
StatsBase.levinson!(x::StridedVector, A::SymmetricToeplitz, b::StridedVector) =
Expand Down
42 changes: 25 additions & 17 deletions src/iterativeLinearSolvers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ function cg{T<:LinAlg.BlasReal}(A::AbstractMatrix{T},
return x, error, iter, flag
end

function cgs{T<:LinAlg.BlasReal}(A::AbstractMatrix{T},
function cgs{T<:Number}(A::AbstractMatrix{T},
x::AbstractVector{T},
b::AbstractVector{T},
M::Preconditioner{T},
Expand Down Expand Up @@ -130,31 +130,37 @@ function cgs{T<:LinAlg.BlasReal}(A::AbstractMatrix{T},

n = length(b)
bnrm2 = norm(b)
if bnrm2 == 0.0 bnrm2 = one(T) end
if bnrm2 == 0
bnrm2 = one(bnrm2)
end

u = zeros(T, n)
p = zeros(T, n)
p̂ = zeros(T, n)
q = zeros(T, n)
û = zeros(T,n)
û = zeros(T, n)
v̂ = zeros(T, n)
ρ = zero(T)
ρ₁ = ρ
r = copy(b)
A_mul_B!(-one(T),A,x,one(T),r)
A_mul_B!(-one(T), A, x, one(T), r)
# r = b - A*x
error = norm(r)/bnrm2

if error < tol return x, error, iter, flag end
if error < tol
return x, error, iter, flag
end

r_tld = copy(r)

for iter = 1:max_it # begin iteration
for iter = 1:max_it # begin iteration

ρ = dot(r_tld,r)
if ρ == 0.0 break end
ρ = dot(r_tld, r)
if ρ == 0
break
end

if iter > 1 # direction vectors
if iter > 1 # direction vectors
β = ρ/ρ₁
for l = 1:n
u[l] = r[l] + β*q[l]
Expand All @@ -168,9 +174,9 @@ function cgs{T<:LinAlg.BlasReal}(A::AbstractMatrix{T},
p̂[:] = p
A_ldiv_B!(M, p̂)
# p̂[:] = M\p
A_mul_B!(one(T),A,p̂,zero(T),v̂) # adjusting scalars
A_mul_B!(one(T), A, p̂, zero(T), v̂) # adjusting scalars
# v̂[:] = A*p̂
α = ρ/dot(r_tld,v̂)
α = ρ/dot(r_tld, v̂)
for l = 1:n
q[l] = u[l] - α*v̂[l]
û[l] = u[l] + q[l]
Expand All @@ -179,21 +185,23 @@ function cgs{T<:LinAlg.BlasReal}(A::AbstractMatrix{T},
# û[:] = M\û

for l = 1:n
x[l] += α*û[l] # update approximation
x[l] += α*û[l] # update approximation
end

A_mul_B!(-α,A,û,one(T),r)
A_mul_B!(-α, A, û, one(T), r)
# r[:] -= α*(A*û)
error = norm(r)/bnrm2 # check convergence
if error <= tol break end
error = norm(r)/bnrm2 # check convergence
if error <= tol
break
end

ρ₁ = ρ

end

if error <= tol # converged
if error <= tol # converged
flag = 0
elseif ρ == 0.0 # breakdown
elseif ρ == 0 # breakdown
flag = -1
else # no convergence
flag = 1
Expand Down
1 change: 1 addition & 0 deletions test/REQUIRE
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
FastTransforms
Loading

0 comments on commit 758640a

Please sign in to comment.