Tutorial

This is a short PyAMG tutorial, and a good place to start familiarizing yourself with the package. First import PyAMG and supporting packages:

import scipy
import numpy
import pyamg
import matplotlib.pyplot as plt

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

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

From this, print out a few properties of the matrix:

print(f'Number of nonzeros: {A.nnz} \n'
      f'Shape of the matrix: {A.shape[0]} x {A.shape[1]} \n'
      f'Format of the matrix: {A.format}\n')

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 can build a multigrid hierarchy directly with

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 seven levels are constructed until a coarse level of 4 unknowns (max_coarse=10) is achieved. The operator complexity (i.e., the sum of the nnz in the operator over all levels divided by the nnz in the fine level operator) is moderately low, and the grid complexity (i.e., the sum of the number of unknowns on all levels divided by the number of unknowns on the fine level) is low. Also, the coarsest level solver is the pseudo-inverse (pinv2).

Now, we can initialize a place-holder for the residual history, a right-hand-side for Ax=b, and an initial guess x0:

residuals = []
b = numpy.random.rand(A.shape[0])
x0 = numpy.random.rand(A.shape[0])

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:

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

which gives (approximately),

.12

in this case. Then with CG acceleration we add accel='cg':

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

which gives a geometric convergence factor of

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

which gives (approximately)

0.05

Plotting the residuals

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

Summary

import scipy
import numpy
import pyamg
import matplotlib.pyplot as plt

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

print(f'Number of nonzeros: {A.nnz} \n'
      f'Shape of the matrix: {A.shape[0]} x {A.shape[1]} \n'
      f'Format of the matrix: {A.format}\n')

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

residuals = []
b = numpy.random.rand(A.shape[0])
x0 = numpy.random.rand(A.shape[0])

x = ml.solve(b=b, x0=x0, tol=1e-10, residuals=residuals)
print((residuals[-1]/residuals[0])**(1.0/(len(residuals)-1)))

x = ml.solve(b=b, x0=x0, tol=1e-10, residuals=residuals, accel='cg')
print((residuals[-1]/residuals[0])**(1.0/(len(residuals)-1)))

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

Result