JuliaLinearAlgebra / AlgebraicMultigrid.jl

Algebraic Multigrid in Julia
Other
117 stars 23 forks source link

v0.7 upgrade, comments, bugfixes and cleanup #43

Closed mohamed82008 closed 6 years ago

mohamed82008 commented 6 years ago

This PR:

mohamed82008 commented 6 years ago

Not all deps are fully upgraded to v0.7 so they are failing the v0.7 and nightly tests. For instance, I had to use the JLD PR https://github.com/JuliaIO/JLD.jl/pull/215 to make sure tests pass on my machine. Not sure what the best fix for this is, maybe allow failures and wait?

codecov-io commented 6 years ago

Codecov Report

Merging #43 into master will decrease coverage by 3.27%. The diff coverage is 81.71%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master      #43      +/-   ##
==========================================
- Coverage   90.85%   87.58%   -3.28%     
==========================================
  Files          10       10              
  Lines         536      604      +68     
==========================================
+ Hits          487      529      +42     
- Misses         49       75      +26
Impacted Files Coverage Δ
src/splitting.jl 100% <100%> (ø) :arrow_up:
src/strength.jl 98.48% <100%> (-0.03%) :arrow_down:
src/aggregate.jl 90.76% <100%> (+1.53%) :arrow_up:
src/smoother.jl 69.01% <47.05%> (-5.32%) :arrow_down:
src/preconditioner.jl 56.25% <53.33%> (-10.42%) :arrow_down:
src/classical.jl 87.96% <80.3%> (-12.04%) :arrow_down:
src/utils.jl 94.82% <81.81%> (-1.33%) :arrow_down:
src/aggregation.jl 91.66% <86.66%> (-2.57%) :arrow_down:
src/multilevel.jl 78.04% <90.56%> (+5.32%) :arrow_up:

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 061a541...94bd0df. Read the comment docs.

ranjanan commented 6 years ago

Would it help to convert all the JLD files to use JLD2 instead?

mohamed82008 commented 6 years ago

I think JLD2 is itself not upgraded yet, but moving to JLD2 seems like the right thing anyways since JLD2 seems more supported than JLD.

mohamed82008 commented 6 years ago

JLD2 also only gives warnings but tests still pass, so we are good now. The only problem now is METADATA. The url there links to this package even though it uses the old AMG.jl. Not sure how to register the new name properly. If they were 2 different repos, this advice https://github.com/JuliaLang/METADATA.jl/issues/4951#issuecomment-208393242 may have helped. Any ideas?

mohamed82008 commented 6 years ago

v0.7 tests fail on Windows with JLD2, same error as in https://github.com/simonster/JLD2.jl/pull/79#issuecomment-400251115

ranjanan commented 6 years ago

Hey, sorry about the late reply: JuliaCon and all. :smile:

Okay so JLD and JLD2 both fail on Windows master, but the package is not directly dependent on them. If you can find any alternative method of saving this data, I'd be open to it.

mohamed82008 commented 6 years ago

I changed the files to simple .jl files that are included in tests. I did this by going to a string representation of the fields of the sparse arrays and then writing the .jl files that construct the sparse matrix.

Two more features to add before this PR is ready:

  1. Fix Jacobi smoother and make parallel
  2. Block GS and Jacobi smoothing
mohamed82008 commented 6 years ago

I implemented a type stable parallel smoothing using KissThreading. This package is not registered yet so we cannot officially depend on it, so I commented out the code for now. I will do one last test of this package then it should be ready to review and hopefully merge.

mohamed82008 commented 6 years ago

I believe this PR is ready for a review @ranjanan.

ranjanan commented 6 years ago

Wow, this has turned out to be quite a great pull request! Thanks! 😄

I am travelling at the moment so please give me a couple days to review this.

ranjanan commented 6 years ago

@mohamed82008 Unfortunately, there seems to be a performance regression, on Julia v0.6 Master:

Time for constructing AMG - Julia
BenchmarkTools.Trial: 
  memory estimate:  1.14 GiB
  allocs estimate:  13412189
  --------------
  minimum time:     2.234 s (34.17% GC)
  median time:      2.496 s (37.90% GC)
  mean time:        2.527 s (37.17% GC)
  maximum time:     2.851 s (38.88% GC)
  --------------
  samples:          3
  evals/sample:     1
Time for running CG - Julia
BenchmarkTools.Trial: 
  memory estimate:  1.25 GiB
  allocs estimate:  2428
  --------------
  minimum time:     4.473 s (2.39% GC)
  median time:      4.874 s (3.46% GC)
  mean time:        4.874 s (3.46% GC)
  maximum time:     5.274 s (4.37% GC)
  --------------
  samples:          2
  evals/sample:     1
sum(abs2, a * c2 - b) = 1.7685279192914618e-12

Your branch:

Time for constructing AMG - Julia
BenchmarkTools.Trial: 
  memory estimate:  3.37 GiB
  allocs estimate:  157650878
  --------------
  minimum time:     3.859 s (27.89% GC)
  median time:      3.889 s (29.20% GC)
  mean time:        3.889 s (29.20% GC)
  maximum time:     3.918 s (30.48% GC)
  --------------
  samples:          2
  evals/sample:     1
Time for running CG - Julia
BenchmarkTools.Trial: 
  memory estimate:  29.05 GiB
  allocs estimate:  1947220855
  --------------
  minimum time:     24.088 s (8.67% GC)
  median time:      24.088 s (8.67% GC)
  mean time:        24.088 s (8.67% GC)
  maximum time:     24.088 s (8.67% GC)
  --------------
  samples:          1
  evals/sample:     1
sum(abs2, a * c2 - b) = 1.7685279192914618e-12
mohamed82008 commented 6 years ago

I will look into it.

mohamed82008 commented 6 years ago

Can you post the exact code you used?

ranjanan commented 6 years ago

1m.jld.zip

using JLD, IterativeSolvers, AlgebraicMultigrid 
using BenchmarkTools
const AMG = AlgebraicMultigrid

function bench()
    a = load("/Users/ranjan/graphs/1m.jld")["G"]
    @show summary(a)
    b = zeros(size(a,1))
    b[98531] = -1
    b[332615] = 1
    ml2 = smoothed_aggregation(a)
    t2 = @benchmark smoothed_aggregation($a)
    println("Time for constructing AMG - Julia")
    display(t2)
    p2 = AMG.aspreconditioner(ml2)
    c2 = cg(a, b, Pl = p2, maxiter = 100_000, tol = 1e-6, verbose = true) 
    t2 = @benchmark x = cg($a, $b, Pl = $p2, maxiter = 100_000, tol = 1e-6) 
    println("Time for running CG - Julia")
    display(t2) 
    @show sum(abs2, a*c2 - b)
    ml
end

bench()
mohamed82008 commented 6 years ago

The performance regression in the operator building comes from making the code identical to the algorithm for unsymmetric matrices. This leads to unnecessary copying and transposing when the matrix is symmetric. Some of these can be avoided for both symmetric and unsymmetric matrices, but others are necessary for unsymmetric matrices so my solution to this is to dispatch on Symmetric to avoid unnecessary copying and transposing when the matrix is symmetric.

mohamed82008 commented 6 years ago

In other words, some extra allocations are not avoidable without being wrong for unsymmetric matrices, or at least not following the known algorithm. One can argue that the direction of being strongly coupled to some other variable is arbitrary and so we can freely swap the semantics of rows and columns, but I am not sure we should do that. I will make a correction that should make the operator building allocate less when wrapping the matrix with the Symmetric wrapper.

mohamed82008 commented 6 years ago

Oh also in v0.6 doing copy(A') allocates twice so I need to correct for that.

mohamed82008 commented 6 years ago

Interestingly the real culprit was failed inference in https://github.com/mohamed82008/AlgebraicMultigrid.jl/blob/76ac8f15e8314174bfdba559095d12a23ea5e2fe/src/smoother.jl#L29. Changing 0 to z solved the problem but I don't understand why this is a problem in the first place. Could it be a Julia bug?

mohamed82008 commented 6 years ago

Oh I see, ifelse needs to return a single type, missed that.

mohamed82008 commented 6 years ago
summary(a) = "868150×868150 SparseMatrixCSC{Float64,Int64}"
Time for constructing AMG - Julia
BenchmarkTools.Trial:
  memory estimate:  1.15 GiB
  allocs estimate:  13412264
  --------------
  minimum time:     1.484 s (34.11% GC)
  median time:      1.538 s (36.10% GC)
  mean time:        1.536 s (35.58% GC)
  maximum time:     1.584 s (35.96% GC)
  --------------
  samples:          4
  evals/sample:     1

Time for running CG - Julia
BenchmarkTools.Trial:
  memory estimate:  33.13 MiB
  allocs estimate:  295
  --------------
  minimum time:     3.011 s (0.01% GC)
  median time:      3.019 s (0.04% GC)
  mean time:        3.019 s (0.04% GC)
  maximum time:     3.027 s (0.07% GC)
  --------------
  samples:          2
  evals/sample:     1
sum(abs2, a * c2 - b) = 1.7681244199514989e-12
1.7681244199514989e-12
ranjanan commented 6 years ago

Sorry I was a little lax on this, I was moving cities. Thanks again for the great PR!