[Contents] [TitleIndex] [WordIndex] 

Overview of sparse matrices in Scilab

Abstract

In this page, we present the management of sparse matrices in Scilab. We present the direct and iterative methods for solving sparse systems of linear equations. We present how these tools can be used to solve the Poisson PDE.

Contents

  1. Overview of sparse matrices in Scilab
    1. Introduction
    2. Key features
    3. Creating a sparse matrix
    4. Solving linear equations with direct methods
    5. Iterative methods for linear systems of equations
    6. The Imsls toolbox
    7. The Spilu toolbox
    8. Read, Write and Plot sparse matrices
    9. Solving the Poisson equation with sparse matrices
    10. Open Source Libraries
    11. Document and tutorials
    12. Conclusion
    13. Bibliography

Introduction

In numerical analysis, a sparse matrix is a matrix with a large number of zeros. Huge sparse matrices often appear in science or engineering when solving partial differential equations. Fortunately, Scilab only stores the nonzero entries of sparse matrices. Hence it can solve sparse linear systems with hundreds of thousands of unknowns within a few seconds.

Scilab provides several features to manage sparse matrices and perform usual linear algebra operations on them. These operations includes all the basic linear algebra including addition, dot product, transpose and the matrix vector product. Moreover, Scilab can solve sparse linear systems of equations, compute various sparse decomposition and can compute eigenvalues of sparse matrices.

Key features

Scilab provides the following tools for sparse matrices:

Scilab manages several file formats read and write sparse matrices.

Creating a sparse matrix

The "sparse" function creates a sparse matrix from a full matrix. In the following session, we create a 5-by-5 sparse matrix of doubles. 

Afull= [
  2  3  0  0  0;
  3  0  4  0  6; 
  0 -1 -3  2  0; 
  0  0  1  0  0; 
  0  4  2  0  1
];
A = sparse(Afull)

The previous script produces the following output. 

-->Afull= [
-->  2  3  0  0  0;
-->  3  0  4  0  6; 
-->  0 -1 -3  2  0; 
-->  0  0  1  0  0; 
-->  0  4  2  0  1
-->];
-->A = sparse(Afull)
 A  =
(    5,    5) sparse matrix
(    1,    1)        2. 
(    1,    2)        3. 
(    2,    1)        3. 
(    2,    3)        4. 
(    2,    5)        6. 
(    3,    2)      - 1. 
(    3,    3)      - 3. 
(    3,    4)        2. 
(    4,    3)        1. 
(    5,    2)        4. 
(    5,    3)        2. 
(    5,    5)        1.

Notice first that only the nonzero entries are stored. Notice also that the entries are stored row-by-row. This is because Scilab internally uses a format which is very similar to the Compressed Sparse Row (CSR) format.

Solving linear equations with direct methods

The sparse backslash operator "\" solves the sparse linear system of equations Ax=b. In the following example, we solve a small 5-by-5 sparse system. 

Afull= [
  2  3  0  0  0;
  3  0  4  0  6; 
  0 -1 -3  2  0; 
  0  0  1  0  0; 
  0  4  2  0  1
];
A = sparse(Afull);
b = [8 ; 45; -3; 3; 19];
x = A\b

The previous script produces the following output. 

-->Afull= [
-->  2  3  0  0  0;
-->  3  0  4  0  6; 
-->  0 -1 -3  2  0; 
-->  0  0  1  0  0; 
-->  0  4  2  0  1
-->];
-->A = sparse(Afull);
-->b = [8 ; 45; -3; 3; 19];
-->x = A\b
 x  =
    1.  
    2.  
    3.  
    4.  
    5.

Depending on the size of the sparse matrix A, backslash operator has a different behavior.

Iterative methods for linear systems of equations

Scilab provides several iterative methods to solve sparse linear systems of equations:

In the following script, we solve a linear system of equations with the pcg function. 

Afull= [ 
     2 -1  0  0  0
    -1  2 -1  0  0
     0 -1  2 -1  0
     0  0 -1  2 -1
     0  0  0 -1  2
];
A = sparse(Afull);
b = [0 ; 0; 0; 0; 6];
x = pcg(A, b)

The previous script produces the following output. 

-->x = pcg(A, b)
 x  =
    1.  
    2.  
    3.  
    4.  
    5.

pcg(A, b)

The Imsls toolbox

The Imsls toolbox provides iterative methods for sparse linear systems of equations. It includes bicg, bicgstab, pcg, cgs, gmres, qmr, and other solvers as well.

To install this module, we can use the statement 

atomsInstall("imsls")

and restart Scilab.

One of the interesting point here is that the matrix-vector products or the preconditionning steps \scivar{M\textbackslash x} can be performed either with full or sparse matrices, or with callback functions. This flexibility makes the module convenient to use in situations when the sparse matrices are not stored in memory, since only the matrix-vector product (or the preconditionning step \scivar{M\textbackslash x}) is required. Moreover, the toolbox provides Matlab-compatible pcg, gmres and qmr solvers:

This contrasts with Scilab's internal functions, where the two last points are completely unsatisfied. Finally, the toolbox provides a complete test suite for these functions, which are using robust argument checking.

The following figure is a benchmark of various solvers on the 176-by-176 sparse Wathen matrix.

imsls-benchmarkWathen.svg

The Spilu toolbox

Spilu is a Scilab toolbox which provides preconditioners based on Incomplete LU (ILU) factorizations. This module is based on a set of Fortran routines from the Sparskit module by Yousef Saaf. More specifically, this module provides some of the preconditioners from the ITSOL sub-module of Sparskit.

The preconditioners which are provided in this toolbox may be used in preconditioned iterative algorithms for solving sparse linear systems of equations. According to Y. Saad, "roughly speaking, a preconditioner is any form of implicit or explicit modification of an original linear system which makes it easier to solve by a given iterative method." Examples of preconditioned iterative algorithms are the Generalized Minimum Residual Method (GMRES) or the Preconditioned Conjugate Gradient (PCG). Hence, the Spilu toolbox is the companion of the Imsls toolbox, which provides these iterative methods.

To install this module, we can use the statement 

atomsInstall('spilu')

and restart Scilab.

The following figure presents the effect of two ILU preconditionners on the PDE225 matrix with the GMRES iterative method. We see that the GMRES algorithm converges much faster when the ILUTP preconditionning is used.

ILU-GMRES.svg

Read, Write and Plot sparse matrices

Scilab can read, write sparse matrices. First, the ReadHBSparse function reads a Harwell-Boeing sparse file. Second, the Matrix Market module can read and write Matrix Market sparse or dense matrices. To install the Matrix Market module, just run: 

atomsInstall("MatrixMarket")

and restart Scilab. The Matrix Market module provides the following functions:

Moreover, Scilab provides the PlotSparse function, which plots the sparsity pattern of the nonzero entries of a sparse matrix.

In the following script, we read a sparse matrix provided in the Umfpack module with the ReadHBSparse function. Then we plot the sparsity pattern with the PlotSparse function. 

umfdir = fullfile(SCI,"modules","umfpack","examples");
filename = fullfile(umfdir,"arc130.rua");
A = ReadHBSparse(filename);
PlotSparse(A,"y+");

The previous script produces the following figure.

sparse-matrixPlot.svg

Solving the Poisson equation with sparse matrices

Sharma and Gobbert analyzed the performance of Scilab for the resolution of sparse linear systems of equations associated with the Poisson equation [4]. We consider the Poisson problem with homogeneous Dirichlet boundary conditions and are interested in the numerical solution based on finite differences. More precisely, we consider the 2 dimensional Partial Differential Equation:

$$\begin{array}{rlll}- \Delta u &=& f &\textrm{ in the domain,} \\ u &=& 0 &\textrm{ on the frontier.} \end{array}$$

where the two dimensionnal Laplace operator is

$$\Delta u = \frac{\partial^2 u}{dx^2} + \frac{\partial^2 u}{dy^2}$$

We consider the domain $$0 \leq x \leq 1$$, $$0 \leq y \leq 1$$. The function f is defined by

$$f(x, y) = -2\pi^2 cos(2\pi x) sin^2(\pi y) - 2\pi^2 sin^2(\pi x) cos(2\pi y) $$

The solution is

$$u(x, y) = sin^2(\pi x) sin^2(\pi y)$$

We use a second order finite difference approximation of the Laplace operator based on a grid of N-by-N points.

To solve this problem, we use the Scibench external module. In order to install this module, we use the statement: 

atomsInstall("scibench")

and restart Scilab.

We use the "poisson.sce" script provided at:

http://forge.scilab.org/index.php/p/scibench/source/tree/HEAD/demos/poisson.sce

The scibench_poissonA(n) statement creates a sparse matrix equation associated with n cells in the X coordinate and n cells in the Y coordinate. Before calling this session, we call the stacksize function in order to let Scilab allocate as much memory as possible. Then we call the PlotSparse function to plot the sparsity pattern of the matrix. 

stacksize("max");
A = scibench_poissonA(50);
PlotSparse(A)

The previous script produces the following plot.

PoissonSparsityPattern.png

The scibench_poisson function solves the 2D Poisson equation. Its calling sequence is 

scibench_poisson ( N , plotgraph , verbose , solver )

where N is the number of cells, plotgraph is a boolean to plot the graphics, verbose is a boolean to print messages and solver is a function which solves the linear equation A*x=b.

The solver argument is designed so that we can customize the linear equation solver which we want to use. For example, if we want to use the sparse backslash operator, all we have to do is to create the mysolverBackslash function as below. 

function u=mysolverBackslash(N, b)
    A = scibench_poissonA(N);
    u = A\b;
endfunction

The following script solves the Poisson equation with N = 50. 

scf();
scibench_poisson(50, %t , %t , mysolverBackslash );

The previous script produces the following output, where h is the dimensionnal spacee step and enorminf is the value of the infinite norm of the error. 

-->scibench_poisson(50, %t , %t , mysolverBackslash );
N =    50
h =  1.9607843137254902e-002
h^2 =  3.8446751249519417e-004
enorminf =  1.2634082165868810e-003
C = enorminf / h^2 =  3.2861247713424775e+000
wall clock time =       0.15 seconds

The previous script also produces the following plot.

plotSolutionPoisson.png

We could use other Scilab functions to solve this problem, including the Umfpack module or iterative solvers (such as the pcg function). Interested readers may read [3] for more details on this topic.

Open Source Libraries

Document and tutorials

A complete review of sparse matrices in Scilab is presented at:

In "CPSC 5006 EL01 Matrix Computations" (HTML), Julien Dompierre, from the Department of Mathematics and Computer Science, Laurentian University, Canada, presents a set of lectures and Scilab scripts for Matrix computations. Some of these lectures are on sparse matrices:

Conclusion

Scilab 5 can manage sparse matrices and solve partial differential equations such as the Poisson equation for example. Indeed, Scilab provides several sparse linear equation solvers, including a sparse backslash operator, iterative methods (e.g. the pre-conditionned conjugate gradient algorithm) and other direct solvers such as the Umfpack or Taucs modules. With these tools we can solve huge systems of linear equations, because Scilab only stores the nonzero entries of these massive matrices.

More details on sparse matrices are presented in [1].

Bibliography

2022-09-08 09:27