Open GiggleLiu opened 6 years ago
@GiggleLiu I added support for LinearMaps.jl in this pull request of the currently active ExpmV.jl repository.
Cool! This interface looks good and the code is qualified. I will definitely have a try!
I'd like to post some benchmark under this issue latter, by comparing its performance in my application with Expokit's realization.
Thanks for your efforts!
First, your implementation is correct!
But somehow less efficient than Expokit
BenchmarkTools.Trial:
memory estimate: 2.24 MiB
allocs estimate: 134
--------------
minimum time: 190.535 ms (0.00% GC)
median time: 194.090 ms (0.00% GC)
mean time: 199.734 ms (1.57% GC)
maximum time: 254.911 ms (24.27% GC)
--------------
samples: 26
evals/sample: 1
BenchmarkTools.Trial:
memory estimate: 31.86 MiB
allocs estimate: 258890
--------------
minimum time: 12.443 s (0.16% GC)
median time: 12.443 s (0.16% GC)
mean time: 12.443 s (0.16% GC)
maximum time: 12.443 s (0.16% GC)
--------------
samples: 1
evals/sample: 1
BenchmarkTools.Trial:
memory estimate: 176.94 MiB
allocs estimate: 3677
--------------
minimum time: 2.393 s (4.12% GC)
median time: 2.464 s (7.00% GC)
mean time: 2.478 s (7.61% GC)
maximum time: 2.577 s (11.43% GC)
--------------
samples: 3
evals/sample: 1
BenchmarkTools.Trial:
memory estimate: 209.85 MiB
allocs estimate: 93982
--------------
minimum time: 3.602 s (4.90% GC)
median time: 3.651 s (6.74% GC)
mean time: 3.651 s (6.74% GC)
maximum time: 3.700 s (8.53% GC)
--------------
samples: 2
evals/sample: 1
The allocation is very high, would you please check the type stability?
The following profile information may be helpful to you.
17171 ./task.jl:85; (::getfield(Atom, Symbol("##108#113")){Dict{String,Any}})()
17171 /home/leo/.julia/packages/Atom/Pab0Z/src/eval.jl:89; macro expansion
17171 /home/leo/.julia/packages/Atom/Pab0Z/src/repl.jl:76; hideprompt(::getfield(Atom, Symbol("##109#114")){String,Int64,String})
17171 /home/leo/.julia/packages/Atom/Pab0Z/src/eval.jl:90; #109
17171 /home/leo/.julia/packages/Atom/Pab0Z/src/eval.jl:46; withpath
17171 /home/leo/.julia/packages/CodeTools/8CjY/src/utils.jl:30; withpath(::getfield(Atom, Symbol("##110#115")){String,Int64,String}, ::String)
17171 /home/leo/.julia/packages/Atom/Pab0Z/src/eval.jl:91; (::getfield(Atom, Symbol("##110#115")){String,Int64,String})()
17171 /home/leo/.julia/packages/CodeTools/8CjY/src/eval.jl:30; include_string(::Module, ::String, ::String, ::Int64)
17171 ./loading.jl:1002; include_string(::Module, ::String, ::String)
17171 ...x64/build/usr/share/julia/stdlib/v1.0/Profile/src/Profile.jl:25; top-level scope
2535 /home/leo/.julia/dev/Yao/src/Blocks/TimeEvolution.jl:23; apply!(::DefaultRegister{1,Complex{Float64}}, ::TimeEvolution{10,Float64,ChainBlock{10,Comp...
2535 ./reducedim.jl:305; mat
2463 ./reduce.jl:314; _mapreduce(::getfield(Yao.Blocks, Symbol("##41#42")), ::typeof(Base.mul_prod), ::IndexLine...
2334 ./reduce.jl:31; mul_prod(::SparseMatrixCSC{Complex{Float64},Int64}, ::SparseMatrixCSC{Complex{Float64},Int...
2334 ...ild/usr/share/julia/stdlib/v1.0/SparseArrays/src/linalg.jl:141; *
2334 ...ild/usr/share/julia/stdlib/v1.0/SparseArrays/src/linalg.jl:151; spmatmul
823 ...ild/usr/share/julia/stdlib/v1.0/SparseArrays/src/linalg.jl:178; #spmatmul#60(::Symbol, ::Function, ::SparseMatrixCSC{Complex{Float64},Int64}, ::SparseMa...
823 ./complex.jl:268; *
562 ./float.jl:399; *
943 ...ild/usr/share/julia/stdlib/v1.0/SparseArrays/src/linalg.jl:186; #spmatmul#60(::Symbol, ::Function, ::SparseMatrixCSC{Complex{Float64},Int64}, ::SparseMa...
679 ./complex.jl:266; +
679 ./float.jl:395; +
14433 /home/leo/jcode/lab/TE_benchmark.jl:21; apply2!(::DefaultRegister{1,Complex{Float64}}, ::TimeEvolution{12,Float64,PutBlock{12,1,XGa...
14433 /home/leo/.julia/packages/ExpmV/GFsI0/src/expmv_fun.jl:26; expmv
14123 /home/leo/.julia/packages/ExpmV/GFsI0/src/expmv_fun.jl:30; #expmv#1(::Nothing, ::String, ::Bool, ::Bool, ::Function, ::Complex{Float64}, ::LuxurySpar...
14123 ...leo/.julia/packages/ExpmV/GFsI0/src/select_taylor_degree.jl:63; macro expansion
13297 ...eo/.julia/packages/ExpmV/GFsI0/src/select_taylor_degree.jl:98; #select_taylor_degree#3(::Int64, ::Int64, ::String, ::Bool, ::Bool, ::Function, ::Luxury...
13296 /home/leo/.julia/packages/ExpmV/GFsI0/src/normAm.jl:30; normAm
8389 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:112; norm1est(::Int64, ::LuxurySparse.PermMatrix{Complex{Float64},Int64,Array{Complex{Float64...
1525 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:33; A_pow_n_B!(::Array{Complex{Float64},2}, ::LuxurySparse.PermMatrix{Complex{Float64},Int64...
1525 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:172; mul!
1525 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:578; generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMat...
680 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:638; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMa...
679 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:480; copy_transpose!
620 .../share/julia/stdlib/v1.0/LinearAlgebra/src/transpose.jl:197; copy_transpose!(::Array{Complex{Float64},2}, ::UnitRange{Int64}, ::UnitRange{Int64}, ...
732 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:646; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMa...
6862 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:35; A_pow_n_B!(::Array{Complex{Float64},2}, ::LuxurySparse.PermMatrix{Complex{Float64},Int64...
6862 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:172; mul!
6862 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:578; generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMat...
3005 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:638; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermM...
2997 ...usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:480; copy_transpose!
2762 .../share/julia/stdlib/v1.0/LinearAlgebra/src/transpose.jl:197; copy_transpose!(::Array{Complex{Float64},2}, ::UnitRange{Int64}, ::UnitRange{Int64},...
780 ./array.jl:771; setindex!
880 ./promotion.jl:425; getindex
1102 /home/leo/.julia/dev/LuxurySparse/src/PermMatrix.jl:46; getindex
882 ./array.jl:731; getindex
3364 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:646; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermM...
1574 ./complex.jl:268; *
1322 ./float.jl:399; *
1530 ./complex.jl:266; +
1530 ./float.jl:395; +
4878 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:182; norm1est(::Int64, ::LuxurySparse.PermMatrix{Complex{Float64},Int64,Array{Complex{Float64...
955 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:42; At_pow_n_B!(::Array{Complex{Float64},2}, ::LuxurySparse.PermMatrix{Complex{Float64},Int6...
954 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:172; mul!
954 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:578; generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMatr...
3923 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:44; At_pow_n_B!(::Array{Complex{Float64},2}, ::LuxurySparse.PermMatrix{Complex{Float64},Int6...
3922 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:172; mul!
3922 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:578; generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMat...
1705 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:638; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermM...
1698 ...usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:480; copy_transpose!
1556 .../share/julia/stdlib/v1.0/LinearAlgebra/src/transpose.jl:197; copy_transpose!(::Array{Complex{Float64},2}, ::UnitRange{Int64}, ::UnitRange{Int64},...
618 /home/leo/.julia/dev/LuxurySparse/src/PermMatrix.jl:46; getindex
1934 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:646; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermM...
962 ./complex.jl:268; *
792 ./float.jl:399; *
786 ./complex.jl:266; +
786 ./float.jl:395; +
Thanks a lot for your thorough testing. There is definitely some issue with allocations and I will take a deeper look at it.
Anyways, from the profile info that you posted it is evident how almost all the time is spent in the mul!
function, which is internal. I don't think that much can be done by reducing allocations.
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 https://github.com/acroy/Expokit.jl/pull/30