JuliaLinearAlgebra / AlgebraicMultigrid.jl

Algebraic Multigrid in Julia
Other
117 stars 22 forks source link

Classical strength should be taking max/maxabs over rows not columns #41

Closed mohamed82008 closed 6 years ago

mohamed82008 commented 6 years ago

The following part of the code: https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl/blob/061a5418f04ffac9eeac4566ec00d378bdf1cf72/src/strength.jl#L16 should be taking the maximum or maximum absolute value over the rows of A not columns. The Python version explains that in the comment section https://github.com/pyamg/pyamg/blob/10e44eaf4934374a3483bc85784349a00a8f749b/pyamg/strength.py#L126.

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: https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl/blob/061a5418f04ffac9eeac4566ec00d378bdf1cf72/src/strength.jl#L13 to S = A'. That way the rest of the strength code is doing the right thing, but the strength matrix is really S'. What is returned from the function should then be made S, S' not S', S as in https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl/blob/061a5418f04ffac9eeac4566ec00d378bdf1cf72/src/strength.jl#L37 as these get assigned into S, T respectively in the calling scope. https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl/blob/061a5418f04ffac9eeac4566ec00d378bdf1cf72/src/classical.jl#L33

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.