First import PyAMG (the others are used later):

import scipy
import numpy
import pyamg
import pylab

Let's look at the 2D Q1 discretization of the Poisson problem. Create a stencil and the matrix

stencil = [[-1,-1,-1],[-1,8,-1],[-1,-1,-1]]
A = pyamg.gallery.stencil_grid(stencil, (1000,1000), dtype=float, format='csr')

Now check out the matrix

A.nnz       # Number of nonzeros
A.shape     # Matrix dimensions
A.tocsr()   # Return A in Compressed Sparse Row (CSR)
A.tocsc()   # Return A in Compressed Sparse Column (CSC)
A.tocoo()   # Return A in Coordinate (COO) format
A.T         # Return the transpose of A

Construct standard aggregation for A. Aggregation is one phase of AMG based on Smoothed Aggregation (SA).

Agg = pyamg.aggregation.standard_aggregation(A)

One can do a lot with this Agg (e.g. visualize, construct full solver). Or we could build a multigrid hierarchy directly

B = numpy.ones((A.shape[0],1))
ml = pyamg.smoothed_aggregation_solver(A,B,max_coarse=10)

This tells PyAMG to create an SA hierarchy using A, a near-null space vector B of ones, and to coarsen until the system has 10 degrees-of-freedom.

Now print the details of the multilevel object:

print ml

which gives

multilevel_solver
Number of Levels:     7
Operator Complexity:  1.125
Grid Complexity:      1.126
Coarse Solver:        'pinv2'
  level   unknowns     nonzeros
    0      1000000      8988004 [88.87%]
    1       111556      1000000 [ 9.89%]
    2        12544       111556 [ 1.10%]
    3         1444        12544 [ 0.12%]
    4          169         1369 [ 0.01%]
    5           25          169 [ 0.00%]
    6            4           16 [ 0.00%]

Here we can see that 7 levels are constructed until a coarse level of 4 (<max_coarse=10) is found. The operator complexity (the total nnz on the fine grid compares to the total nnz on all levels) is moderates low, and the grid complexity (total dof on the finest level compared with the fine level) is low. Also the coarsest level solve set at a pseudo-inverse (pinv2).

Now we can start a places holder for the residuals, initialize a right-hand-side to Ax=b, and give the solver an initial guess x0:

residuals=[]
b = scipy.rand(A.shape[0],1)
x0 = scipy.rand(A.shape[0],1)

Now, let's finally solve to a tolerance of 1e-10

x = ml.solve(b=b,x0=x0,tol=1e-10,residuals=residuals)

How well did we do? Let's look at the geometric convergence factor:

(residuals[-1]/residuals[0])**(1.0/len(residuals))

which gives,

.16

in this case. Not bad. Now with CG acceleration we do

x = ml.solve(b=b,x0=x0,tol=1e-10,residuals=residuals,accel='cg')

which gives a geometric convergence factor of

(residuals[-1]/residuals[0])**(1.0/len(residuals))

0.072

Plotting the residuals

pylab.semilogy(residuals/residuals[0],'o-')
pylab.xlabel('iterations')
pylab.ylabel('normalized residual')
pylab.show()

Summary

import scipy
import numpy
import pyamg
import pylab
stencil = [[-1,-1,-1],[-1,8,-1],[-1,-1,-1]]
A = pyamg.gallery.stencil_grid(stencil, (1000,1000), dtype=float, format='csr')
B = numpy.ones((A.shape[0],1))
ml = pyamg.smoothed_aggregation_solver(A,B,max_coarse=10)
print ml
residuals=[]
b = scipy.rand(A.shape[0],1)
x0 = scipy.rand(A.shape[0],1)
x = ml.solve(b=b,x0=x0,tol=1e-10,residuals=residuals)
(residuals[-1]/residuals[0])**(1.0/len(residuals))
x = ml.solve(b=b,x0=x0,tol=1e-10,residuals=residuals,accel='cg')
(residuals[-1]/residuals[0])**(1.0/len(residuals))
pylab.semilogy(residuals/residuals[0],'o-')
pylab.xlabel('iterations')
pylab.ylabel('normalized residual')
pylab.show()

Result