Link to figure in 2nd paragraph; mention PyTrilinos; insert new code.
PyAMG is attributed to Bell et al.
"restriction array" -> "restriction matrix"

Although the code uses arrays, matrices would perhaps make more
sense, and simplify the code somewhat. But I don't have time right now
to go back and make the changes, and ensure correctness. Tests are
passing right now, with the arrays-based version.
A little more information to make Fig. 3 more intuitively
understandable.
Relate R and A lists to the R_k and A_k matrices.
changed len() to .size

This is actually less likely to give wrong results.
Get dimensionality consistently by the length of the shape tuple.