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

Overview of sparse matrices in Scilab
 Introduction
 Key features
 Creating a sparse matrix
 Solving linear equations with direct methods
 Iterative methods for linear systems of equations
 The Imsls toolbox
 The Spilu toolbox
 Read, Write and Plot sparse matrices
 Solving the Poisson equation with sparse matrices
 Open Source Libraries
 Document and tutorials
 Conclusion
 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:
 basic arithmetic on sparse matrices,
 sparse LU decomposition and resolution of linear equations from this decomposition,
 sparse Cholesky decomposition and resolution of linear equations from this decomposition,
 iterative methods of linear systems of equations,
 incomplete LU (ILU) decompositions,
 sparse LU decomposition and resolution of linear equations from this decomposition (based on the Umfpack library),
 sparse Cholesky decomposition and resolution of linear equations from this decomposition (based on the Taucs library),
 eigenvalues of sparse matrices.
Scilab manages several file formats read and write sparse matrices.
 The ReadHBSparse function reads matrices in the HarwellBoeing format.
 The Matrix Market external module (available in ATOMS) provides functions to read and write matrices in the Matrix Market format.
Creating a sparse matrix
The "sparse" function creates a sparse matrix from a full matrix. In the following session, we create a 5by5 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 rowbyrow. 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 5by5 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.
 If the sparse matrix A is square, then the sparse backslash uses a sparse LU decomposition and uses this decomposition to solve the Ax=b equation.
 If the sparse matrix A is not square, then the normal equations is first used to create a square system of linear equations. Then the sparse LU decomposition is used to solve the system.
Iterative methods for linear systems of equations
Scilab provides several iterative methods to solve sparse linear systems of equations:
 pcg: Preconditionned Conjugate Gradient,
 gmres : Generalized Minimum Residual method,
 qmr : Quasi Minimal Residual method with preconditionning.
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 matrixvector 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 matrixvector product (or the preconditionning step \scivar{M\textbackslash x}) is required. Moreover, the toolbox provides Matlabcompatible pcg, gmres and qmr solvers:
 the order of the arguments are the same as in Matlab,
 the default values of the Matlab functions are the same as in Matlab,
 the headers of the callback functions are the same as in Matlab.
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 176by176 sparse Wathen matrix.
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 submodule 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.
Read, Write and Plot sparse matrices
Scilab can read, write sparse matrices. First, the ReadHBSparse function reads a HarwellBoeing 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:
 mminfo: Extracts size and storage information
 mmread: Reads a Matrix Market file
 mmwrite: Writes a sparse or dense matrix
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.
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:
where the two dimensionnal Laplace operator is
We consider the domain , . The function f is defined by
The solution is
We use a second order finite difference approximation of the Laplace operator based on a grid of NbyN 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.
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.9607843137254902e002 h^2 = 3.8446751249519417e004 enorminf = 1.2634082165868810e003 C = enorminf / h^2 = 3.2861247713424775e+000 wall clock time = 0.15 seconds
The previous script also produces the following plot.
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
The sparse LU decomposition is based on the Sparse package written by Kenneth S. Kundert and Alberto SangiovanniVincentelli, http://www.netlib.org/sparse/.
 The sparse Cholesky decomposition is based on the package developed by Esmond Ng and Barry Peyton at ORNL and a multiple minimundegree ordering package by Joseph Liu at University of Waterloo.
The Umfpack library was developped by Timothy Davis, from the University of Florida. http://www.cise.ufl.edu/research/sparse/umfpack
The Taucs library was developed by Sivan Toledo from TelAviv University, http://www.tau.ac.il/~stoledo/taucs/
Document and tutorials
A complete review of sparse matrices in Scilab is presented at:
"Introduction to sparse matrices in Scilab", Consortium Scilab  Digiteo, Michael Baudin, 20082011 (HTML)
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:
Krylov methods and GMRES Scilab diary. (HTML)
gmres0.sci Scilab function. (HTML)
Conjugate gradient method Scilab diary. (HTML)
SteepestDescent.sci Scilab function. (HTML)
ConjugateGradients.sci Scilab function. (HTML)
Vertex reordoring Scilab diary. (HTML)
Breadth First Search algorithm : BFS_reordering.sci Scilab function. (HTML)
CuthillMcKee algorithm CMK_reordering.sci Scilab function. (HTML)
Basic iterative methods Scilab diary. (HTML)
Jacobi.sci Scilab function. (HTML)
Jacobi_MatrixForm.sci Scilab function. (HTML)
GaussSeidel.sci Scilab function. (HTML)
SOR.sci Scilab function. (HTML)
Sparse storage Scilab diary. (HTML)
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 preconditionned 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
[1] "Introduction to sparse matrices in Scilab", Consortium Scilab  Digiteo, Michael Baudin, 20082011 (HTML)
[2] "Sparse matrix  Wikipedia, The Free Encyclopedia", "http://en.wikipedia.org/w/index.php?title=Sparse_matrix&oldid=364961772"
[3] "Solving Poisson PDE with Sparse Matrices", Consortium Scilab  Digiteo, Michael Baudin, 2011 Solving Poisson PDE with Sparse Matrices
[4] "A Comparative Evaluation Of Matlab, Octave, Freemat, And Scilab For Research And Teaching", 2010, Neeraj Sharma and Matthias K. Gobbert, http://userpages.umbc.edu/~gobbert/papers/SharmaGobbertTR2010.pdf
[5] "Proposition d'ARC  SpaSciAl, Sparse linear Algebra in Scilab", 2002, http://graal.enslyon.fr/~jylexcel/scilabsparse/arc/arc.ps.gz
[6] "Scilab Sparse meeting", 07/12/2001, http://graal.enslyon.fr/~jylexcel/scilabsparse/meeting07/
[7] "Les matrices creuses dans Scilab", Serge Steer, 07/12/2001, http://graal.enslyon.fr/~jylexcel/scilabsparse/meeting07/steer.ps.
[8] "Scilin : resoudre des systemes lineaires avec le logiciel Scilab", Jocelyne Erhel, 07/12/2001, http://graal.enslyon.fr/~jylexcel/scilabsparse/meeting07/scilin0901.pdf
[9] "Direct solvers, Scilab to Scilab //", JeanYves L'Excellent, 07/12/2001, http://graal.enslyon.fr/~jylexcel/scilabsparse/meeting07/jy.ps
[10] "An ordering package for Scilab", JeanPaul BOUFFLET, 07/12/2001, http://graal.enslyon.fr/~jylexcel/scilabsparse/meeting07/sl_meeting_071201.ps.gz