Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Improved tridiagonal and new Woodbury matrix algebra

- Converts Tridiagonal to being Lapack-compatible
- Wraps Lapack's major tridiagonal routines
- Adds a new "Woodbury" type for solving equations using the Woodbury matrix identity
- Moves the functionality into appropriate functions in base
- Provides a set of tests for this functionality
  • Loading branch information...
commit cbaad366a2a58cc22bfd31838b99327c8329098f 1 parent 2a011e3
@timholy timholy authored
View
6 base/export.jl
@@ -81,15 +81,20 @@ export
SubOrDArray,
SubString,
TransformedString,
+ Tridiagonal,
VecOrMat,
Vector,
VersionNumber,
WeakKeyDict,
+ Woodbury,
Zip,
Stat,
Factorization,
Cholesky,
LU,
+ LUTridiagonal,
+ LDLT,
+ LDLTTridiagonal,
QR,
QRP,
@@ -634,6 +639,7 @@ export
randsym,
rank,
rref,
+ solve,
svd,
svdvals,
trace,
View
45 base/factorizations.jl
@@ -281,3 +281,48 @@ end
##ToDo: Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly
## Add rcond methods for Cholesky, LU, QR and QRP types
## Lower priority: Add LQ, QL and RQ factorizations
+
+
+#### Factorizations for Tridiagonal ####
+type LDLTTridiagonal{T} <: Factorization{T}
+ D::Vector{T}
+ E::Vector{T}
+end
+function LDLTTridiagonal{T<:LapackScalar}(A::Tridiagonal{T})
+ D = copy(A.d)
+ E = copy(A.dl)
+ _jl_lapack_pttrf(D, E)
+ LDLTTridiagonal(D, E)
+end
+LDLT(A::Tridiagonal) = LDLTTridiagonal(A)
+
+(\){T<:LapackScalar}(C::LDLTTridiagonal{T}, B::StridedVecOrMat{T}) =
+ _jl_lapack_pttrs(C.D, C.E, copy(B))
+
+type LUTridiagonal{T} <: Factorization{T}
+ lu::Tridiagonal{T}
+ ipiv::Vector{Int32}
+ function LUTridiagonal(lu::Tridiagonal{T}, ipiv::Vector{Int32})
+ m, n = size(lu)
+ m == numel(ipiv) ? new(lu, ipiv) : error("LU: dimension mismatch")
+ end
+end
+show(io, lu::LUTridiagonal) = print(io, "LU decomposition of ", summary(lu.lu))
+
+function LU{T<:LapackScalar}(A::Tridiagonal{T})
+ lu, ipiv = _jl_lapack_gttrf(copy(A))
+ LUTridiagonal{T}(lu, ipiv)
+end
+
+function lu(A::Tridiagonal)
+ error("lu(A) is not defined when A is Tridiagonal. Use LU(A) instead.")
+end
+
+function det(lu::LUTridiagonal)
+ prod(lu.lu.d) * (bool(sum(lu.ipiv .!= 1:n) % 2) ? -1 : 1)
+end
+
+det(A::Tridiagonal) = det(LU(A))
+
+(\){T<:LapackScalar}(lu::LUTridiagonal{T}, B::StridedVecOrMat{T}) =
+ _jl_lapack_gttrs('N', lu.lu, lu.ipiv, copy(B))
View
4 base/linalg.jl
@@ -1,4 +1,6 @@
-## linalg.jl: Basic Linear Algebra interface specifications ##
+## linalg.jl: Basic Linear Algebra interface specifications and
+## specialized matrix types
+
#
# This file mostly contains commented functions which are supposed
# to be defined in type-specific linalg_<type>.jl files.
View
184 base/linalg_lapack.jl
@@ -885,3 +885,187 @@ end
expm{T<:Union(Float32,Float64,Complex64,Complex128)}(A::StridedMatrix{T}) = expm!(copy(A))
expm{T<:Integer}(A::StridedMatrix{T}) = expm!(float(A))
+
+#### Tridiagonal matrix routines ####
+function \{T<:LapackScalar}(M::Tridiagonal{T}, rhs::StridedVecOrMat{T})
+ if stride(rhs, 1) == 1
+ x = copy(rhs)
+ Mc = copy(M)
+ Mlu, x = _jl_lapack_gtsv(Mc, x)
+ return x
+ end
+ solve(M, rhs) # use the Julia "fallback"
+end
+
+eig(M::Tridiagonal) = _jl_lapack_stev('V', copy(M))
+
+# Decompositions
+for (gttrf, pttrf, elty) in
+ ((:dgttrf_,:dpttrf_,:Float64),
+ (:sgttrf_,:spttrf_,:Float32),
+ (:zgttrf_,:zpttrf_,:Complex128),
+ (:cgttrf_,:cpttrf_,:Complex64))
+ @eval begin
+ function _jl_lapack_gttrf(M::Tridiagonal{$elty})
+ info = zero(Int32)
+ n = int32(length(M.d))
+ ipiv = Array(Int32, n)
+ ccall(dlsym(_jl_liblapack, $string(gttrf)),
+ Void,
+ (Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+ Ptr{Int32}, Ptr{Int32}),
+ &n, M.dl, M.d, M.du, M.dutmp, ipiv, &info)
+ if info != 0 throw(LapackException(info)) end
+ M, ipiv
+ end
+ function _jl_lapack_pttrf(D::Vector{$elty}, E::Vector{$elty})
+ info = zero(Int32)
+ n = int32(length(D))
+ if length(E) != n-1
+ error("subdiagonal must be one element shorter than diagonal")
+ end
+ ccall(dlsym(_jl_liblapack, $string(pttrf)),
+ Void,
+ (Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{Int32}),
+ &n, D, E, &info)
+ if info != 0 throw(LapackException(info)) end
+ D, E
+ end
+ end
+end
+# Direct solvers
+for (gtsv, ptsv, elty) in
+ ((:dgtsv_,:dptsv_,:Float64),
+ (:sgtsv_,:sptsv,:Float32),
+ (:zgtsv_,:zptsv,:Complex128),
+ (:cgtsv_,:cptsv,:Complex64))
+ @eval begin
+ function _jl_lapack_gtsv(M::Tridiagonal{$elty}, B::StridedVecOrMat{$elty})
+ if stride(B,1) != 1
+ error("_jl_lapack_gtsv: matrix columns must have contiguous elements");
+ end
+ info = zero(Int32)
+ n = int32(length(M.d))
+ nrhs = int32(size(B, 2))
+ ldb = int32(stride(B, 2))
+ ccall(dlsym(_jl_liblapack, $string(gtsv)),
+ Void,
+ (Ptr{Int32}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+ Ptr{Int32}, Ptr{Int32}),
+ &n, &nrhs, M.dl, M.d, M.du, B, &ldb, &info)
+ if info != 0 throw(LapackException(info)) end
+ M, B
+ end
+ function _jl_lapack_ptsv(M::Tridiagonal{$elty}, B::StridedVecOrMat{$elty})
+ if stride(B,1) != 1
+ error("_jl_lapack_ptsv: matrix columns must have contiguous elements");
+ end
+ info = zero(Int32)
+ n = int32(length(M.d))
+ nrhs = int32(size(B, 2))
+ ldb = int32(stride(B, 2))
+ ccall(dlsym(_jl_liblapack, $string(ptsv)),
+ Void,
+ (Ptr{Int32}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+ Ptr{Int32}, Ptr{Int32}),
+ &n, &nrhs, M.d, M.dl, B, &ldb, &info)
+ if info != 0 throw(LapackException(info)) end
+ M, B
+ end
+ end
+end
+# Solvers using decompositions
+for (gttrs, pttrs, elty) in
+ ((:dgttrs_,:dpttrs_,:Float64),
+ (:sgttrs_,:spttrs,:Float32),
+ (:zgttrs_,:zpttrs,:Complex128),
+ (:cgttrs_,:cpttrs,:Complex64))
+ @eval begin
+ function _jl_lapack_gttrs(trans::LapackChar, M::Tridiagonal{$elty}, ipiv::Vector{Int32}, B::StridedVecOrMat{$elty})
+ if stride(B,1) != 1
+ error("_jl_lapack_gttrs: matrix columns must have contiguous elements");
+ end
+ info = zero(Int32)
+ n = int32(length(M.d))
+ nrhs = int32(size(B, 2))
+ ldb = int32(stride(B, 2))
+ ccall(dlsym(_jl_liblapack, $string(gttrs)),
+ Void,
+ (Ptr{Uint8}, Ptr{Int32}, Ptr{Int32},
+ Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+ Ptr{Int32}, Ptr{$elty}, Ptr{Int32}, Ptr{Int32}),
+ &trans, &n, &nrhs, M.dl, M.d, M.du, M.dutmp, ipiv, B, &ldb, &info)
+ if info != 0 throw(LapackException(info)) end
+ B
+ end
+ function _jl_lapack_pttrs(D::Vector{$elty}, E::Vector{$elty}, B::StridedVecOrMat{$elty})
+ if stride(B,1) != 1
+ error("_jl_lapack_pttrs: matrix columns must have contiguous elements");
+ end
+ info = zero(Int32)
+ n = int32(length(D))
+ if length(E) != n-1
+ error("subdiagonal must be one element shorter than diagonal")
+ end
+ nrhs = int32(size(B, 2))
+ ldb = int32(stride(B, 2))
+ ccall(dlsym(_jl_liblapack, $string(pttrs)),
+ Void,
+ (Ptr{Int32}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+ Ptr{Int32}, Ptr{Int32}),
+ &n, &nrhs, D, E, B, &ldb, &info)
+ if info != 0 throw(LapackException(info)) end
+ B
+ end
+ end
+end
+# Eigenvalue-eigenvector (symmetric only)
+for (stev, elty) in
+ ((:dstev_,:Float64),
+ (:sstev_,:Float32),
+ (:zstev_,:Complex128),
+ (:cstev_,:Complex64))
+ @eval begin
+ function _jl_lapack_stev(Z::Array, M::Tridiagonal{$elty})
+ n = int32(length(M.d))
+ if isempty(Z)
+ job = 'N'
+ ldz = 1
+ work = Array($elty, 0)
+ Ztmp = work
+ else
+ if stride(Z,1) != 1
+ error("_jl_lapack_stev: eigenvector matrix columns must have contiguous elements");
+ end
+ if size(Z, 1) != n
+ error("_jl_lapack_stev: eigenvector matrix columns are not of the correct size")
+ end
+ Ztmp = Z
+ job = 'V'
+ ldz = int32(stride(Z, 2))
+ work = Array($elty, max(1, 2*n-2))
+ end
+ info = zero(Int32)
+ ccall(dlsym(_jl_liblapack, $string(stev)),
+ Void,
+ (Ptr{Uint8}, Ptr{Int32},
+ Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+ Ptr{Int32}, Ptr{$elty}, Ptr{Int32}),
+ &job, &n, M.d, M.dl, Ztmp, &ldz, work, &info)
+ if info != 0 throw(LapackException(info)) end
+ M.d
+ end
+ end
+end
+function _jl_lapack_stev(job::LapackChar, M::Tridiagonal)
+ if job == 'N' || job == 'n'
+ Z = []
+ elseif job == 'V' || job == 'v'
+ n = length(M.d)
+ Z = Array(eltype(M), n, n)
+ else
+ error("Job type not recognized")
+ end
+ D = _jl_lapack_stev(Z, M)
+ return D, Z
+end
View
299 base/linalg_specialized.jl
@@ -0,0 +1,299 @@
+#### Specialized matrix types ####
+
+# Some of these also have important routines defined in factorizations.jl,
+# linalg_lapack.jl, etc.
+
+#### Tridiagonal matrices ####
+type Tridiagonal{T} <: AbstractMatrix{T}
+ dl::Vector{T} # sub-diagonal
+ d::Vector{T} # diagonal
+ du::Vector{T} # sup-diagonal
+ dutmp::Vector{T} # scratch space for vector RHS solver, sup-diagonal
+ rhstmp::Vector{T} # scratch space, rhs
+
+ function Tridiagonal(N::Int)
+ dutmp = Array(T, N-1)
+ rhstmp = Array(T, N)
+ new(dutmp, rhstmp, dutmp, dutmp, rhstmp) # first three will be overwritten
+ end
+end
+function Tridiagonal{T<:Number}(dl::Vector{T}, d::Vector{T}, du::Vector{T})
+ N = length(d)
+ if length(dl) != N-1 || length(du) != N-1
+ error("The sub- and super-diagonals must have length N-1")
+ end
+ M = Tridiagonal{T}(N)
+ M.dl = copy(dl)
+ M.d = copy(d)
+ M.du = copy(du)
+ return M
+end
+function Tridiagonal{Tl<:Number, Td<:Number, Tu<:Number}(dl::Vector{Tl}, d::Vector{Td}, du::Vector{Tu})
+ R = promote(Tl, Td, Tu)
+ Tridiagonal(convert(Vector{R}, dl), convert(Vector{R}, d), convert(Vector{R}, du))
+end
+size(M::Tridiagonal) = (length(M.d), length(M.d))
+function show(io, M::Tridiagonal)
+ println(io, summary(M), ":")
+ print(io, " sub: ")
+ print_matrix(io, (M.dl)')
+ print(io, "\ndiag: ")
+ print_matrix(io, (M.d)')
+ print(io, "\n sup: ")
+ print_matrix(io, (M.du)')
+end
+full{T}(M::Tridiagonal{T}) = convert(Matrix{T}, M)
+function convert{T}(::Type{Matrix{T}}, M::Tridiagonal{T})
+ A = zeros(T, size(M))
+ for i = 1:length(M.d)
+ A[i,i] = M.d[i]
+ end
+ for i = 1:length(M.d)-1
+ A[i+1,i] = M.dl[i]
+ A[i,i+1] = M.du[i]
+ end
+ return A
+end
+function similar(M::Tridiagonal, T, dims::Dims)
+ if length(dims) != 2 || dims[1] != dims[2]
+ error("Tridiagonal matrices must be square")
+ end
+ return Tridiagonal{T}(dims[1])
+end
+copy(M::Tridiagonal) = Tridiagonal(M.dl, M.d, M.du)
+
+# Operations on Tridiagonal matrices
+round(M::Tridiagonal) = Tridiagonal(round(M.dl), round(M.d), round(M.du))
+iround(M::Tridiagonal) = Tridiagonal(iround(M.dl), iround(M.d), iround(M.du))
+
+## Solvers
+# Allocation-free variants
+# Note that solve is non-aliasing, so you can use the same array for
+# input and output
+function solve(x::AbstractArray, xrng::Ranges{Int}, M::Tridiagonal, rhs::AbstractArray, rhsrng::Ranges{Int})
+ d = M.d
+ N = length(d)
+ if length(xrng) != N || length(rhsrng) != N
+ error("dimension mismatch")
+ end
+ dl = M.dl
+ du = M.du
+ dutmp = M.dutmp
+ rhstmp = M.rhstmp
+ xstart = first(xrng)
+ xstride = step(xrng)
+ rhsstart = first(rhsrng)
+ rhsstride = step(rhsrng)
+ # Forward sweep
+ denom = d[1]
+ dulast = du[1] / denom
+ dutmp[1] = dulast
+ rhslast = rhs[rhsstart] / denom
+ rhstmp[1] = rhslast
+ irhs = rhsstart+rhsstride
+ for i in 2:N-1
+ dltmp = dl[i-1]
+ denom = d[i] - dltmp*dulast
+ dulast = du[i] / denom
+ dutmp[i] = dulast
+ rhslast = (rhs[irhs] - dltmp*rhslast)/denom
+ rhstmp[i] = rhslast
+ irhs += rhsstride
+ end
+ dltmp = dl[N-1]
+ denom = d[N] - dltmp*dulast
+ xlast = (rhs[irhs] - dltmp*rhslast)/denom
+ # Backward sweep
+ ix = xstart + (N-2)*xstride
+ x[ix+xstride] = xlast
+ for i in N-1:-1:1
+ xlast = rhstmp[i] - dutmp[i]*xlast
+ x[ix] = xlast
+ ix -= xstride
+ end
+ return x
+end
+solve(x::StridedVector, M::Tridiagonal, rhs::StridedVector) = solve(x, 1:length(x), M, rhs, 1:length(rhs))
+function solve(M::Tridiagonal, rhs::StridedVector)
+ x = similar(rhs)
+ solve(x, M, rhs)
+end
+function solve(X::StridedMatrix, M::Tridiagonal, B::StridedMatrix)
+ if size(B, 1) != size(M, 1)
+ error("dimension mismatch")
+ end
+ if size(X) != size(B)
+ error("dimension mismatch in output")
+ end
+ m, n = size(B)
+ r = 1:m
+ for j = 1:n
+ r.start = (j-1)*m+1
+ solve(X, r, M, B, r)
+ end
+ return X
+end
+function solve(M::Tridiagonal, B::StridedMatrix)
+ X = similar(B)
+ solve(X, M, B)
+end
+
+# User-friendly solver
+\(M::Tridiagonal, rhs::Union(StridedVector,StridedMatrix)) = solve(M, rhs)
+
+# Tridiagonal multiplication
+function mult(x::AbstractArray, xrng::Ranges{Int}, M::Tridiagonal, v::AbstractArray, vrng::Ranges{Int})
+ dl = M.dl
+ d = M.d
+ du = M.du
+ N = length(d)
+ xi = first(xrng)
+ xstride = step(xrng)
+ vi = first(vrng)
+ vstride = step(vrng)
+ x[xi] = d[1]*v[vi] + du[1]*v[vi+vstride]
+ xi += xstride
+ for i = 2:N-1
+ x[xi] = dl[i-1]*v[vi] + d[i]*v[vi+vstride] + du[i]*v[vi+2*vstride]
+ xi += xstride
+ vi += vstride
+ end
+ x[xi] = dl[N-1]*v[vi] + d[N]*v[vi+vstride]
+ return x
+end
+mult(x::StridedVector, M::Tridiagonal, v::StridedVector) = mult(x, 1:length(x), M, v, 1:length(v))
+function mult(X::StridedMatrix, M::Tridiagonal, B::StridedMatrix)
+ if size(B, 1) != size(M, 1)
+ error("dimension mismatch")
+ end
+ if size(X) != size(B)
+ error("dimension mismatch in output")
+ end
+ m, n = size(B)
+ r = 1:m
+ for j = 1:n
+ r.start = (j-1)*m+1
+ solve(X, r, M, B, r)
+ end
+ return X
+end
+mult(X::StridedMatrix, M1::Tridiagonal, M2::Tridiagonal) = mult(X, M1, full(M2))
+function *(M::Tridiagonal, B::Union(StridedVector,StridedMatrix))
+ X = similar(B)
+ mult(X, M, B)
+end
+
+
+
+#### Woodbury matrices ####
+# This type provides support for the Woodbury matrix identity
+type Woodbury{T} <: AbstractMatrix{T}
+ A
+ U::Matrix{T}
+ C
+ Cp
+ V::Matrix{T}
+ tmpN1::Vector{T}
+ tmpN2::Vector{T}
+ tmpk1::Vector{T}
+ tmpk2::Vector{T}
+
+ function Woodbury(A::AbstractMatrix{T}, U::Matrix{T}, C, V::Matrix{T})
+ N = size(A, 1)
+ k = size(U, 2)
+ if size(A, 2) != N || size(U, 1) != N || size(V, 1) != k || size(V, 2) != N
+ error("Sizes do not match")
+ end
+ if k > 1
+ if size(C, 1) != k || size(C, 2) != k
+ error("Size of C is incorrect")
+ end
+ end
+ Cp = inv(inv(C) + V*(A\U))
+ # temporary space for allocation-free solver
+ tmpN1 = Array(T, N)
+ tmpN2 = Array(T, N)
+ tmpk1 = Array(T, k)
+ tmpk2 = Array(T, k)
+ # don't copy A, it could be huge
+ new(A, copy(U), copy(C), Cp, copy(V), tmpN1, tmpN2, tmpk1, tmpk2)
+ end
+end
+Woodbury{T}(A::AbstractMatrix{T}, U::Matrix{T}, C, V::Matrix{T}) = Woodbury{T}(A, U, C, V)
+Woodbury{T}(A::AbstractMatrix{T}, U::Vector{T}, C, V::Matrix{T}) = Woodbury{T}(A, reshape(U, length(U), 1), C, V)
+
+size(W::Woodbury) = size(W.A)
+function show(io, W::Woodbury)
+ println(io, summary(W), ":")
+ print(io, "A: ", W.A)
+ print(io, "\nU:\n")
+ print_matrix(io, W.U)
+ if isa(W.C, Matrix)
+ print(io, "\nC:\n")
+ print_matrix(io, W.C)
+ else
+ print(io, "\nC: ", W.C)
+ end
+ print(io, "\nV:\n")
+ print_matrix(io, W.V)
+end
+full{T}(W::Woodbury{T}) = convert(Matrix{T}, W)
+convert{T}(::Type{Matrix{T}}, W::Woodbury{T}) = full(W.A) + W.U*W.C*W.V
+function similar(W::Woodbury, T, dims::Dims)
+ if length(dims) != 2 || dims[1] != dims[2]
+ error("Woodbury matrices must be square")
+ end
+ n = size(W, 1)
+ k = size(W.U, 2)
+ return Woodbury{T}(similar(W.A), Array(T, n, k), Array(T, k, k), Array(T, k, n))
+end
+copy(W::Woodbury) = Woodbury(W.A, W.U, W.C, W.V)
+
+## Woodbury matrix routines ##
+
+function *(W::Woodbury, B::StridedVecOrMat)
+ return W.A*B + W.U*(W.C*(W.V*B))
+end
+function \(W::Woodbury, R::StridedVecOrMat)
+ AinvR = W.A\R
+ return AinvR - W.A\(W.U*(W.Cp*(W.V*AinvR)))
+end
+# Allocation-free solver for arbitrary strides (requires that W.A has a
+# non-aliasing "solve" routine, e.g., is Tridiagonal)
+function solve(x::AbstractArray, xrng::Ranges{Int}, W::Woodbury, rhs::AbstractArray, rhsrng::Ranges{Int})
+ solve(W.tmpN1, 1:length(W.tmpN1), W.A, rhs, rhsrng)
+ A_mul_B(W.tmpk1, W.V, W.tmpN1)
+ A_mul_B(W.tmpk2, W.Cp, W.tmpk1)
+ A_mul_B(W.tmpN2, W.U, W.tmpk2)
+ solve(W.tmpN2, W.A, W.tmpN2)
+ indx = first(xrng)
+ xinc = step(xrng)
+ for i = 1:length(W.tmpN2)
+ x[indx] = W.tmpN1[i] - W.tmpN2[i]
+ indx += xinc
+ end
+end
+solve(x::AbstractVector, W::Woodbury, rhs::AbstractVector) = solve(x, 1:length(x), W, rhs, 1:length(rhs))
+function solve(W::Woodbury, rhs::AbstractVector)
+ x = similar(rhs)
+ solve(x, W, rhs)
+end
+function solve(X::StridedMatrix, W::Woodbury, B::StridedMatrix)
+ if size(B, 1) != size(W, 1)
+ error("dimension mismatch")
+ end
+ if size(X) != size(B)
+ error("dimension mismatch in output")
+ end
+ m, n = size(B)
+ r = 1:m
+ for j = 1:n
+ r.start = (j-1)*m+1
+ solve(X, r, W, B, r)
+ end
+ return X
+end
+function solve(W::Woodbury, B::StridedMatrix)
+ X = similar(B)
+ solve(X, W, B)
+end
View
1  base/sysimg.jl
@@ -130,6 +130,7 @@ include("build_h.jl")
# linear algebra
include("linalg.jl")
include("linalg_dense.jl")
+include("linalg_specialized.jl")
include("linalg_blas.jl")
include("linalg_lapack.jl")
include("factorizations.jl")
View
61 extras/linalg_sparse.jl
@@ -572,64 +572,3 @@ diagmm{Tv,Ti,T}(A::SparseMatrixCSC{Tv,Ti}, b::Vector{T}) =
diagmm{T,Tv,Ti}(b::Vector{T}, A::SparseMatrixCSC{Tv,Ti}) =
diagmm!(SparseMatrixCSC(size(A,1),size(A,2),Ti[],Ti[],promote_type(Tv,T)[]), b, A)
-
-
-# Tridiagonal solver
-# Allocation-free variants
-function solve(x::Array, xstart::Int, xstride::Int, M::Tridiagonal, d::Array, dstart::Int, dstride::Int)
- # Grab refs to members (for efficiency)
- a = M.a
- b = M.b
- c = M.c
- cp = M.cp
- dp = M.dp
- N = length(b)
- # Forward sweep
- cp[1] = c[1] / b[1]
- dp[1] = d[dstart] / b[1]
- id = dstart+dstride
- for i = 2:N
- atmp = a[i]
- temp = b[i] - atmp*cp[i-1]
- cp[i] = c[i] / temp
- dp[i] = (d[id] - atmp*dp[i-1])/temp
- id += dstride
- end
- # Backward sweep
- ix = xstart + (N-2)*xstride
- x[ix+xstride] = dp[N]
- for i = N-1:-1:1
- x[ix] = dp[i] - cp[i]*x[ix+xstride]
- ix -= xstride
- end
-end
-function solve(x::Vector, M::Tridiagonal, d::Vector)
- if length(d) != length(M.b)
- error("Size mismatch between matrix and rhs")
- end
- solve(x, 1, 1, M, d, 1, 1)
-end
-# User-friendly solver
-function \(M::Tridiagonal, d::Vector)
- x = similar(d)
- solve(x, M, d)
- return x
-end
-
-# Tridiagonal multiplication
-function mult(x::Vector, M::Tridiagonal, v::Vector)
- a = M.a
- b = M.b
- c = M.c
- N = length(b)
- x[1] = b[1]*v[1] + c[1]*v[2]
- for i = 2:N-1
- x[i] = a[i]*v[i-1] + b[i]*v[i] + c[i]*v[i+1]
- end
- x[N] = a[N]*v[N-1] + b[N]*v[n]
-end
-function *(M::Tridiagonal, v::Vector)
- x = similar(v)
- mult(x, M, v)
- return x
-end
View
45 extras/sparse.jl
@@ -1023,48 +1023,3 @@ function assign(S::SparseAccumulator, v, i::Integer)
end
return S
end
-
-type Tridiagonal{T<:FloatingPoint} <: AbstractMatrix{T}
- a::Vector{T} # sub-diagonal
- b::Vector{T} # diagonal
- c::Vector{T} # sup-diagonal
- cp::Vector{T} # scratch space, sup-diagonal
- dp::Vector{T} # scratch space, rhs
-
- function Tridiagonal(N::Int)
- cp = Array(T, N)
- dp = Array(T, N)
- new(cp, cp, cp, cp, dp) # first three will be overwritten
- end
-end
-function Tridiagonal{T}(a::Vector{T}, b::Vector{T}, c::Vector{T})
- N = length(b)
- if length(a) != N || length(c) != N
- error("All three vectors must have the same length")
- end
- M = Tridiagonal{T}(N)
- M.a = copy(a)
- M.b = copy(b)
- M.c = copy(c)
- return M
-end
-size(M::Tridiagonal) = (length(M.b), length(M.b))
-function show(io, M::Tridiagonal)
- println(io, summary(M), ":")
- print_matrix(io, vcat((M.a)', (M.b)', (M.c)'))
-# println(io, " sub: ", (M.a)')
-# println(io, "diag: ", (M.b)')
-# println(io, " sup: ", (M.c)')
-end
-full{T}(M::Tridiagonal{T}) = convert(Matrix{T}, M)
-function convert{T}(::Type{Matrix{T}}, M::Tridiagonal{T})
- A = zeros(T, size(M))
- for i = 1:length(M.b)
- A[i,i] = M.b[i]
- end
- for i = 1:length(M.b)-1
- A[i+1,i] = M.a[i+1]
- A[i,i+1] = M.c[i]
- end
- return A
-end
View
56 test/lapack.jl
@@ -140,4 +140,60 @@ eA3 = [-1.50964415879218 -5.6325707998812 -4.934938326092;
0.367879439109187 1.47151775849686 1.10363831732856;
0.135335281175235 0.406005843524598 0.541341126763207]'
@assert norm((expm(A3) - eA3) ./ eA3) < 50000*eps()
+
+
+# basic tridiagonal operations
+n = 5
+d = 1 + rand(n)
+dl = -rand(n-1)
+du = -rand(n-1)
+T = Tridiagonal(dl, d, du)
+@assert size(T, 1) == n
+@assert size(T) == (n, n)
+F = diagm(d)
+for i = 1:n-1
+ F[i,i+1] = du[i]
+ F[i+1,i] = dl[i]
+end
+@assert full(T) == F
+
+# tridiagonal linear algebra
+Eps = sqrt(eps())
+v = randn(n)
+@assert T*v == F*v
+invFv = F\v
+@assert norm(T\v - invFv) < Eps
+@assert norm(solve(T,v) - invFv) < Eps
+B = randn(n,2)
+@assert norm(solve(T, B) - F\B) < Eps
+Tlu = LU(T)
+x = Tlu\v
+@assert norm(x - invFv) < Eps
+
+# symmetric tridiagonal
+Ts = Tridiagonal(dl, d, dl)
+Fs = full(Ts)
+invFsv = Fs\v
+Tldlt = LDLT(Ts)
+x = Tldlt\v
+@assert norm(x - invFsv) < Eps
+
+# eigenvalues/eigenvectors of symmetric tridiagonal
+DT, VT = eig(Ts)
+D, V = eig(Fs)
+@assert norm(DT - D) < Eps
+@assert norm(VT - V) < Eps
+
+
+# Woodbury
+U = randn(n,2)
+V = randn(2,n)
+C = randn(2,2)
+W = Woodbury(T, U, C, V)
+F = full(W)
+
+@assert norm(W*v - F*v) < Eps
+@assert norm(W\v - F\v) < Eps
+
+
end
View
1  test/sparse.jl
@@ -41,5 +41,4 @@ for i = 1 : 10
a = sprand(5, 4, 0.5)
@assert all([a[1:2,1:2] a[1:2,3:4]; a[3:5,1] [a[3:4,2:4]; a[5,2:4]]] == a)
end
-
end # cd
Please sign in to comment.
Something went wrong with that request. Please try again.