diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index bc3b3ea..76bdcbd 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -14,7 +14,7 @@ jobs: fail-fast: false matrix: version: - - '1.0' + - '1.6' - '1' - 'nightly' os: diff --git a/Project.toml b/Project.toml index 6ceb1bc..7dbad62 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,15 @@ name = "ToeplitzMatrices" uuid = "c751599d-da0a-543b-9d20-d0a503d91d24" -version = "0.7.0" +version = "0.8.0" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" +DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" +IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] AbstractFFTs = "0.4, 0.5, 1" StatsBase = "0.32, 0.33" -julia = "1" +julia = "1.6" diff --git a/src/ToeplitzMatrices.jl b/src/ToeplitzMatrices.jl index 704e910..f83de90 100644 --- a/src/ToeplitzMatrices.jl +++ b/src/ToeplitzMatrices.jl @@ -1,28 +1,86 @@ module ToeplitzMatrices using StatsBase - import Base: convert, *, \, getindex, print_matrix, size, Matrix, +, -, copy, similar, sqrt, copyto!, adjoint, transpose +using Base.Threads: @threads + import LinearAlgebra: Cholesky, DimensionMismatch, cholesky, cholesky!, eigvals, inv, ldiv!, mul!, pinv, rmul!, tril, triu -using LinearAlgebra: LinearAlgebra, Adjoint, Factorization, factorize +using LinearAlgebra +using LinearAlgebra: LinearAlgebra, Adjoint, Factorization, factorize, checksquare using AbstractFFTs using AbstractFFTs: Plan +using DSP: conv + flipdim(A, d) = reverse(A, dims=d) export Toeplitz, SymmetricToeplitz, Circulant, TriangularToeplitz, Hankel, chan, strang +using IterativeSolvers -include("iterativeLinearSolvers.jl") +const DEFAULT_MAXITER = 1000 # default maximum iteration count for iterative solvers +const MIN_FFT_LENGTH = 512 # minimum vector length below which the dense, non-fft based multiplication is used +const DEFAULT_TOL = 1e-6 # Abstract abstract type AbstractToeplitz{T<:Number} <: AbstractMatrix{T} end +function ldiv!(x::AbstractVector, A::AbstractToeplitz, b::AbstractVector; + verbose::Bool = false, + maxiter = DEFAULT_MAXITER, + abstol::Real = max(DEFAULT_TOL, zero(real(eltype(b)))), + reltol::Real = max(DEFAULT_TOL, sqrt(eps(real(eltype(b)))))) + n, m = size(A) + F = n + m - 1 < MIN_FFT_LENGTH ? A : factorize(A) # if we are going to use FFT-based multiplication, factorize once before iterative method + if n == m + P = factorize(strang(A)) + P = sqrt(abs(P)) + IterativeSolvers.gmres!(x, F, b, Pl = P, Pr = P', maxiter = maxiter, + abstol = abstol, reltol = reltol, + log = false, verbose = verbose) + else # if the system is rectangular, use lsqr + IterativeSolvers.lsqr!(x, F, b, maxiter = maxiter, + atol = reltol, btol = reltol, + log = false, verbose = verbose) # this has a default tolerance of 1e-6 + end +end + +function ldiv!(X::AbstractMatrix, A::AbstractToeplitz, B::AbstractMatrix) + for j in 1:size(B, 2) + ldiv!(view(X, :, j), A, view(B, :, j)) + end + return X +end +function (\)(A::AbstractToeplitz, b::AbstractVector) + T = promote_type(eltype(A), eltype(b)) + x = zeros(T, size(A, 2)) + ldiv!(x, A, b) +end +function (\)(A::AbstractToeplitz, B::AbstractMatrix) + T = promote_type(eltype(A), eltype(B)) + X = zeros(T, size(A, 2), size(B, 2)) + ldiv!(X, A, B) +end +# unclear if having these methods is of much use, since they copy +function ldiv!(A::AbstractToeplitz, b::AbstractVector) + eltype(A) == eltype(b) || throw(TypeError("storing result in b only allowed if eltype(A) == eltype(b)")) + T = promote_type(eltype(A), eltype(b)) + n = checksquare(A) + x = zeros(T, n) + copyto!(b, ldiv!(x, A, b)) +end +function ldiv!(A::AbstractToeplitz, B::AbstractMatrix) + for j in 1:size(B, 2) + ldiv!(A, view(B, :, j)) + end + return B +end + """ ToeplitzFactorization @@ -32,18 +90,50 @@ struct ToeplitzFactorization{T<:Number,A<:AbstractToeplitz{T},S<:Number,P<:Plan{ vcvr_dft::Vector{S} tmp::Vector{S} dft::P + n::Int + m::Int +end + +# doing this non-lazily simplifies implementation of mul!, ldiv! for adjoints +# of Toeplitz factorizations significantly Base.adjoint(A::ToeplitzFactorization) = Adjoint(A) +Base.adjoint(T::ToeplitzFactorization) = adjoint!(copy(T)) +function Base.copy(T::ToeplitzFactorization) + vcvr_dft = copy(T.vcvr_dft) + dft = plan_fft!(vcvr_dft) + typeof(T)(vcvr_dft, copy(T.tmp), dft, T.n, T.m) +end +# calculates the adjoint but reuses the temporary memory of T +function adjoint!(T::ToeplitzFactorization) + @. T.vcvr_dft = conj(T.vcvr_dft) + typeof(T)(T.vcvr_dft, T.tmp, T.dft, T.m, T.n) # switching n and m +end +# Base.copy() +Base.size(A::AbstractToeplitz) = (size(A, 1), size(A, 2)) +Base.size(A::ToeplitzFactorization) = (A.n, A.m) + +Base.length(A::ToeplitzFactorization) = A.n * A.m +function Base.size(A::ToeplitzFactorization, i::Int) + if i == 1 + A.n + elseif i == 2 + A.m + elseif i > 2 + 1 + else + throw(DomainError("dimension i cannot be non-positive")) + end end -size(A::AbstractToeplitz) = (size(A, 1), size(A, 2)) -function getindex(A::AbstractToeplitz, i::Integer) + +function Base.getindex(A::AbstractToeplitz, i::Integer) return A[mod(i - 1, size(A, 1)) + 1, div(i - 1, size(A, 1)) + 1] end -convert(::Type{AbstractMatrix{T}}, S::AbstractToeplitz) where {T} = convert(AbstractToeplitz{T}, S) -convert(::Type{AbstractArray{T}}, S::AbstractToeplitz) where {T} = convert(AbstractToeplitz{T}, S) +Base.convert(::Type{AbstractMatrix{T}}, S::AbstractToeplitz) where {T} = convert(AbstractToeplitz{T}, S) +Base.convert(::Type{AbstractArray{T}}, S::AbstractToeplitz) where {T} = convert(AbstractToeplitz{T}, S) # Convert an abstract Toeplitz matrix to a full matrix -function Matrix(A::AbstractToeplitz{T}) where T +function Base.Matrix(A::AbstractToeplitz{T}) where T m, n = size(A) Af = Matrix{T}(undef, m, n) for j = 1:n @@ -54,11 +144,11 @@ function Matrix(A::AbstractToeplitz{T}) where T return Af end -convert(::Type{Matrix}, A::AbstractToeplitz) = Matrix(A) +Base.convert(::Type{Matrix}, A::AbstractToeplitz) = Matrix(A) # Fast application of a general Toeplitz matrix to a column vector via FFT -function mul!( - y::StridedVector, A::AbstractToeplitz, x::StridedVector, α::Number, β::Number +function LinearAlgebra.mul!( + y::AbstractVector, A::AbstractToeplitz, x::AbstractVector, α::Number, β::Number ) m, n = size(A) if length(y) != m @@ -74,7 +164,7 @@ function mul!( # Small case: don't use FFT N = m + n - 1 - if N < 512 + if N < MIN_FFT_LENGTH # Scale/initialize y if iszero(β) fill!(y, 0) @@ -95,8 +185,8 @@ function mul!( return y end -function mul!( - y::StridedVector, A::ToeplitzFactorization, x::StridedVector, α::Number, β::Number +function LinearAlgebra.mul!( + y::AbstractVector, A::ToeplitzFactorization, x::AbstractVector, α::Number, β::Number ) n = length(x) m = length(y) @@ -133,17 +223,26 @@ function mul!( end end end - return y end +function Base.:*(A::ToeplitzFactorization, x::AbstractVector) + T = promote_type(eltype(A), eltype(x)) + y = zeros(T, size(A, 1)) + mul!(y, A, x) +end +function Base.:*(A::ToeplitzFactorization, X::AbstractMatrix) + T = promote_type(eltype(A), eltype(X)) + Y = zeros(T, size(A, 1), size(X, 2)) + mul!(Y, A, X) +end # Application of a general Toeplitz matrix to a general matrix -function mul!( +function LinearAlgebra.mul!( C::StridedMatrix, A::AbstractToeplitz, B::StridedMatrix, α::Number, β::Number ) return mul!(C, factorize(A), B, α, β) end -function mul!( +function LinearAlgebra.mul!( C::StridedMatrix, A::ToeplitzFactorization, B::StridedMatrix, α::Number, β::Number ) l = size(B, 2) @@ -156,27 +255,6 @@ function mul!( return C end -# Left division of a general matrix B by a general Toeplitz matrix A, i.e. the solution x of Ax=B. -function ldiv!(A::AbstractToeplitz, B::StridedMatrix) - if size(A, 1) != size(A, 2) - error("Division: Rectangular case is not supported.") - end - for j = 1:size(B, 2) - ldiv!(A, view(B, :, j)) - end - 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) - copyto!(bb, b) - ldiv!(A, bb) -end - # General Toeplitz matrix """ Toeplitz @@ -220,13 +298,21 @@ end function LinearAlgebra.factorize(A::Toeplitz) T = eltype(A) - m, n = size(A) + n, m = size(A) S = promote_type(T, Complex{Float32}) - tmp = Vector{S}(undef, m + n - 1) + tmp = Vector{S}(undef, n + m - 1) copyto!(tmp, A.vc) - copyto!(tmp, m + 1, Iterators.reverse(A.vr), 1, n - 1) + copyto!(tmp, n + 1, Iterators.reverse(A.vr), 1, m - 1) dft = plan_fft!(tmp) - return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) + return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, n, m) +end +function Base.inv(A::Toeplitz) + n = checksquare(A) + if MIN_FFT_LENGTH > 2n - 1 + inv(Matrix(A)) + else + A \ (one(eltype(A))*I)(n) + end end convert(::Type{AbstractToeplitz{T}}, A::Toeplitz) where {T} = convert(Toeplitz{T}, A) @@ -254,8 +340,8 @@ end function getindex(A::Toeplitz, i::Integer, j::Integer) m = size(A,1) n = size(A,2) - if i > m || j > n - error("BoundsError()") + @boundscheck if i > m || j > n + error(BoundsError("index ($i, $j) out of bounds for A of size $(size(A))")) end if i >= j @@ -293,11 +379,6 @@ function triu(A::Toeplitz, k = 0) return Al end -function ldiv!(A::Toeplitz, b::StridedVector) - preconditioner = factorize(strang(A)) - copyto!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, preconditioner, 1000, 100eps())[1]) -end - # Symmetric """ SymmetricToeplitz @@ -323,7 +404,7 @@ function LinearAlgebra.factorize(A::SymmetricToeplitz) @inbounds tmp[m + 1] = zero(T) copyto!(tmp, m + 2, Iterators.reverse(vc), 1, m - 1) dft = plan_fft!(tmp) - return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) + return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, m, m) end convert(::Type{AbstractToeplitz{T}}, A::SymmetricToeplitz) where {T} = convert(SymmetricToeplitz{T},A) @@ -332,6 +413,8 @@ convert(::Type{SymmetricToeplitz{T}}, A::SymmetricToeplitz) where {T} = Symmetri adjoint(A::SymmetricToeplitz) = SymmetricToeplitz(conj(A.vr), conj(A.vc)) adjoint(A::SymmetricToeplitz{<:Real}) = A transpose(A::SymmetricToeplitz) = A +LinearAlgebra.issymmetric(::SymmetricToeplitz) = true +LinearAlgebra.ishermitian(::SymmetricToeplitz{<:Real}) = true function size(A::SymmetricToeplitz, dim::Int) if 1 <= dim <= 2 @@ -343,8 +426,26 @@ end getindex(A::SymmetricToeplitz, i::Integer, j::Integer) = A.vc[abs(i - j) + 1] -ldiv!(A::SymmetricToeplitz, b::StridedVector) = - copyto!(b, IterativeLinearSolvers.cg(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1]) +function ldiv!(x::AbstractVector, A::SymmetricToeplitz, b::AbstractVector; + isposdef::Bool = false, verbose::Bool = false, + maxiter::Int = DEFAULT_MAXITER, + abstol::Real = max(DEFAULT_TOL, zero(real(eltype(b)))), + reltol::Real = max(DEFAULT_TOL, sqrt(eps(real(eltype(b)))))) + n, m = size(A) + F = n + m - 1 < MIN_FFT_LENGTH ? A : factorize(A) # if we are going to use FFT-based multiplication, factorize once before iterative method + P = factorize(strang(A)) # since this is circulant, we always factorize it + if isposdef + IterativeSolvers.cg!(x, F, b, Pl = P, maxiter = maxiter, + abstol = abstol, reltol = reltol, + verbose = verbose, log = false) + else + P = sqrt(abs(P)) # left and right preconditioner tends to have slightly better conditioning + IterativeSolvers.gmres!(x, F, b, Pl = P, Pr = P', maxiter = maxiter, + abstol = abstol, reltol = reltol, + verbose = verbose, log = false) + # IterativeSolvers.minres!(x, F, b, maxiter = DEFAULT_MAXITER, abstol = abstol, reltol = reltol, verbose = verbose, log = false) # currently does not support preconditioners + end +end # Circulant """ @@ -371,24 +472,26 @@ Create a circulant matrix from the first column of matrix `A`. Circulant(A::AbstractMatrix) = Circulant{eltype(A)}(A) Circulant{T}(A::AbstractMatrix) where {T<:Number} = Circulant{T}(A[:,1]) -const CirculantFactorization{T<:Number} = ToeplitzFactorization{T,Circulant{T}} +const CirculantFactorization{T<:Number} = ToeplitzFactorization{T, Circulant{T}} function LinearAlgebra.factorize(C::Circulant) T = eltype(C) vc = C.vc + n = length(vc) S = promote_type(T, Complex{Float32}) - tmp = Vector{S}(undef, length(vc)) + tmp = Vector{S}(undef, n) copyto!(tmp, vc) dft = plan_fft!(tmp) - return ToeplitzFactorization{T,typeof(C),S,typeof(dft)}(dft * tmp, similar(tmp), dft) + return ToeplitzFactorization{T,typeof(C),S,typeof(dft)}(dft * tmp, zero(tmp), dft, n, n) end convert(::Type{AbstractToeplitz{T}}, A::Circulant) where {T} = convert(Circulant{T}, A) convert(::Type{Circulant{T}}, A::Circulant) where {T} = Circulant(convert(Vector{T}, A.vc)) - function size(C::Circulant, dim::Integer) if 1 <= dim <= 2 return length(C.vc) + elseif dim > 2 + return 1 else error("arraysize: dimension out of range") end @@ -402,6 +505,7 @@ function getindex(C::Circulant, i::Integer, j::Integer) return C.vc[mod(i - j, length(C.vc)) + 1] end +# IDEA: have 3-arg ldiv! for Circulant LinearAlgebra.ldiv!(C::Circulant, b::AbstractVector) = ldiv!(factorize(C), b) function LinearAlgebra.ldiv!(C::CirculantFactorization, b::AbstractVector) n = length(b) @@ -414,46 +518,49 @@ function LinearAlgebra.ldiv!(C::CirculantFactorization, b::AbstractVector) end dft = C.dft @inbounds begin - for i in 1:n + for i in 1:n # IDEA @simd tmp[i] = b[i] end dft * tmp - for i in 1:n + for i in 1:n # IDEA @simd tmp[i] /= vcvr_dft[i] end - dft \ tmp + dft \ tmp # ? T = eltype(C) - for i in 1:n + for i in 1:n # IDEA @simd b[i] = maybereal(T, tmp[i]) end end return b end -function Base.inv(C::Circulant) - F = factorize(C) - vdft = map(inv, F.vcvr_dft) - vc = F.dft \ vdft - return Circulant(maybereal(eltype(C), vc)) -end -function Base.inv(C::CirculantFactorization) - vdft = map(inv, C.vcvr_dft) - return CirculantFactorization(vdft, similar(vdft), C.dft) -end - +""" +``` + strang(A::AbstractMatrix) +``` +Computes circulant preconditioner for `A` according to Gil Strang’s 'Proposal for Toeplitz Matrix Calculations' (Studies in Applied Mathematics, 74, pp. 171–176, 1986.), +for use with Conjugate Gradients for the solution of positive definite systems. +""" function strang(A::AbstractMatrix{T}) where T - n = size(A, 1) + n = checksquare(A) v = Vector{T}(undef, n) n2 = div(n, 2) for i = 1:n if i <= n2 + 1 - v[i] = A[i,1] + v[i] = A[i, 1] else v[i] = A[1, n - i + 2] end end return Circulant(v) end + +""" +``` + chan(A::AbstractMatrix) +``` +Computes circulant preconditioner for `A` according to Raymond H. Chan and Man-Chung Yeung's 'Circulant Preconditioners for Toeplitz Matrices with Positive Continuous Generating Functions'. +""" function chan(A::AbstractMatrix{T}) where T n = size(A, 1) v = Vector{T}(undef, n) @@ -476,10 +583,21 @@ end LinearAlgebra.eigvals(C::Circulant) = eigvals(factorize(C)) LinearAlgebra.eigvals(C::CirculantFactorization) = copy(C.vcvr_dft) -sqrt(C::Circulant) = sqrt(factorize(C)) -function Base.sqrt(C::CirculantFactorization) - vc = C.dft \ sqrt.(C.vcvr_dft) +Base.sqrt(C::Union{Circulant, CirculantFactorization}) = apply(sqrt, C) +Base.abs(C::Union{Circulant, CirculantFactorization}) = apply(abs, C) +Base.inv(C::Union{Circulant, CirculantFactorization}) = apply(inv, C) + +# helper function to apply matrix functions to Circulant matrices +# using efficient diagonalization +function apply(f, C::CirculantFactorization) + vcvr_dft = f.(C.vcvr_dft) + return typeof(C)(vcvr_dft, copy(C.tmp), C.dft, C.n, C.m) +end +function apply(f, C::Circulant) + F = apply(f, factorize(C)) + vc = F.dft \ F.vcvr_dft return Circulant(maybereal(eltype(C), vc)) + end copy(C::Circulant) = Circulant(copy(C.vc)) @@ -520,33 +638,9 @@ function Base.:*(A::CirculantFactorization, B::CirculantFactorization) return Circulant(maybereal(Base.promote_eltype(A, B), vc)) end -Base.:*(A::Adjoint{<:Circulant}, B::Circulant) = factorize(parent(A))' * factorize(B) -Base.:*(A::Adjoint{<:Circulant}, B::CirculantFactorization) = factorize(parent(A))' * B -Base.:*(A::Adjoint{<:CirculantFactorization}, B::Circulant) = A * factorize(B) -function Base.:*(A::Adjoint{<:CirculantFactorization}, B::CirculantFactorization) - C = parent(A) - C_vcvr_dft = C.vcvr_dft - B_vcvr_dft = B.vcvr_dft - m = length(C_vcvr_dft) - n = length(B_vcvr_dft) - if m != n - throw(DimensionMismatch( - "size of matrix A, $(m)x$(m), does not match size of matrix B, $(n)x$(n)" - )) - end - - tmp = map(C_vcvr_dft, B_vcvr_dft) do c, b - conj(c) * b - end - vc = C.dft \ tmp - - return Circulant(maybereal(Base.promote_eltype(C, B), vc)) -end - (*)(scalar::Number, C::Circulant) = Circulant(scalar*C.vc) (*)(C::Circulant,scalar::Number) = Circulant(scalar*C.vc) - # Triangular struct TriangularToeplitz{T<:Number} <: AbstractToeplitz{T} ve::Vector{T} @@ -582,7 +676,7 @@ function LinearAlgebra.factorize(A::TriangularToeplitz) copyto!(tmp, n + 1, Iterators.reverse(ve), 1, n - 1) end dft = plan_fft!(tmp) - return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) + return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft, n, n) end function convert(::Type{Toeplitz}, A::TriangularToeplitz) @@ -645,7 +739,7 @@ function smallinv(A::TriangularToeplitz{T}) where T end b[k] = -tmp/A.ve[1] end - return TriangularToeplitz(b, symbol(A.uplo)) + return TriangularToeplitz(b, Symbol(A.uplo)) end function inv(A::TriangularToeplitz{T}) where T @@ -653,21 +747,15 @@ function inv(A::TriangularToeplitz{T}) where T if n <= 64 return smallinv(A) end - np2 = nextpow2(n) + np2 = nextpow(2, n) if n != np2 - return TriangularToeplitz(inv(TriangularToeplitz([A.ve, zeros(T, np2 - n)], - symbol(A.uplo))).ve[1:n], symbol(A.uplo)) + return TriangularToeplitz(inv(TriangularToeplitz(vcat(A.ve, zeros(T, np2 - n)), + Symbol(A.uplo))).ve[1:n], Symbol(A.uplo)) end nd2 = div(n, 2) - a1 = inv(TriangularToeplitz(A.ve[1:nd2], symbol(A.uplo))).ve - return TriangularToeplitz([a1, -(TriangularToeplitz(a1, symbol(A.uplo)) * - (Toeplitz(A.ve[nd2 + 1:end], A.ve[nd2 + 1:-1:2]) * a1))], symbol(A.uplo)) -end - -# ldiv!(A::TriangularToeplitz,b::StridedVector) = inv(A)*b -function ldiv!(A::TriangularToeplitz, b::StridedVector) - preconditioner = factorize(chan(A)) - copyto!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, preconditioner, 1000, 100eps())[1]) + a1 = inv(TriangularToeplitz(A.ve[1:nd2], Symbol(A.uplo))).ve + return TriangularToeplitz(vcat(a1, -(TriangularToeplitz(a1, Symbol(A.uplo)) * + (Toeplitz(A.ve[nd2 + 1:end], A.ve[nd2 + 1:-1:2]) * a1))), Symbol(A.uplo)) end # extend levinson @@ -775,8 +863,6 @@ convert(::Type{AbstractMatrix{T}}, A::Hankel{T}) where {T<:Number} = A convert(::Type{AbstractMatrix{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A) convert(::Type{Hankel{T}}, A::Hankel) where {T<:Number} = _Hankel(convert(Toeplitz{T}, A.T)) - - # Size size(H::Hankel,k...) = size(H.T,k...) @@ -789,8 +875,8 @@ getindex(A::Hankel, i::Integer, j::Integer) = A.T[i,end-j+1] # Fast application of a general Hankel matrix to a general matrix *(A::Hankel, B::AbstractMatrix) = A.T * flipdim(B, 1) -## BigFloat support +# BigFloat support (*)(A::Toeplitz{T}, b::AbstractVector) where {T<:BigFloat} = irfft( rfft([ A.vc; diff --git a/src/iterativeLinearSolvers.jl b/src/iterativeLinearSolvers.jl deleted file mode 100644 index e856ccd..0000000 --- a/src/iterativeLinearSolvers.jl +++ /dev/null @@ -1,217 +0,0 @@ -module IterativeLinearSolvers - using LinearAlgebra - import LinearAlgebra: Factorization -# Included from https://github.com/andreasnoack/IterativeLinearSolvers.jl -# Eventually, use IterativeSolvers.jl - -Preconditioner{T} = Union{AbstractMatrix{T}, Factorization{T}} - -function cg(A::AbstractMatrix{T}, - x::AbstractVector{T}, - b::AbstractVector{T}, - M::Preconditioner{T}, - max_it::Integer, - tol::Real) where T<:LinearAlgebra.BlasReal -# -- Iterative template routine -- -# Univ. of Tennessee and Oak Ridge National Laboratory -# October 1, 1993 -# Details of this algorithm are described in "Templates for the -# Solution of Linear Systems: Building Blocks for Iterative -# Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra, -# Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications, -# 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps). -# -# [x, error, iter, flag] = cg(A, x, b, M, max_it, tol) -# -# cg.m solves the symmetric positive definite linear system Ax=b -# using the Conjugate Gradient method with preconditioning. -# -# input A REAL symmetric positive definite matrix -# x REAL initial guess vector -# b REAL right hand side vector -# M REAL preconditioner matrix -# max_it INTEGER maximum number of iterations -# tol REAL error tolerance -# -# output x REAL solution vector -# error REAL error norm -# iter INTEGER number of iterations performed -# flag INTEGER: 0 = solution found to tolerance -# 1 = no convergence given max_it - - n = length(b) - flag = 0 # initialization - iter = 0 - - bnrm2 = norm(b) - if bnrm2 == 0.0 bnrm2 = one(T) end - - local ρ₁ - z = zeros(T, n) - q = zeros(T, n) - p = zeros(T, n) - # r = copy(b) - # mul!(r,A,x,-one(T),one(T)) - r = b - A*x - error = norm(r)/bnrm2 - if error < tol - return - end - - for iter_inner = 1:max_it # begin iteration - iter = iter_inner - - z[:] = r - ldiv!(M, z) - # z[:] = M\r - ρ = dot(r,z) - - if iter > 1 # direction vector - β = ρ/ρ₁ - for l = 1:n - p[l] = z[l] + β*p[l] - end - else - p[:] = z - end - - # mul!(q,A,p,one(T),zero(T)) - q[:] = A*p - α = ρ / dot(p,q) - for l = 1:n - x[l] += α*p[l] # update approximation vector - r[l] -= α*q[l] # compute residual - end - - error = norm(r)/bnrm2 # check convergence - if error <= tol - break - end - - ρ₁ = ρ - - end - - if error > tol flag = 1 end # no convergence - return x, error, iter, flag -end - -function cgs(A::AbstractMatrix{T}, - x::AbstractVector{T}, - b::AbstractVector{T}, - M::Preconditioner{T}, - max_it::Integer, - tol::Real) where T<:Number -# -- Iterative template routine -- -# Univ. of Tennessee and Oak Ridge National Laboratory -# October 1, 1993 -# Details of this algorithm are described in "Templates for the -# Solution of Linear Systems: Building Blocks for Iterative -# Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra, -# Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications, -# 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps). -# -# [x, error, iter, flag] = cgs(A, x, b, M, max_it, tol) -# -# cgs.m solves the linear system Ax=b using the -# Conjugate Gradient Squared Method with preconditioning. -# -# input A REAL matrix -# x REAL initial guess vector -# b REAL right hand side vector -# M REAL preconditioner -# max_it INTEGER maximum number of iterations -# tol REAL error tolerance -# -# output x REAL solution vector -# error REAL error norm -# iter INTEGER number of iterations performed -# flag INTEGER: 0 = solution found to tolerance -# 1 = no convergence given max_it - - iter = 0 # initialization - flag = 0 - - n = length(b) - bnrm2 = norm(b) - 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) - v̂ = zeros(T, n) - ρ = zero(T) - ρ₁ = ρ - r = copy(b) - mul!(r, A, x, -one(T), one(T)) - # r = b - A*x - error = norm(r)/bnrm2 - - if error < tol - return x, error, iter, flag - end - - r_tld = copy(r) - - for iter_inner = 1:max_it # begin iteration - iter = iter_inner - - ρ = dot(r_tld, r) - if ρ == 0 - break - end - - if iter > 1 # direction vectors - β = ρ/ρ₁ - for l = 1:n - u[l] = r[l] + β*q[l] - p[l] = u[l] + β*(q[l] + β*p[l]) - end - else - u[:] = r - p[:] = u - end - - p̂[:] = p - ldiv!(M, p̂) - # p̂[:] = M\p - mul!(v̂, A, p̂, one(T), zero(T)) # adjusting scalars - # v̂[:] = A*p̂ - α = ρ/dot(r_tld, v̂) - for l = 1:n - q[l] = u[l] - α*v̂[l] - û[l] = u[l] + q[l] - end - ldiv!(M, û) - # û[:] = M\û - - for l = 1:n - x[l] += α*û[l] # update approximation - end - - mul!(r, A, û, -α, one(T)) - # r[:] -= α*(A*û) - error = norm(r)/bnrm2 # check convergence - if error <= tol - break - end - - ρ₁ = ρ - - end - - if error <= tol # converged - flag = 0 - elseif ρ == 0 # breakdown - flag = -1 - else # no convergence - flag = 1 - end - return x, error, iter, flag -end - -end diff --git a/test/Project.toml b/test/Project.toml index a9a5921..692d529 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -10,4 +10,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" AbstractFFTs = "0.4, 0.5, 1" FFTW = "1" StatsBase = "0.32, 0.33" -julia = "1" +julia = "1.6" diff --git a/test/runtests.jl b/test/runtests.jl index 69f2108..d548b4a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,56 +12,224 @@ using ToeplitzMatrices, StatsBase, Test, LinearAlgebra using Base: copyto! using FFTW: fft -ns = 101 -nl = 2000 +const atol = 1e-6 +const rtol = 1e-6 +function isapprox_helper(x, y, verbose::Bool = true) + b = isapprox(x, y, atol = atol, rtol = rtol) + if ~b && verbose + println("x !≈ y") + println("norm(x-y): ", norm(x-y)) + println("norm(x-y)/norm(y): ", norm(x-y)/norm(y)) + end + return b +end + +ns = 17 +nl = 33 +k = 1 -xs = randn(ns, 5) -xl = randn(nl, 5) -vc = LinRange(1,3,3) # for testing with AbstractVector +xs = randn(ns, k) +vc = LinRange(1, 3, 3) # for testing with AbstractVector vv = Vector(vc) vr = [1, 5.] -cases = [ - (Toeplitz(0.9.^(0:ns-1), 0.4.^(0:ns-1)), - Toeplitz(0.9.^(0:nl-1), 0.4.^(0:nl-1)), - "Real general square"), - (Toeplitz(complex(0.9.^(0:ns-1)), complex(0.4.^(0:ns-1))), - Toeplitz(complex(0.9.^(0:nl-1)), complex(0.4.^(0:nl-1))), - "Complex general square"), - (Circulant(0.9.^(0:ns - 1)), - Circulant(0.9.^(0:nl - 1)), - "Real circulant"), - (Circulant(complex(0.9.^(0:ns - 1))), - Circulant(complex(0.9.^(0:nl - 1))), - "Complex circulant"), - (TriangularToeplitz(0.9.^(0:ns - 1), :U), - TriangularToeplitz(0.9.^(0:nl - 1), :U), - "Real upper triangular"), - (TriangularToeplitz(complex(0.9.^(0:ns - 1)), :U), - TriangularToeplitz(complex(0.9.^(0:nl - 1)), :U), - "Complex upper triangular"), - (TriangularToeplitz(0.9.^(0:ns - 1), :L), - TriangularToeplitz(0.9.^(0:nl - 1), :L), - "Real lower triangular"), - (TriangularToeplitz(complex(0.9.^(0:ns - 1)), :L), - TriangularToeplitz(complex(0.9.^(0:nl - 1)), :L), - "Complex lower triangular"), -] - -for (As, Al, st) in cases - @testset "Toeplitz: $st" begin - @test As * xs ≈ Matrix(As) * xs - @test As'* xs ≈ Matrix(As)' * xs - @test Al * xl ≈ Matrix(Al) * xl - @test Al'* xl ≈ Matrix(Al)' * xl - @test [As[n] for n in 1:length(As)] == vec(As) - @test [Al[n] for n in 1:length(Al)] == vec(Al) - @test ldiv!(As, LinearAlgebra.copy_oftype(xs, eltype(As))) ≈ Matrix(As) \ xs - @test ldiv!(Al, LinearAlgebra.copy_oftype(xl, eltype(Al))) ≈ Matrix(Al) \ xl - @test Matrix(As') == Matrix(As)' - @test Matrix(transpose(As)) == transpose(Matrix(As)) - end -end +diag_val = 2 # elevating diagonal encourages well-conditioned spectrum +sizes = (17, 33, 513, 1024) +for n in sizes + @testset "n = $n" begin + exp_95 = 0.9.^(0:n-1) + exp_40 = 0.4.^(0:n-1) + exp_40[1] = exp_95[1] = diag_val + + rs1 = randn(n) / n # this scaling encourages good conditioning and diagonal dominance + rs2 = randn(n) / n + rs1[1] = rs2[1] = diag_val # making sure the first elements have the same value + + cases = [ + (Toeplitz(exp_95, exp_40), + "Real square"), + (Toeplitz(complex.(exp_95), complex.(exp_40)), + "Complex square"), + (Toeplitz(rs1, rs2), + "Real random square"), + (Toeplitz(complex(rs1), complex(rs2)), + "Complex random square"), + (Circulant(exp_95), + "Real circulant"), + (Circulant(complex.(exp_95)), + "Complex circulant"), + (Circulant(rs1), + "Real random circulant"), + (Circulant(complex.(rs1)), + "Complex random circulant"), + (TriangularToeplitz(exp_95, :U), + "Real upper triangular"), + (TriangularToeplitz(complex.(exp_95), :U), + "Complex upper triangular"), + (TriangularToeplitz(exp_95, :L), + "Real lower triangular"), + (TriangularToeplitz(complex.(exp_95), :L), + "Complex lower triangular"), + ] + + for (A, st) in cases + x = randn(eltype(A), n) + X = randn(eltype(A), n, k) # for multiplication and solver tests + @testset "Toeplitz: $st" begin + M = Matrix(A) + @test A * x ≈ M * x + @test A'* x ≈ M' * x + @test A * X ≈ M * X + @test A'* X ≈ M' * X + @test [A[n] for n in 1:length(A)] == vec(A) + @test isapprox_helper(A \ x, M \ x) + @test isapprox_helper(A \ X, M \ X) + @test isapprox_helper(ldiv!(A, LinearAlgebra.copy_oftype(X, eltype(A))), M \ X) + @test Matrix(A') == M' + @test adjoint(factorize(A)) * X ≈ M' * X + @test Matrix(transpose(A)) == transpose(M) + @test size(A) == size(M) + @test size(A, 3) == size(M, 3) + F = factorize(A) + @test size(F) == size(M) + @test size(F, 3) == size(M, 3) + @test A * A ≈ M * M + @test isapprox_helper(inv(A), inv(M)) + end # testset matrices + end # loop over matrices + + @testset "Rectangular" begin + for m in sizes + @testset "m = $m" begin + exp_40_m = 0.4.^(0:m-1) + exp_40_m[1] = diag_val + rectangular_cases = [ + (Toeplitz(exp_95, exp_40_m), + "Real general rectangular"), + (Toeplitz(complex.(exp_95), complex.(exp_40_m)), + "Complex general rectangular"), + ] + + for (A, st) in rectangular_cases + xn, xm = randn(eltype(A), n, k), randn(eltype(A), m, k) + @testset "Toeplitz: $st" begin + M = Matrix(A) + @test A * xm ≈ M * xm + @test A'* xn ≈ M' * xn + @test [A[n] for n in 1:length(A)] == vec(A) + @test isapprox(A \ xn, M \ xn, atol = 1e-4, rtol = 1e-4) + @test Matrix(A') == M' + @test Matrix(transpose(A)) == transpose(M) + end + end + end # testset m + end # loop m + end + + @testset "Symmetric Toeplitz" begin + Xs = randn(n, k) + xs = randn(n) + + As = SymmetricToeplitz(exp_95) + Ms = Matrix(As) + @test isapprox_helper(As * xs, Ms * xs) + @test isapprox_helper(As * Xs, Ms * Xs) + @test isapprox_helper(ldiv!(As, copy(xs)), Ms \ xs) + @test isapprox_helper(ldiv!(As, copy(Xs)), Ms \ Xs) + @test StatsBase.levinson(As, xs) ≈ Ms \ xs # this should be exact to numerical precision given good conditioning + @test Matrix(As') ≈ Ms' + @test Matrix(transpose(As)) ≈ transpose(Ms) + y_cg = ldiv!(zero(xs), As, xs, isposdef = true) + @test isapprox_helper(y_cg, As \ xs) # testing cg-based solve + + Ab = SymmetricToeplitz(rs1) # this is still positive definite + Mb = Matrix(Ab) + @test isapprox_helper(Ab * xs, Mb * xs) + @test isapprox_helper(Ab * Xs, Mb * Xs) + @test isapprox_helper(ldiv!(Ab, copy(xs)), Mb \ xs) + @test isapprox_helper(ldiv!(Ab, copy(Xs)), Mb \ Xs) + @test StatsBase.levinson(Ab, xs) ≈ Mb \ xs + @test Matrix(Ab') ≈ Mb' + @test Matrix(transpose(Ab)) ≈ transpose(Mb) + + @test Matrix(SymmetricToeplitz(vc)) == Matrix(SymmetricToeplitz(vv)) + end + + @testset "Hankel" begin + @testset "Real square" begin + H = Hankel([1.0,2,3,4,5],[5.0,6,7,8,9]) + @test Matrix(H) == [1 2 3 4 5; + 2 3 4 5 6; + 3 4 5 6 7; + 4 5 6 7 8; + 5 6 7 8 9] + + @test convert(Hankel{Float64}, H) == H + @test convert(AbstractMatrix{Float64}, H) == H + @test convert(AbstractArray{Float64}, H) == H + + @test H[2,2] == 3 + @test H[7] == 3 + @test diag(H) == [1,3,5,7,9] + + x = ones(5) + @test Matrix(H)*x ≈ H*x + + Hs = Hankel(0.9.^(n-1:-1:0), 0.4.^(0:n-1)) + xs = randn(n) + @test Hs * xs[:,1] ≈ Matrix(Hs) * xs[:,1] + @test Hs * xs ≈ Matrix(Hs) * xs + @test Matrix(Hankel(reverse(vc),vr)) == Matrix(Hankel(reverse(vv),vr)) + end + + @testset "Complex square" begin + H = Hankel(complex([1.0,2,3,4,5]), complex([5.0,6,7,8,0])) + x = ones(5) + @test Matrix(H)*x ≈ H*x + + Hs = Hankel(complex(0.9.^(n-1:-1:0)), complex(0.4.^(0:n-1))) + xs = randn(n) + @test Hs * xs[:,1] ≈ Matrix(Hs) * xs[:,1] + @test Hs * xs ≈ Matrix(Hs) * xs + end + + for m in sizes + xm = randn(m) + exp_95_n = reverse(exp_95) + exp_40_m = 0.4.^(0:m-1) + exp_40_m[1] = diag_val + @testset "m = $m" begin + @testset "Real rectangular" begin + Hs = Hankel(exp_95_n, exp_40_m) + @test Hs * xm[:,1] ≈ Matrix(Hs) * xm[:,1] + @test Hs * xm ≈ Matrix(Hs) * xm + end + + @testset "Complex rectangular" begin + Hs = Hankel(complex.(exp_95_n), complex.(exp_40_m)) + @test Hs * xm[:,1] ≈ Matrix(Hs) * xm[:,1] + @test Hs * xm ≈ Matrix(Hs) * xm + end + end # testset m + end # loop m + + @testset "Convert" begin + H = Hankel([1.0,2,3,4,5],[5.0,6,7,8,0]) + @test Hankel(H) == Hankel{Float64}(H) == H + @test convert(Hankel,H) == convert(Hankel{Float64},H) == + convert(AbstractArray,H) == convert(AbstractArray{Float64},H) == H + @test convert(Array, H) ≈ Matrix(H) + + A = [1.0 2; 3 4] + @test Hankel(A) == [1 3; 3 4] + T = Toeplitz([1.0,2,3,4,5],[1.0,6,7,8,0]) + @test Hankel(T) == Hankel([1.0,2,3,4,5],[5.0,4,3,2,1]) + @test Hankel(T) ≠ ToeplitzMatrices._Hankel(T) + end + end + + end # testset n +end # loop over n @testset "Mixed types" begin @test eltype(Toeplitz([1, 2], [1, 2])) === Int @@ -72,105 +240,6 @@ end @test Matrix(TriangularToeplitz(vc,:U)) == Matrix(TriangularToeplitz(vv,:U)) end -@testset "Real general rectangular" begin - Ar1 = Toeplitz(0.9.^(0:nl-1), 0.4.^(0:ns-1)) - Ar2 = Toeplitz(0.9.^(0:ns-1), 0.4.^(0:nl-1)) - @test Ar1 * xs ≈ Matrix(Ar1) * xs - @test Ar2 * xl ≈ Matrix(Ar2) * xl -end - -@testset "Complex general rectangular" begin - Ar1 = Toeplitz(complex(0.9.^(0:nl-1)), complex(0.4.^(0:ns-1))) - Ar2 = Toeplitz(complex(0.9.^(0:ns-1)), complex(0.4.^(0:nl-1))) - @test Ar1 * xs ≈ Matrix(Ar1) * xs - @test Ar2 * xl ≈ Matrix(Ar2) * xl -end - -@testset "Symmetric Toeplitz" begin - As = SymmetricToeplitz(0.9.^(0:ns-1)) - Ab = SymmetricToeplitz(abs.(randn(ns))) - Al = SymmetricToeplitz(0.9.^(0:nl-1)) - @test As * xs ≈ Matrix(As) * xs - @test Ab * xs ≈ Matrix(Ab) * xs - @test Al * xl ≈ Matrix(Al) * xl - @test ldiv!(As, copy(xs)) ≈ Matrix(As) \ xs - @test ldiv!(Ab, copy(xs)) ≈ Matrix(Ab) \ xs - @test ldiv!(Al, copy(xl)) ≈ Matrix(Al) \ xl - @test StatsBase.levinson(As, xs) ≈ Matrix(As) \ xs - @test StatsBase.levinson(Ab, xs) ≈ Matrix(Ab) \ xs - @test Matrix(SymmetricToeplitz(vc)) == Matrix(SymmetricToeplitz(vv)) -end - -@testset "Hankel" begin - @testset "Real square" begin - H = Hankel([1.0,2,3,4,5],[5.0,6,7,8,9]) - @test Matrix(H) == [1 2 3 4 5; - 2 3 4 5 6; - 3 4 5 6 7; - 4 5 6 7 8; - 5 6 7 8 9] - - @test convert(Hankel{Float64}, H) == H - @test convert(AbstractMatrix{Float64}, H) == H - @test convert(AbstractArray{Float64}, H) == H - - @test H[2,2] == 3 - @test H[7] == 3 - @test diag(H) == [1,3,5,7,9] - - x = ones(5) - @test Matrix(H)*x ≈ H*x - - Hs = Hankel(0.9.^(ns-1:-1:0), 0.4.^(0:ns-1)) - Hl = Hankel(0.9.^(nl-1:-1:0), 0.4.^(0:nl-1)) - @test Hs * xs[:,1] ≈ Matrix(Hs) * xs[:,1] - @test Hs * xs ≈ Matrix(Hs) * xs - @test Hl * xl ≈ Matrix(Hl) * xl - @test Matrix(Hankel(reverse(vc),vr)) == Matrix(Hankel(reverse(vv),vr)) - end - - @testset "Complex square" begin - H = Hankel(complex([1.0,2,3,4,5]), complex([5.0,6,7,8,0])) - x = ones(5) - @test Matrix(H)*x ≈ H*x - - Hs = Hankel(complex(0.9.^(ns-1:-1:0)), complex(0.4.^(0:ns-1))) - Hl = Hankel(complex(0.9.^(nl-1:-1:0)), complex(0.4.^(0:nl-1))) - @test Hs * xs[:,1] ≈ Matrix(Hs) * xs[:,1] - @test Hs * xs ≈ Matrix(Hs) * xs - @test Hl * xl ≈ Matrix(Hl) * xl - end - - @testset "Real rectangular" begin - Hs = Hankel(0.9.^(ns-1:-1:0), 0.4.^(0:nl-1)) - Hl = Hankel(0.9.^(nl-1:-1:0), 0.4.^(0:ns-1)) - @test Hs * xl[:,1] ≈ Matrix(Hs) * xl[:,1] - @test Hs * xl ≈ Matrix(Hs) * xl - @test Hl * xs ≈ Matrix(Hl) * xs - end - - @testset "Complex rectangular" begin - Hs = Hankel(complex(0.9.^(ns-1:-1:0)), complex(0.4.^(0:nl-1))) - Hl = Hankel(complex(0.9.^(nl-1:-1:0)), complex(0.4.^(0:ns-1))) - @test Hs * xl[:,1] ≈ Matrix(Hs) * xl[:,1] - @test Hs * xl ≈ Matrix(Hs) * xl - @test Hl * xs ≈ Matrix(Hl) * xs - end - - @testset "Convert" begin - H = Hankel([1.0,2,3,4,5],[5.0,6,7,8,0]) - @test Hankel(H) == Hankel{Float64}(H) == H - @test convert(Hankel,H) == convert(Hankel{Float64},H) == - convert(AbstractArray,H) == convert(AbstractArray{Float64},H) == H - - A = [1.0 2; 3 4] - @test Hankel(A) == [1 3; 3 4] - T = Toeplitz([1.0,2,3,4,5],[1.0,6,7,8,0]) - @test Hankel(T) == Hankel([1.0,2,3,4,5],[5.0,4,3,2,1]) - @test Hankel(T) ≠ ToeplitzMatrices._Hankel(T) - end -end - @testset "Convert" begin T = Toeplitz(ones(2),ones(2)) @@ -196,13 +265,14 @@ end @test isa(convert(ToeplitzMatrices.AbstractToeplitz{ComplexF64},T),Circulant{ComplexF64}) @test isa(convert(ToeplitzMatrices.Circulant{ComplexF64},T),Circulant{ComplexF64}) - T = TriangularToeplitz(ones(2),:U) + T = TriangularToeplitz(ones(2), :U) @test isa(convert(Matrix{ComplexF64},T),Matrix{ComplexF64}) @test isa(convert(AbstractMatrix{ComplexF64},T),TriangularToeplitz{ComplexF64}) @test isa(convert(AbstractArray{ComplexF64},T),TriangularToeplitz{ComplexF64}) @test isa(convert(ToeplitzMatrices.AbstractToeplitz{ComplexF64},T),TriangularToeplitz{ComplexF64}) @test isa(convert(ToeplitzMatrices.TriangularToeplitz{ComplexF64},T),TriangularToeplitz{ComplexF64}) + @test isa(convert(Toeplitz, T), Toeplitz) T = Hankel(ones(2),ones(2)) @@ -316,6 +386,13 @@ end @test C isa Circulant @test C*C ≈ C3 + lambda_abs = eigvals(abs(C1)) + abs_lambda = abs.(eigvals(C1)) + @test all(<(1e-12), imag(lambda_abs)) + lambda_abs = real.(lambda_abs) + @test all(>(0), lambda_abs) + @test sort!(lambda_abs) ≈ sort!(abs_lambda) + C = copy(C1) @test C isa Circulant C2 = similar(C1)