The tag name "0.0.3" is not of the appropriate SemVer form (vX.Y.Z).
cc: @ranjanan

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Do we need compat as a dependency any more? It should be possible to drop it.

Also, why do we need such a specific version of IterativeSolvers.jl. Perhaps we want 0.8.3+ or some such thing. Isn't 0.8 sufficient unless versions prior to 0.8.3 are broken?

Would it be possible to get some code like below working? The first example with the function f is meant to show that this definition of hessian_vector_product can work. The second example shows that this fails with g, which uses AlgebraicMultigrid. If Hessian-vector products could be computed efficiently this way, it would be really useful.

using Test
import AlgebraicMultigrid
import ForwardDiff
import LinearAlgebra
import SparseArrays
import Zygote

hessian_vector_product(f, x, v) = ForwardDiff.jacobian(s->Zygote.gradient(f, x + s[1] * v)[1], [0.0])[:]

n = 4
A = randn(n, n)
hessian = A + A'
f(x) = LinearAlgebra.dot(x, A * x) 
x = randn(n)
v = randn(n)
hvp1 = hessian_vector_product(f, x, v)
hvp2 = hessian * v
@test hvp1 ≈ hvp2#the hessian_vector_product plausibly works!

function g(x)
	k = x[1:n + 1]
	B = SparseArrays.spdiagm(0=>k[1:end - 1] + k[2:end], -1=>-k[2:end - 1], 1=>-k[2:end - 1])
	ml = AlgebraicMultigrid.ruge_stuben(B)
	return sum(AlgebraicMultigrid.solve(ml, x[N + 2:end]))
end
x = randn(2 * n + 1)
v = randn(2 * n + 1)
hessian_vector_product(g, x, v)#seems to fail during the coarse solve in AlgebraicMultigrid

The "Coarse Solver" entry (printed to console when inspecting a MultiLevel solver) contains numerical values for pinv. This takes a lot of space and is hard to read (particularly in a Jupyter notebook).

The default example:

using AlgebraicMultigrid
A = poisson(1000)
ml = ruge_stuben(A)

shows:

Multilevel Solver
-----------------
Operator Complexity: 1.986
Grid Complexity: 1.986
No. of Levels: 8
Coarse Solver: AlgebraicMultigrid.Pinv{Float64}([93.98601398601403 80.55944055944069 67.1328671328673 53.7062937062939 40.279720279720436 26.853146853146935 13.42657342657344; 80.55944055944065 178.76523476523496 148.97102897102926 119.17682317682355 89.38261738261767 59.58841158841175 29.79420579420582; 67.13286713286729 148.9710289710293 230.80919080919125 184.64735264735313 138.48551448551484 92.32367632367648 46.1618381618382; 53.70629370629386 119.17682317682356 184.64735264735316 250.11788211788266 187.588411588412 125.0589410589412 62.529470529470565; 40.279720279720436 89.38261738261772 138.48551448551498 187.58841158841216 236.69130869130905 157.7942057942059 78.89710289710294; 26.853146853147013 59.58841158841189 92.32367632367674 125.05894105894147 157.79420579420605 190.5294705294706 95.2647352647353; 13.426573426573487 29.794205794205915 46.16183816183832 62.52947052947072 78.89710289710304 95.26473526473532 111.63236763236768])
Level     Unknowns     NonZeros
-----     --------     --------
    1         1000         2998 [50.35%]
    2          500         1498 [25.16%]
    3          250          748 [12.56%]
    4          125          373 [ 6.26%]
    5           62          184 [ 3.09%]
    6           31           91 [ 1.53%]
    7           15           43 [ 0.72%]
    8            7           19 [ 0.32%]

Related source code:

One quick fix is to print typeof(ml.coarse_solver), instead of ml.coarse_solver, so that the summary is easier to read:

Multilevel Solver
-----------------
Operator Complexity: 1.986
Grid Complexity: 1.986
No. of Levels: 8
Coarse Solver: AlgebraicMultigrid.Pinv{Float64}
Level     Unknowns     NonZeros
-----     --------     --------
    1         1000         2998 [50.35%]
    2          500         1498 [25.16%]
    3          250          748 [12.56%]
    4          125          373 [ 6.26%]
    5           62          184 [ 3.09%]
    6           31           91 [ 1.53%]
    7           15           43 [ 0.72%]
    8            7           19 [ 0.32%]

This is also consistent with PyAMG which prints Coarse Solver: 'pinv2' without numerical values.

Package version

[email protected]

LoadError: UndefVarError: SparseMatrixCSC not defined

It appears that the verbose argument in solve is ignored. I would suggest dropping the argument instead of implementing it, since the log of residuals is plenty.

For a test problem (Laplace equation on a 200^3 grid), ruge_stuben allocates ~8 GB of memory. For comparison, LinearSolve.KrylovJL_CG allocates ~90 MB.

Is this normal / can it be optimized?

Thanks!

As described in https://discourse.julialang.org/t/ann-plans-for-removing-packages-that-do-not-yet-support-1-0-from-the-general-registry/ we are planning on removing packages that do not support 1.0 from the General registry. This package has been detected to not support 1.0 and is thus slated to be removed. The removal of packages from the registry will happen approximately a month after this issue is open.

To transition to the new Pkg system using Project.toml, see https://github.com/JuliaRegistries/Registrator.jl#transitioning-from-require-to-projecttoml.
To then tag a new version of the package, see https://github.com/JuliaRegistries/Registrator.jl#via-the-github-app.

If you believe this package has erronously been detected as not supporting 1.0 or have any other questions, don't hesistate to can be discussed here or in the thread linked at the top of this post.

Hi - thank you for developing this, I'll be very happy to finally retire PyAMG.jl.

One question: the order of arguments currently is

solve(ml, b, maxiter, cycle, tol; verbose = false, log = false)

I think most users would really prefer to call this with

solve(ml, b; maxiter=..., cycle=..., tol=..., verbose = false, log = false)

but if you prefer to keep non-kw-args, then putting tol in front seems the more natural option.

I couldn't find any support from literature as to why the column scaling by maximum (positive) value in the column is done before taking the absolute value. Did I miss something or is this a bug? It would be the same if the largest entry by magnitude in each column was a positive one anyways. But if one negative entry has a high magnitude, then the resulting strength matrix will have some strength entries > 1. This is also inconsistent with the Classical strength calculation function that takes abs value before scaling. Swapping the 2 did not change my test results. If you approve this change, I will make a PR next.

Need a ton more tests for SA

Hello everyone,

I want to use the AlgebraicMultigrid in julia - Version 1.0.2 (2018-11-08) to solve the equations, and I run the example like that

julia> using AlgebraicMultigrid

julia> A = poisson(1000);

julia> ml = ruge_stuben(A)

but I get the following error

julia> ml = ruge_stuben(A)
Error showing value of type AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64},GaussSeidel{SymmetricSweep},GaussSeidel{SymmetricSweep},SparseArrays.SparseMatrixCSC{Float64,Int64},LinearAlgebra.Adjoint{Float64,SparseArrays.SparseMatrixCSC{Float64,Int64}},SparseArrays.SparseMatrixCSC{Float64,Int64},AlgebraicMultigrid.MultiLevelWorkspace{Array{Float64,1},1}}:
ERROR: MethodError: no method matching round(::Float64, ::Int64)
Closest candidates are:
  round(::Float64, ::RoundingMode{:Nearest}) at float.jl:368
  round(::Float64, ::RoundingMode{:Up}) at float.jl:366
  round(::Float64, ::RoundingMode{:Down}) at float.jl:364
  ...
Stacktrace:
 [1] show(::IOContext{REPL.Terminals.TTYTerminal}, ::AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64},GaussSeidel{SymmetricSweep},GaussSeidel{SymmetricSweep},SparseArrays.SparseMatrixCSC{Float64,Int64},LinearAlgebra.Adjoint{Float64,SparseArrays.SparseMatrixCSC{Float64,Int64}},SparseArrays.SparseMatrixCSC{Float64,Int64},AlgebraicMultigrid.MultiLevelWorkspace{Array{Float64,1},1}}) at ./printf.jl:159
...
...

Can this error be solved? Thank you!

I think it would be valuable to add more different matrices to the tests (e.g. non-symmetric). We can pull these via MatrixDepot.jl.

Makes smoothed aggregation unstable when using diagonal weighting on the Jacobi prolongation

Hello,

Sorry if this is not the right place to post questions.

I have recently installed julia 1.5.3 and AlgebraicMultigrid package. I have previously deleted earlier versions of julia (included hidden files), so the installation is clean. When I execute
ml = AlgebraicMultigrid.smoothed_aggregation(AUU, strength=SymmetricStrength(), max_levels=10, max_coarse=10)
I get the following error
ERROR: MethodError: no method matching smoothed_aggregation(::SparseMatrixCSC{Float64,Int64}; strength=SymmetricStrength{Float64}(0.0), max_levels=10, max_coarse=10) Closest candidates are: smoothed_aggregation(::TA) where {T, V, TA<:SparseMatrixCSC{T,V}} at /home/jcarpio/.julia/packages/AlgebraicMultigrid/RU7pA/src/aggregation.jl:1 got unsupported keyword arguments "strength", "max_levels", "max_coarse" smoothed_aggregation(::TA, ::Type{Val{bs}}) where {T, V, bs, TA<:SparseMatrixCSC{T,V}} at /home/jcarpio/.julia/packages/AlgebraicMultigrid/RU7pA/src/aggregation.jl:1 got unsupported keyword arguments "strength", "max_levels", "max_coarse" smoothed_aggregation(::TA, ::Type{Val{bs}}, ::Any) where {T, V, bs, TA<:SparseMatrixCSC{T,V}} at /home/jcarpio/.julia/packages/AlgebraicMultigrid/RU7pA/src/aggregation.jl:1 got unsupported keyword arguments "strength", "max_levels", "max_coarse"

Strangely, the code works perfectly in other machines with julia 1.5.3 (but which also have some earlier versions installed), but not in a AMD EPYC 7601 32-Core Processor.

Am I doing something wrong?

I am trying to migrate to AlgebraicMultigrid.jl from PyAMG.jl which is apparently being abandoned based on the discussion here. The first thing I want to do is get the tests working for my project which relies on a Ruge-Stuben preconditioner. The tests involve solving small linear systems, and I ran into a problem with AlgebraicMultigrid.jl for these systems. The same problem occurs with a slight modification of the example code from the README, and seems to arise when there is only 1 level in the MultiLevel object. Here's a minimal example that illustrates the issue:

using AlgebraicMultigrid
function f(n)
	A = poisson(n)
	ml = ruge_stuben(A)
	return solve(ml, A * ones(n))
end
for i = 10:-1:1
	print("starting $(2^i) ")
	f(2^i)
	println("done with $(2^i)")
end

and here is the output I see

starting 1024 done with 1024
starting 512 done with 512
starting 256 done with 256
starting 128 done with 128
starting 64 done with 64
starting 32 done with 32
starting 16 done with 16
starting 8 ERROR: LoadError: BoundsError: attempt to access 0-element Array{Array{Float64,1},1} at index [1]
Stacktrace:
 [1] getindex at ./array.jl:731 [inlined]
 [2] #solve!#12(::Int64, ::Float64, ::Bool, ::Bool, ::Bool, ::Function, ::Array{Float64,1}, ::AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64},GaussSeidel{SymmetricSweep},GaussSeidel{SymmetricSweep},SparseArrays.SparseMatrixCSC{Float64,Int64},LinearAlgebra.Adjoint{Float64,SparseArrays.SparseMatrixCSC{Float64,Int64}},SparseArrays.SparseMatrixCSC{Float64,Int64},AlgebraicMultigrid.MultiLevelWorkspace{Array{Float64,1},1}}, ::Array{Float64,1}, ::AlgebraicMultigrid.V) at /Users/omalled/.julia/packages/AlgebraicMultigrid/p0Evn/src/multilevel.jl:161
 [3] #solve#11 at /Users/omalled/.julia/packages/AlgebraicMultigrid/p0Evn/src/multilevel.jl:151 [inlined]
 [4] solve at /Users/omalled/.julia/packages/AlgebraicMultigrid/p0Evn/src/multilevel.jl:138 [inlined]
 [5] f(::Int64) at /Users/omalled/blah/amg/runtests.jl:5
 [6] top-level scope at /Users/omalled/blah/amg/runtests.jl:9 [inlined]
 [7] top-level scope at ./none:0
 [8] include at ./boot.jl:317 [inlined]
 [9] include_relative(::Module, ::String) at ./loading.jl:1038
 [10] include(::Module, ::String) at ./sysimg.jl:29
 [11] include(::String) at ./client.jl:398
 [12] top-level scope at none:0
in expression starting at /Users/omalled/blah/amg/runtests.jl:7

I am using julia v0.7.0 and AlgebraicMultigrid v0.2.0.

mat.jls.zip
vec.jls.zip

Ref Circuitscape/Circuitscape.jl#129

This NeurIPS 2021 paper seems interesting and relevant. https://openreview.net/pdf?id=xXYjxli-2i

When all nodes aren't aggregated, we need to construct a matrix excluding that node. That doesn't currently happen, which leads to a BoundsError.

Hi! First, let me congratulate you for this package, I have been using it and it performs very well! It's faster than PETSc GAMG in many cases according to my experience.

I am wondering if it is possible to pass a user-defined near null space to the SA-AMG routines. This is important e.g. in linear elasticity and other vector-valued problems.

By looking at this line

it looks like the near null space is hardcoded to be a vector of ones. Thus, I guess that the answer to my question is NO at the moment.

I would like to know if somebody is working on this feature. If not, I can try to address it myself and file a PR.

@ranjanan Can we upgrade this package to use the usual tooling like Circuitscape.jl and other packages?

The current implementation calls n iterations of n iterationsForwardSweep() and BackwardSweep() each. (n^2). Should only do n times.

First of all, I'd like to express my appreciation for the very nice preconditioners! They are almost magic, bringing down the number of iterations required to solve Poisson's problem from almost maxiter to about 5!

I discovered the following error when testing my package on AppVeyor, which also tests on Win32:

MethodError: no method matching AlgebraicMultigrid.Solver(::AlgebraicMultigrid.Classical{Float64}, ::AlgebraicMultigrid.RS, ::AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, ::AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, ::Int32, ::Int32)
Closest candidates are:
  AlgebraicMultigrid.Solver(::S, ::T, ::P, ::PS, !Matched::Int64, !Matched::Int64) where {S, T, P, PS} at C:\Users\appveyor\.julia\packages\AlgebraicMultigrid\xJj9G\src\classical.jl:2

The reason is that max_levels and max_coarse are hard-coded to be Int64:
https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl/blob/master/src/classical.jl#L6
Shouldn't Int be enough?

The following part of the code:

This is currently not the case which means that we are computing the wrong strength matrix S when A is not symmetric. Of course, in CSC matrices, it is more natural to do column-based iteration, so the simplest fix is to change this:

In the Python version, S gives row access and T = S' gives column access of S. In the Julia version, S's columns are Python's S's rows and T's columns are Python's S's columns. So we get row and column access of Python's S using only column access in Julia. This also makes the remaining code compatible, e.g. we are scaling S by columns which does the same thing as scaling it by rows in Python, etc.

Since this package overloads Base.length, without specifying ml::MultiLevel as the argument type, you get weird LoadErrors if you have this package loaded and you call length on a object which has no method defined. E.g. length(nothing) results in

ERROR: LoadError: type Nothing has no field levels

I am trying the example provided on the github page. When ml = ruge_stuben(A) is executed, it throws out the error:

Error showing value of type AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64},GaussSeidel{SymmetricSweep},GaussSeidel{SymmetricSweep},SparseArrays.SparseMatrixCSC{Float64,Int64},LinearAlgebra.Adjoint{Float64,SparseArrays.SparseMatrixCSC{Float64,Int64}},SparseArrays.SparseMatrixCSC{Float64,Int64},AlgebraicMultigrid.MultiLevelWorkspace{Array{Float64,1},1}}:
ERROR: MethodError: no method matching round(::Float64, ::Int64)
Closest candidates are:
round(::Float64, ::RoundingMode{:Nearest}) at float.jl:370
round(::Float64, ::RoundingMode{:Up}) at float.jl:368
round(::Float64, ::RoundingMode{:Down}) at float.jl:366
...
Stacktrace:
[1] show(::IOContext{REPL.Terminals.TTYTerminal}, ::AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64},GaussSeidel{SymmetricSweep},GaussSeidel{SymmetricSweep},SparseArrays.SparseMatrixCSC{Float64,Int64},LinearAlgebra.Adjoint{Float64,SparseArrays.SparseMatrixCSC{Float64,Int64}},SparseArrays.SparseMatrixCSC{Float64,Int64},AlgebraicMultigrid.MultiLevelWorkspace{Array{Float64,1},1}}) at /Users/huizhang/.julia/packages/AlgebraicMultigrid/xJj9G/src/multilevel.jl:82
[2] show(::IOContext{REPL.Terminals.TTYTerminal}, ::MIME{Symbol("text/plain")}, ::AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64},GaussSeidel{SymmetricSweep},GaussSeidel{SymmetricSweep},SparseArrays.SparseMatrixCSC{Float64,Int64},LinearAlgebra.Adjoint{Float64,SparseArrays.SparseMatrixCSC{Float64,Int64}},SparseArrays.SparseMatrixCSC{Float64,Int64},AlgebraicMultigrid.MultiLevelWorkspace{Array{Float64,1},1}}) at ./sysimg.jl:194
[3] display(::REPL.REPLDisplay, ::MIME{Symbol("text/plain")}, ::Any) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:131
[4] display(::REPL.REPLDisplay, ::Any) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:135
[5] display(::Any) at ./multimedia.jl:287
[6] #invokelatest#1 at ./essentials.jl:742 [inlined]
[7] invokelatest at ./essentials.jl:741 [inlined]
[8] print_response(::IO, ::Any, ::Any, ::Bool, ::Bool, ::Any) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:155
[9] print_response(::REPL.AbstractREPL, ::Any, ::Any, ::Bool, ::Bool) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:140
[10] (::getfield(REPL, Symbol("#do_respond#38")){Bool,getfield(REPL, Symbol("##48#57")){REPL.LineEditREPL,REPL.REPLHistoryProvider},REPL.LineEditREPL,REPL.LineEdit.Prompt})(::Any, ::Any, ::Any) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:714
[11] #invokelatest#1 at ./essentials.jl:742 [inlined]
[12] invokelatest at ./essentials.jl:741 [inlined]
[13] run_interface(::REPL.Terminals.TextTerminal, ::REPL.LineEdit.ModalInterface, ::REPL.LineEdit.MIState) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/REPL/src/LineEdit.jl:2273
[14] run_frontend(::REPL.LineEditREPL, ::REPL.REPLBackendRef) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:1035
[15] run_repl(::REPL.AbstractREPL, ::Any) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:192
[16] (::getfield(Base, Symbol("##734#736")){Bool,Bool,Bool,Bool})(::Module) at ./client.jl:362
[17] #invokelatest#1 at ./essentials.jl:742 [inlined]
[18] invokelatest at ./essentials.jl:741 [inlined]
[19] run_main_repl(::Bool, ::Bool, ::Bool, ::Bool, ::Bool) at ./client.jl:346
[20] exec_options(::Base.JLOptions) at ./client.jl:284
[21] _start() at ./client.jl:436

I haven't done a rigorous profile, but eyeballing the CPU usage, it seems that the _solve step is single-threaded. Do you plan on adding multithreading? Thanks!

To reproduce

Download this FEMLAB/poisson3Db file and run:

using MatrixMarket
import AlgebraicMultigrid as AMG

A = MatrixMarket.mmread("poisson3Db.mtx")
b = MatrixMarket.mmread("poisson3Db_b.mtx") |> vec

ml = AMG.ruge_stuben(A)
x_amg = AMG._solve(ml, b, verbose=true, maxiter=20)

The residual diverges:

Norm of residual at iteration      1 is 1.8735e+06
Norm of residual at iteration      2 is 2.1074e+06
Norm of residual at iteration      3 is 3.6537e+06
Norm of residual at iteration      4 is 5.8408e+06
Norm of residual at iteration      5 is 9.5234e+06
Norm of residual at iteration      6 is 1.6343e+07
Norm of residual at iteration      7 is 2.9776e+07
Norm of residual at iteration      8 is 5.7516e+07
Norm of residual at iteration      9 is 1.1839e+08
Norm of residual at iteration     10 is 2.7924e+08
Norm of residual at iteration     11 is 9.3488e+08
Norm of residual at iteration     12 is 4.4738e+09
Norm of residual at iteration     13 is 2.4183e+10
Norm of residual at iteration     14 is 1.3371e+11
Norm of residual at iteration     15 is 7.4188e+11
Norm of residual at iteration     16 is 4.1185e+12
Norm of residual at iteration     17 is 2.2864e+13
Norm of residual at iteration     18 is 1.2693e+14
Norm of residual at iteration     19 is 7.0469e+14
Norm of residual at iteration     20 is 3.9122e+15

For a good coarse grid transfer, the column sum of restriction operator R (or equivalently the row sum of prolongation P) should be approximately equal to 1. However, the resulting operator contains sum as large as 6.

R_sum = vec(sum(ml.levels[1].R, dims=1))

maximum(R_sum)  # 6.5, way too big
sum(R_sum .> 1.5)  # 9174

using Plots
plot(R_sum)  # show all sums

Compare with PyAMG default

PyAMG with default settings is able to converge:

import pyamg
from scipy.io import mmread

A = mmread("poisson3Db.mtx").tocsr()
b = mmread("poisson3Db_b.mtx").flatten()

ml = pyamg.ruge_stuben_solver(A)

residuals = []
x_amg = ml.solve(b, residuals=residuals, maxiter=20)
print(residuals)

[1873512.4673384393,
 233255.5159887513,
 109241.08890436463,
 70342.75287986077,
 51906.30512772909,
 41392.049922233826,
 34744.68012516487,
 30219.93538009633,
 26944.726717410056,
 24444.551602043706,
 22448.232978979628,
 20794.899486606348,
 19385.86164700615,
 18158.491431757277,
 17071.605472539315,
 16097.097721476017,
 15215.053123024978,
 14410.837321621537,
 13673.326931758267,
 12993.808948353657,
 12365.279421869809]

The interpolation operators also looks good:

import numpy as np
import matplotlib.pyplot as plt

R_sum = np.squeeze(np.asarray(ml.levels[0].R.sum(axis=0)))
R_sum.max()  #  1.26
np.sum(R_sum > 1.5)  # 0

plt.plot(R_sum, marker='.', linewidth=0, alpha=0.2)  # all nicely bounded

By default, both PyAMG and AMG.jl use theta=0.25 for strength of connection, and symmetric Gauss-Seidel as smoother, so their behaviors should have been similar. Maybe some subtle problems in the interpolation code. Haven't fully figured it out yet.

As described in https://discourse.julialang.org/t/ann-plans-for-removing-packages-that-do-not-yet-support-1-0-from-the-general-registry/ we are planning on removing packages that do not support 1.0 from the General registry. This package has been detected to not support 1.0 and is thus slated to be removed. The removal of packages from the registry will happen approximately a month after this issue is open.

To transition to the new Pkg system using Project.toml, see https://github.com/JuliaRegistries/Registrator.jl#transitioning-from-require-to-projecttoml.
To then tag a new version of the package, see https://github.com/JuliaRegistries/Registrator.jl#via-the-github-app.

If you believe this package has erroneously been detected as not supporting 1.0 or have any other questions, don't hesitate to discuss here or in the thread linked at the top of this post.

It can be helpful to allow custom interpolations to be passed into the solver.

The hierarchy adds one level extra.

A = float.(poisson(10^6))
smoothed_aggregation(A)

Multilevel Solver
-----------------
Operator Complexity: 1.5
Grid Complexity: 1.5
No. of Levels: 11
Coarse Solver: AMG.Pinv()
Level     Unknowns     NonZeros
-----     --------     --------
    1      1000000      2999998 [66.67%]
    2       333334      1000000 [22.22%]
    3       111112       333334 [ 7.41%]
    4        37038       111112 [ 2.47%]
    5        12346        37036 [ 0.82%]
    6         4116        12346 [ 0.27%]
    7         1372         4114 [ 0.09%]
    8          458         1372 [ 0.03%]
    9          153          457 [ 0.01%]
   10           51          151 [ 0.00%]
   11           17           49 [ 0.00%]

The following gives me an error:

ml = ruge_stuben(poisson(27000)+24.0*I)

This results in:
ERROR: ArgumentError: reducing over an empty collection is not allowed
It occurs for sufficiently large values on the diagonal, I'm guessing something goes wrong in the coarsening.

Just an issue to reference in the coming PR.

using Circuitscape, AMG
a = model_problem(Float32, 10)
ml = ruge_stuben(a)
b = rand(Float32, 100)
solve(ml, b)

errors out:

ERROR: MethodError: no method matching gs!(::SparseMatrixCSC{Float64,Int64}, ::Array{Float32,1}, ::Array{Float32,1}, ::Int64, ::Int64, ::Int64)
Closest candidates are:
  gs!(::SparseMatrixCSC{T,Ti}, ::Array{T,1}, ::Array{T,1}, ::Any, ::Any, ::Any) where {T, Ti} at /Users/ranjan/.julia/v0.6/AMG/src/smoother.jl:36
Stacktrace:
 [1] smoother!(::AMG.GaussSeidel{AMG.SymmetricSweep}, ::AMG.ForwardSweep, ::SparseMatrixCSC{Float64,Int64}, ::Array{Float32,1}, ::Array{Float32,1}) at /Users/ranjan/.julia/v0.6/AMG/src/smoother.jl:21
 [2] smoother!(::AMG.GaussSeidel{AMG.SymmetricSweep}, ::AMG.SymmetricSweep, ::SparseMatrixCSC{Float64,Int64}, ::Array{Float32,1}, ::Array{Float32,1}) at /Users/ranjan/.julia/v0.6/AMG/src/smoother.jl:26
 [3] __solve(::AMG.V, ::AMG.MultiLevel{AMG.Pinv,AMG.GaussSeidel{AMG.SymmetricSweep},AMG.GaussSeidel{AMG.SymmetricSweep},Float64,Int64}, ::Array{Float32,1}, ::Array{Float32,1}, ::Int64) at /Users/ranjan/.julia/v0.6/AMG/src/multilevel.jl:97
 [4] #solve#10(::Bool, ::Bool, ::Function, ::AMG.MultiLevel{AMG.Pinv,AMG.GaussSeidel{AMG.SymmetricSweep},AMG.GaussSeidel{AMG.SymmetricSweep},Float64,Int64}, ::Array{Float32,1}, ::Int64, ::AMG.V, ::Float64) at /Users/ranjan/.julia/v0.6/AMG/src/multilevel.jl:83
 [5] solve(::AMG.MultiLevel{AMG.Pinv,AMG.GaussSeidel{AMG.SymmetricSweep},AMG.GaussSeidel{AMG.SymmetricSweep},Float64,Int64}, ::Array{Float32,1}) at /Users/ranjan/.julia/v0.6/AMG/src/multilevel.jl:69

Same with smoothed_aggregation

Must relax operator B with 4 iterations of SymmetricSweep. Maybe can skip based on accuracy needs?

What is the reference for the direct interpolation used in this package? I see that it is an extension to non Z matrices but is there a paper or book that describes the rationale?

Displaying the MultiLevel object result in an error. The Base.show seems to contain some code Julia 0.6 not jet ported to julia 0.7.

using AlgebraicMultigrid
A = poisson(1000);
ml = ruge_stuben(A);
display(ml)

The fix is easy and I will make a PR.

ERROR: MethodError: no method matching round(::Float64, ::Int64)
Closest candidates are:
  round(::Float64, ::RoundingMode{:Nearest}) at float.jl:368
  round(::Float64, ::RoundingMode{:Up}) at float.jl:366
  round(::Float64, ::RoundingMode{:Down}) at float.jl:364
  ...
Stacktrace:
 [1] show(::IOContext{REPL.Terminals.TTYTerminal}, ::AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64},GaussSeidel{SymmetricSweep},GaussSeidel{SymmetricSweep},SparseMatrixCSC{Float64,Int64},Adjoint{Float64,SparseMatrixCSC{Float64,Int64}},SparseMatrixCSC{Float64,Int64},AlgebraicMultigrid.MultiLevelWorkspace{Array{Float64,1},1}}) at ./printf.jl:159
 [2] show(::IOContext{REPL.Terminals.TTYTerminal}, ::MIME{Symbol("text/plain")}, ::AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64},GaussSeidel{SymmetricSweep},GaussSeidel{SymmetricSweep},SparseMatrixCSC{Float64,Int64},Adjoint{Float64,SparseMatrixCSC{Float64,Int64}},SparseMatrixCSC{Float64,Int64},AlgebraicMultigrid.MultiLevelWorkspace{Array{Float64,1},1}}) at ./sysimg.jl:194
 [3] display(::REPL.REPLDisplay, ::MIME{Symbol("text/plain")}, ::Any) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/REPL/src/REPL.jl:131
 [4] display(::REPL.REPLDisplay, ::Any) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/REPL/src/REPL.jl:135
 [5] display(::Any) at ./multimedia.jl:287
 [6] top-level scope at none:0

The only available CoarseSolver is Pinv, which gets very slow if the coarsest-level A is large, for example in a two-grid setting instead of multigrid V-cycle. I wrote a code demo to use UMFPACK sparse LU to solve the coarse equation. It shows significant speed-up in both AMG setup phase and solve phase. I can make a PR if useful.

Current performance with dense Pinv

import Random: seed!
using SparseArrays
using LinearAlgebra
using AlgebraicMultigrid
import AlgebraicMultigrid as AMG
using BenchmarkTools

# Poisson matrix also makes the point, but the second-level A with Ruge-Stuben coarsening is a bit too dense,
# so the gap between sparse and dense solves is not significant enough.
# Here use a very sparse random matrix that lead to a sparse-enough second level.
# A = poisson((64, 64))
seed!(0)
A = sprand(8000, 8000, 1e-4) * 0.1 + I  #  add to diagonal to avoid singular matrix
b = A * ones(size(A)[1])

ruge_stuben(A, max_levels=10)  # fast if coarsest A is very small
@time ml = ruge_stuben(A, max_levels=2)  # takes 5 seconds, mostly spends on `inv(Matrix(ml.final_A))`
@btime AMG._solve(ml, b, maxiter=10);  # takes 16 ms

Use Sparse LU to improve performance

import SuiteSparse.UMFPACK: UmfpackLU
import AlgebraicMultigrid: CoarseSolver

# note: UMFPACK only supports Float64 or ComplexF64
# not sure how to best restrict the types
struct Splu <: CoarseSolver
    LU::UmfpackLU
    Splu(A) = new(lu(A))
end

Base.show(io::IO, p::Splu) = print(io, "Splu")

function (p::Splu)(x, b)
    x[:] = p.LU \ b
end

ml_sp = ruge_stuben(A, max_levels=2, coarse_solver=Splu)  
@btime ruge_stuben(A, max_levels=2, coarse_solver=Splu)  # takes 6 ms, ~1000x faster than before!
@btime AMG._solve(ml_sp, b, maxiter=10);  # takes 9 ms, 2x faster than before

• Side note 1: PyAMG also supports splu as coarse solver. Its code logic is to run the LU factorization in the first iteration of the solve phase, and reuse the LU factors in later iterations. I think a more natural logic is to compute the LU factors in setup phase, so it is easier to separately benchmark factor vs solve.
• Side note 2: (In case someone is interested in doing so) The sparse LU apply stage breaks Zygote autodiff FluxML/Zygote.jl#1168
• Side note 3: Might also be useful to call Pardiso.jl instead of Umfpack.

The second part of the above condition (lambda[i] == 1 && Tj[Tp[i]] == i) will always be false since we dropped the diagonal elements of S earlier so even if there exists one off-diagonal element in the ith column of S, i.e. lambda[i] == 1, this element will not be a diagonal element, i.e. Tj[Tp[i]] == i is false.

The line ground = round(op, digits = 3) below should be ground = round(g, digits = 3) instead:

Package version

[email protected]

I've tried to incorporate AMG into one of my codes, where -- unfortunately -- it fails quite spectacularly:

import JLD, AMG
A, f = JLD.load("AMG_example.jld", "A", "f")
x = AMG.solve(AMG.ruge_stuben(A), f, 1000, AMG.V(), 1e-10)
@show vecnorm(A \ f - x, Inf)
@show vecnorm(A * x - f, Inf)
@show any(isnan.(A\f))
@show vecnorm(A - A', Inf)
@show minimum(eigvals(Symmetric(full(A))))

generates the output

vecnorm(A \ f - x, Inf) = Inf
vecnorm(A * x - f, Inf) = NaN
any(isnan.(A \ f)) = false
vecnorm(A - A', Inf) = 0.0
minimum(eigvals(Symmetric(full(A)))) = 0.009999999999991717

The file needed to reproduce this can be obtained from HERE

The matrix A here is essentially H + \epsilon I. where \epsilon ~ 0.01, and H can be thought of as a weighted graph-laplacian of a two-dimensional network.