marcusps / expmv.jl Goto Github PK

The code currently computes the 1 norm using norm(A^m,1) instead of estimating more efficiently. This ends up being the bottleneck in the calculation of not so sparse matrices.

The MATLAB code uses normest1(), a proprietary function. However, the algorithm for normest1 is described in

Nicholas J. Higham and Françoise Tisseur (2000) A block algorithm for matrix 1-norm estimation, with an application to 1-norm pseudospectra. SIAM Journal On Matrix Analysis And Applications, 21 (4). pp. 1185-1201. ISSN 1095-7162 (preprint)

and it should not be too hard to reimplemented it, thereby avoiding this performance bottle neck.

Any chance of an update to be compatible with Julia 1.x? I'd be willing to help with this; #10 already updates to Julia 0.6 if that can be merged, and then one would update to Julia 0.7.

For larger random matrices I see the error given in the subject. For example,

const n = 500;
A = sprand(n,n, 0.4);
v = eye(n,1)[:];
ExpmV.expmv(1.,A,v)

gives more or less reliably

ERROR: DimensionMismatch("")
 in * at linalg/sparse.jl:112
 in normAm at /Users/acr/Development/Julia/ExpmV.jl/src/normAm.jl:56
 in select_taylor_degree at /Users/acr/Development/Julia/ExpmV.jl/src/select_taylor_degree.jl:66
 in expmv at /Users/acr/Development/Julia/ExpmV.jl/src/expmv_fun.jl:51

Just noticed that you use global variables in the benchmark. Currently, Julia isn't very good optimizing code which involves globals (basically since their type can change any time). Maybe you can wrap everything in a function?

Hey, would you consider finishing up the last pieces and registering? What exactly needs to be done for that? Where would you need help? I am considering using this as part of DifferentialEquations.jl to build some exponential integrators.

My problem of interest is imaginary time evolution of a Schrodinger-like equation, so I need to discretize time and re-normalize the wavefunction after each time step in order to avoid the trivial zero-solution.

In point 5 of their preprint, Mohy and Higham talk about efficiently computing e^{tA}B over a time interval. I think this can be useful in my case. Is it implemented in ExpmV?

LinearMaps.jl can map a matrix vector multiplication to a linear operator.

Is it possible to support expmv fallback for a general linear operator? I need it very much in my current project.

see this pr acroy/Expokit.jl#30

It looks like may be using norm(matrix). In Julia 0.7, this will compute the Frobenius norm (vecnorm in Julia 0.6), due to JuliaLang/julia#27401. If you want the induced/operator norm as in Julia 0.6, use opnorm(matrix) instead, or Compat.opnorm(matrix) to work in 0.6 and 0.7 (JuliaLang/Compat.jl#577).

Note that, for testing purposes, rather than @test norm(A - B) ≤ tol, it is usually preferred to do @test A ≈ B or @test A ≈ B rtol=... (which uses isapprox).