Rename chbv -> chbmv; improved in-place treatment

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

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

src/chbv.jl → src/chbmv.jl

@@ -1,4 +1,4 @@
-export chbv, chbv!
+export chbmv, chbmv!
 
 const α0 = 0.183216998528140087E-11
 
@@ -23,74 +23,84 @@ 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(θ)
 
 """
-    chbv(A, vec)
+    chbmv(A, vec[, fac])
 
 Calculate matrix exponential acting on some vector using the Chebyshev method.
 
 # Input
 
 - `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
 
-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.
 
     Roger B. Sidje ([email protected])
     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}, 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 chbv!{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 # chbv!
+end # chbmv!

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