PySparse code examples

Main

Modules

Examples

Project Page

Download

Building a 1D Poisson matrix

This example illustrates how to create sparse matrices and how to manipulate their elements. Both functions create the one-dimensional Poisson matrix with twos on the diagonal and minus ones on the off-diagonals. poisson1d creates the matrix in general linked list format and poisson1d_sym creates it in symmetric linked list format.

def poisson1d(n):
    L = spmatrix.ll_mat(n, n)
    for i in range(n):
        L[i,i] = 2
        if i > 0:
            L[i,i-1] = -1
        if i < n-1:
            L[i,i+1] = -1
    return L

def poisson1d_sym(n):
    L = spmatrix.ll_mat_sym(n)
    for i in range(n):
        L[i,i] = 2
        if i > 0:
            L[i,i-1] = -1
    return L

Building a 2D Poisson matrix

def poisson2d_sym_blk(n):
    L = spmatrix.ll_mat_sym(n*n)
    I = spmatrix.ll_mat_sym(n)
    for i in range(n):
        I[i,i] = -1
    P = spmatrix.ll_mat_sym(n)
    for i in range(n):
        P[i,i] = 4
        if i > 0: P[i,i-1] = -1
    for i in range(0, n*n, n):
        L[i:i+n,i:i+n] = P
        if i > 0: L[i:i+n,i-n:i] = I
    return L

Solving a 2D-Poisson system iteratively

Here a 2D-Posson system is solved using the Conjugate Gradient method together with a SSOR preconditioner.

import Numeric, itsolvers, poisson, precon

A = poisson.poisson2d_sym(300)   # the system matrix
A = A.to_sss()                   # convert matrix to SSS format for better performance
n = A.shape[0]
K = precon.ssor(A)               # construct SSOR preconditioner
x = Numeric.zeros(n, 'd')        # the initial guess
b = Numeric.ones(n, 'd')         # the right hand side
info, iter, relres = itsolvers.pcg(A, b, x, 1e-6, 1000, K)
print 'Linear system solved in %d iterations to a residual of  %e' % (iter, relres)

Solving a Poisson system using a direct method

Here the 2D-Poisson system is solved by factorising the system matrix, followed by back- and forward substitution.

import Numeric, poisson, superlu

A = poisson.poisson2d_sym(300)
n = A.shape[0]
x = Numeric.zeros(n, 'd')        # the solution
b = Numeric.ones(n, 'd')         # the right hand side
Alu = superlu.factorize(A.to_csr(), diag_pivot_thresh=0.0)
Alu.solve(b, x)

Computing the five smallest eigenvalues of a 2D Poisson matrix

In this example a 2D-Poisson matrix of order 200² is constructed together with an appropriate SSOR preconditioner. The eigensolutions are computed using the JDSYM eigenvalue solver. The target value is set to 0.0, so the algorithm will converge to the smallest eigenvalues.

import itsolvers, poisson, precon, jdsym

A = poisson.poisson2d_sym(200).to_sss()
K = precon.ssor(A)
k_conv, lmbd, Q, it  = \
    jdsym.jdsym(A, None, K, 5, 0.0, 1e-10, 150, itsolvers.qmrs, clvl=1)