[pyamg-user] Coarsening by Compatible Relaxation

[Florian Wilhelm]

Having read O.E.Livne's paper "Coarsening by compatible relaxation", I
stumbled upon PyAMG which implements Livne's algorithm. Is it possible
right know in PyAMG to combine CR coarsening with a multilevel solver?
According to Livne's conclusions it should be possible to use it as a
"stand-alone" module that can be incorporated into existing AMG
software/algorithms. As far as I can tell this is not possible in
general because for instance classical AMG's prolongation/
interpolation makes some assumption about CF connections that CR
generated coarse grids may not fulfill. This means that in order to
use it, I'd have to add an interpolation operator. Am I right about
that?

--
You received this message because you are subscribed to the Google Groups "pyamg-user" group.
To post to this group, send email to pyamg...@googlegroups.com.
To unsubscribe from this group, send email to pyamg-user+...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/pyamg-user?hl=en.

[Luke]

Florian,

CR is indeed possible in a multilevel solver. In fact, we have an
example in Examples/CompatibleRelaxation of creating a splitting.
The idea is that
splitting = CR(A)
is called and fed to a classical AMG solver. As you point out, there
may be some limitations with interpolation in classical AMG. We have
direct interpolation implemented in pyamg/classical/interpolate.py,
but it would be straightforward to extend this to more general
settings. (In the solver below, direction interpolation works fine in
a two level method).

Also, two cautions with the state of CR in PyAMG right now. The first
is that only 'habituated' CR is implemented. There is a placeholder
for a 'concurrent' version. The second is that we haven't profiled
this code heavily so there may be performance issues. For example, I
suspect that the biorthogonalization would benefit from some C++
optimization.

I've added to two tickets for these issues. We're targeting a new
release soon and I'll try to put these in.

Let me know if you run into problems, but the framework is all there
to analyze a multilevel CR. I'm pasting an example called demo2.py.
To run it you need to make one change to classical.py so that CR can
be called. Add this:
elif CF in ['CR']:
import cr
splitting = cr.CR(A,method='habituated',maxiter=20)
after the first if CF statement here:
http://code.google.com/p/pyamg/source/browse/trunk/pyamg/classical/classical.py#109

One issue with the example is that it hangs for more than two levels.
I'm not sure why this is happening, but the two level method works
fine.

Luke

%%%%
import numpy
import scipy

from pyamg.gallery import poisson
from pyamg.classical import CR
from pyamg.vis import vis_splitting
from pyamg import ruge_stuben_solver

print '----setting up problem'
n = 25
A=poisson((n,n)).tocsr()

xx = numpy.arange(0,n,dtype=float)
x,y = numpy.meshgrid(xx,xx)
Verts = numpy.concatenate([[x.ravel()],[y.ravel()]],axis=0).T

print '----setting up RS solver'
mls = ruge_stuben_solver(A, max_levels=2, max_coarse=5, CF='CR')
print mls

print '----running RS solver'
resvec = []
b = numpy.zeros((A.shape[0],))
z = x = scipy.rand(A.shape[0])
z = mls.solve(b, x0=z, maxiter=200, tol=1e-8, residuals=resvec)
avgconvergencefactor = (resvec[-1]/resvec[0])**(1.0/len(resvec))
allconvergencefactor = numpy.array(resvec[1::])/
numpy.array(resvec[0:-1:])
print avgconvergencefactor
print allconvergencefactor

print '----setting up print'
from pylab import figure, axis, scatter, show
splitting = mls.levels[0].splitting
C_nodes = splitting == 1
F_nodes = splitting == 0
figure(figsize=(6,6))
axis('equal')
scatter(Verts[:,0][C_nodes], Verts[:,1][C_nodes], c='r', s=100.0)
#plot C-nodes in red
scatter(Verts[:,0][F_nodes], Verts[:,1][F_nodes], c='b', s=100.0)
#plot F-nodes in blue
show()
%%%%

[Florian Wilhelm]

Luke,

thanks for your reply! I added the lines in classical.py and ran the
script, which lead to following error message:

----setting up problem
----setting up RS solver
Traceback (most recent call last):
File "cr-solver.py", line 18, in <module>

mls = ruge_stuben_solver(A, max_levels=2, max_coarse=5, CF='CR')

File "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/
lib/python2.6/site-packages/pyamg/classical/classical.py", line 76, in
ruge_stuben_solver
extend_hierarchy(levels, strength, CF)
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/
lib/python2.6/site-packages/pyamg/classical/classical.py", line 116,
in extend_hierarchy
P = direct_interpolation(A, C, splitting)
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/
lib/python2.6/site-packages/pyamg/classical/interpolate.py", line 55,
in direct_interpolation
C.indptr, C.indices, splitting, Pp)
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/
lib/python2.6/site-packages/pyamg/amg_core/amg_core.py", line 83, in
rs_direct_interpolation_pass1
return _amg_core.rs_direct_interpolation_pass1(*args)
TypeError: Array of type 'int' required. Array of type 'long' given

I'm running pyamg revision 659 and python 2.6 + pyamg requirements
installed on MacOS X 10.6 via macports.
When printing the args given to rs_direct_interpolation_pass1 I found
that the splitting array was missing dtype=int32.
Following changes then fixed the problem for me:

elif CF in ['CR']:
import cr

import numpy
splitting =
numpy.array(cr.CR(A,method='habituated',maxiter=20),dtype='int32')

Florian

> after the first if CF statement here:http://code.google.com/p/pyamg/source/browse/trunk/pyamg/classical/cl...

[Luke Olson]

Florian,

Ok, I couldn't reproduce it on my mac, but I changed int to 'intc' in cr.py and committed it. See if that fixes it.
http://code.google.com/p/pyamg/source/diff?spec=svn660&r=660&format=side&path=/trunk/pyamg/classical/cr.py&old_path=/trunk/pyamg/classical/cr.py&old=545

I'll hold off on the modification to classical.py until we get multilevel cr working.

Luke

[Florian Wilhelm]

Luke,

yeap, that works for me now without my modifications. Thanks.

Regards,

Florian

On 29 Apr., 18:12, Luke Olson <luke.ol...@gmail.com> wrote:
> Florian,
>

> Ok, I couldn't reproduce it on my mac, but I changed int to 'intc' in cr.py and committed it.  See if that fixes it.http://code.google.com/p/pyamg/source/diff?spec=svn660&r=660&format=s...