From cb9457afa4889757489b7f3e27622c8edde184a6 Mon Sep 17 00:00:00 2001 From: Alexander Croy Date: Fri, 24 Nov 2017 14:33:53 +0100 Subject: [PATCH 1/3] Rename chbv to chbmv for consistency. --- README.md | 2 +- src/Expokit.jl | 2 +- src/{chbv.jl => chbmv.jl} | 24 ++++++++++++------------ test/runtests.jl | 28 ++++++++++++++-------------- 4 files changed, 28 insertions(+), 28 deletions(-) rename src/{chbv.jl => chbmv.jl} (88%) diff --git a/README.md b/README.md index 01fbb8d..04cd120 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ for large sparse matrices. For more details about the methods see *R.B. Sidje, ACM Trans. Math. Softw., 24(1):130-156, 1998* (or [its preprint](http://www.maths.uq.edu.au/expokit/paper.pdf)). -**Note:** Apart from `expmv` (which is called `expv` in EXPOKIT) also `phimv`, `padm` and `chbv` are available. +**Note:** Apart from `expmv` (which is called `expv` in EXPOKIT) also `phimv`, `padm` and `chbmv` are available. ## Usage ```julia diff --git a/src/Expokit.jl b/src/Expokit.jl index 1f67b92..420bcdf 100644 --- a/src/Expokit.jl +++ b/src/Expokit.jl @@ -19,7 +19,7 @@ const expm! = Base.LinAlg.expm! include("expmv.jl") include("padm.jl") -include("chbv.jl") +include("chbmv.jl") include("phimv.jl") end # module diff --git a/src/chbv.jl b/src/chbmv.jl similarity index 88% rename from src/chbv.jl rename to src/chbmv.jl index f9188e8..e3aa69c 100644 --- a/src/chbv.jl +++ b/src/chbmv.jl @@ -1,4 +1,4 @@ -export chbv, chbv! +export chbmv, chbmv! const α0 = 0.183216998528140087E-11 @@ -40,7 +40,7 @@ const θconj = [0.562314417475317895E+01 + im*0.119406921611247440E+01, -0.889777151877331107E+01 + im*0.166309842834712071E+02] """ - chbv(A, vec) + chbmv(A, vec) Calculate matrix exponential acting on some vector using the Chebyshev method. @@ -51,16 +51,16 @@ Calculate matrix exponential acting on some vector using the Chebyshev method. # Notes -This Julia implementation is based on Expokit's CHBV Matlab code by +This Julia implementation is based on Expokit's chbv Matlab code by Roger B. Sidje, see below. ---- +--- y = chbv( H, x ) - CHBV computes the direct action of the matrix exponential on + chbv computes the direct action of the matrix exponential on a vector: y = exp(H) * x. It uses the partial fraction expansion of the uniform rational Chebyshev approximation of type (14,14). - About 14-digit accuracy is expected if the matrix H is symmetric + About 14-digit accuracy is expected if the matrix H is symmetric negative definite. The algorithm may behave poorly otherwise. See also PADM, EXPOKIT. @@ -68,14 +68,14 @@ Roger B. Sidje, see below. EXPOKIT: Software Package for Computing Matrix Exponentials. ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 """ -function chbv{T}(A, vec::Vector{T}) +function chbmv{T}(A, vec::Vector{T}) result = convert(Vector{promote_type(eltype(A), T)}, copy(vec)) - return chbv!(result, A, vec) + return chbmv!(result, A, vec) end -chbv!{T}(A, vec::Vector{T}) = chbv!(vec, A, copy(vec)) +chbmv!{T}(A, vec::Vector{T}) = chbmv!(vec, A, copy(vec)) -function chbv!{T<:Real}(w::Vector{T}, A, vec::Vector{T}) +function chbmv!{T<:Real}(w::Vector{T}, A, vec::Vector{T}) p = min(length(θ), length(α)) scale!(copy!(w, vec), α0) @inbounds for i = 1:p @@ -84,7 +84,7 @@ function chbv!{T<:Real}(w::Vector{T}, A, vec::Vector{T}) return w end -function chbv!{T<:Complex}(w::Vector{T}, A, vec::Vector{T}) +function chbmv!{T<:Complex}(w::Vector{T}, A, vec::Vector{T}) p = min(length(θ), length(α)) scale!(copy!(w, vec), α0) t = [θ; θconj] @@ -93,4 +93,4 @@ function chbv!{T<:Complex}(w::Vector{T}, A, vec::Vector{T}) w .+= (A - t[i]*I) \ (a[i] * vec) end return w -end # chbv! +end # chbmv! diff --git a/test/runtests.jl b/test/runtests.jl index 3b377f0..6e19fc8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -106,17 +106,17 @@ res, t1, t2 = test_padm2(1000) println("residuum: $res\n") @test res < 1e-6 -function test_chbv(n::Int) +function test_chbmv(n::Int) p = 0.1 D = diagm(-rand(n)) T = sprandn(n, n, p) H = T * D * T.' # random negative semidefinite symmetric matrix vec = randn(n) - w1 = chbv(H, vec); + w1 = chbmv(H, vec); tic() - w1 = chbv(H, vec) + w1 = chbmv(H, vec) t1 = toc() w2 = expm(full(H)) * vec @@ -127,17 +127,17 @@ function test_chbv(n::Int) return norm(w1-w2)/norm(w2), t1, t2 end -function test_chbv2(n::Int) +function test_chbmv2(n::Int) p = 0.1 D = diagm(-rand(n)) T = sprandn(n, n, p) + sprandn(n, n, p)*im H = T * D * T' # random negative semidefinite hermitian matrix vec = randn(n) + randn(n)*im - w1 = chbv(H, vec); + w1 = chbmv(H, vec); tic() - w1 = chbv(H, vec) + w1 = chbmv(H, vec) t1 = toc() w2 = expm(full(H)) * vec @@ -148,23 +148,23 @@ function test_chbv2(n::Int) return norm(w1-w2)/norm(w2), t1, t2 end -println("testing real n=100 (first chbv, then expm)") -res, t1, t2 = test_chbv(100) +println("testing real n=100 (first chbmv, then expm)") +res, t1, t2 = test_chbmv(100) println("residuum: $res\n") @test res < 1e-6 -println("testing complex n=100 (first chbv, then expm)") -res, t1, t2 = test_chbv2(100) +println("testing complex n=100 (first chbmv, then expm)") +res, t1, t2 = test_chbmv2(100) println("residuum: $res\n") @test res < 1e-6 -println("testing real n=1000 (first chbv, then expm)") -res, t1, t2 = test_chbv(1000) +println("testing real n=1000 (first chbmv, then expm)") +res, t1, t2 = test_chbmv(1000) println("residuum: $res\n") @test res < 1e-6 -println("testing complex n=1000 (first chbv, then expm)") -res, t1, t2 = test_chbv2(1000) +println("testing complex n=1000 (first chbmv, then expm)") +res, t1, t2 = test_chbmv2(1000) println("residuum: $res\n") @test res < 1e-6 From e02c24541fcd5c1efefef3ad4c5184636baea9c8 Mon Sep 17 00:00:00 2001 From: Alexander Croy Date: Fri, 24 Nov 2017 15:20:46 +0100 Subject: [PATCH 2/3] Use A_ldiv_B! and improve in-place behavior. --- src/chbmv.jl | 50 +++++++++++++++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 21 deletions(-) diff --git a/src/chbmv.jl b/src/chbmv.jl index e3aa69c..56021e4 100644 --- a/src/chbmv.jl +++ b/src/chbmv.jl @@ -23,21 +23,8 @@ const θ = [0.562314417475317895E+01 - im*0.119406921611247440E+01, -0.370327340957595652E+01 - im*0.136563731924991884E+02, -0.889777151877331107E+01 - im*0.166309842834712071E+02] -const αconj = [-0.557503973136501826E+02 - im*0.204295038779771857E+03, - 0.938666838877006739E+02 + im*0.912874896775456363E+02, - -0.469965415550370835E+02 - im*0.116167609985818103E+02, - 0.961424200626061065E+01 - im*0.264195613880262669E+01, - -0.752722063978321642E+00 + im*0.670367365566377770E+00, - 0.188781253158648576E-01 - im*0.343696176445802414E-01, - -0.143086431411801849E-03 + im*0.287221133228814096E-03] - -const θconj = [0.562314417475317895E+01 + im*0.119406921611247440E+01, - 0.508934679728216110E+01 + im*0.358882439228376881E+01, - 0.399337136365302569E+01 + im*0.600483209099604664E+01, - 0.226978543095856366E+01 + im*0.846173881758693369E+01, - -0.208756929753827868E+00 + im*0.109912615662209418E+02, - -0.370327340957595652E+01 + im*0.136563731924991884E+02, - -0.889777151877331107E+01 + im*0.166309842834712071E+02] +const αconj = conj(α) +const θconj = conj(θ) """ chbmv(A, vec) @@ -75,22 +62,43 @@ end chbmv!{T}(A, vec::Vector{T}) = chbmv!(vec, A, copy(vec)) -function chbmv!{T<:Real}(w::Vector{T}, A, vec::Vector{T}) +function chbmv!{T<:Real}(w::Vector{T}, A, vec::Vector{T}, fac=factorize) p = min(length(θ), length(α)) scale!(copy!(w, vec), α0) + tmp_v = similar(vec, Complex128) + B = A - θ[1]*I # make copy of A with diagonal elements + @inbounds for i = 1:p - w .+= real((A - θ[i]*I) \ (α[i] * vec)) + scale!(copy!(tmp_v, vec), α[i]) + copy!(B,A) + for n=1:size(B,1) + B[n,n] = A[n,n] - θ[i] + end + A_ldiv_B!( fac(B), tmp_v) + w .+= real.( tmp_v ) end + return w -end +end # chbmv -function chbmv!{T<:Complex}(w::Vector{T}, A, vec::Vector{T}) +function chbmv!{T<:Complex}(w::Vector{T}, A, vec::Vector{T}, fac=factorize) p = min(length(θ), length(α)) - scale!(copy!(w, vec), α0) t = [θ; θconj] a = 0.5 * [α; αconj] + + scale!(copy!(w, vec), α0) + tmp_v = similar(vec) + B = A - θ[1]*I # make copy of A with diagonal elements + @inbounds for i = 1:2*p - w .+= (A - t[i]*I) \ (a[i] * vec) + scale!(copy!(tmp_v, vec), a[i]) + copy!(B,A) + for n=1:size(B,1) + B[n,n] = A[n,n] - t[i] + end + A_ldiv_B!( fac(B), tmp_v) + w .+= tmp_v end + return w end # chbmv! From e3b5d017607d9f664fdc108fa6755c9c1ece083c Mon Sep 17 00:00:00 2001 From: Alexander Croy Date: Sun, 26 Nov 2017 17:59:45 +0100 Subject: [PATCH 3/3] Update docstring of chbmv to include optional argument fac. --- src/chbmv.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/chbmv.jl b/src/chbmv.jl index 56021e4..e856899 100644 --- a/src/chbmv.jl +++ b/src/chbmv.jl @@ -27,7 +27,7 @@ const αconj = conj(α) const θconj = conj(θ) """ - chbmv(A, vec) + chbmv(A, vec[, fac]) Calculate matrix exponential acting on some vector using the Chebyshev method. @@ -35,6 +35,8 @@ Calculate matrix exponential acting on some vector using the Chebyshev method. - `A` -- matrix which can be dense or sparse - `vec` -- vector on which the matrix exponential of `A` is applied +- `fac` -- optional factorization function (default: `factorize`); + `fac` may change its input to avoid allocation # Notes