Matlab

SPARSE MATRICES IN MATLAB: DESIGN AND IMPLEMENTATION
JOHN R. GILBERT , CLEVE MOLERy , AND ROBERT SCHREIBERz

Dedicated to Gene Golub on the occasion of his 60th birthday.

Abstract. We have extended the matrix computation language and environment Matlab to include sparse matrix storage and operations. The only change to the outward appearance of the Matlab language is a pair of commands tocreate full or sparse matrices. Nearly all the operations of Matlab now apply equally to full or sparse matrices, without any explicit action by the user. The sparse data structure represents a matrix in space proportional to the number of nonzero entries, and most of the operations compute sparse results in time proportional to the number of arithmetic operations on nonzeros. Key words. Matlab,mathematical software, matrix computation, sparse matrix algorithms. AMS subject classi cations. 65{04, 65F05, 65F20, 65F50, 68N15, 68R10.

1. Introduction. Matlab is an interactive environment and programming language for numeric scienti c computation 18]. One of its distinguishing features is the use of matrices as the only data type. In Matlab, a matrix is a rectangular array of real or complexnumbers. All quantities, even loop variables and character strings, are represented as matrices, although matrices with only one row, or one column, or one element are sometimes treated specially. The part of Matlab that involves computational linear algebra on dense matrices is based on direct adaptations of subroutines from Linpack and Eispack 5, 23]. An m n real matrix is stored as a fullarray of mn oating point numbers. The computational complexity of basic operations such as addition or transposition is proportional to mn. The complexity of more complicated operations such as triangular factorization is proportional to mn2 . This has limited the applicability of Matlab to problems involving matrices of order a few hundred on contemporary workstations and perhaps a few thousand oncontemporary supercomputers. We have now added sparse matrix storage and operations to Matlab. This report describes our design and implementation. Sparse matrices are widely used in scienti c computation, especially in largescale optimization, structural and circuit analysis, computational uid dynamics, and, generally, the numerical solution of partial di erential equations. Several e ectiveFortran subroutine packages for solving sparse linear systems are available, including Sparspak 11], the Yale Sparse Matrix Package 9], and some of the routines in the Harwell Subroutine Library 25]. Our work was facilitated by our knowledge of the techniques used in the Fortran sparse matrix packages, but we have not directly adapted any of their code. Matlab
Xerox Palo Alto Research Center, 3333Coyote Hill Road, Palo Alto, California 94304. The MathWorks, 325 Lin eld Place, Menlo Park, California 94025. Research Institute for Advanced Computer Science, MS T045-1, NASA Ames Research Center, Mo ett Field, CA 94035. This author's work was supported by the NAS Systems Division and DARPA via Cooperative Agreement NCC 2-387 between NASA and the University Space Research Association (USRA).Copyright c 1991 by Xerox Corporation, Research Institute for Advanced Computer Science, and The MathWorks Incorporated. All rights reserved. 1
y z

Operations with the 4096 by 4096 discrete Laplacian.

Table 1

Sparse Full Memory 0.25 megabyte 128 megabytes Compute Dx 0.2 seconds 30 seconds Solve Dx = b 10 seconds > 12 hours is written in C and we wished to take advantage of the data structuresand other programming features of C that would not be used in a simple translation of Fortran code. We also wanted to implement the full range of matrix operations that Matlab provides the Fortran packages do not generally have routines for simply adding or transposing sparse matrices, for example. And, nally, we wanted to incorporate some recent algorithmic ideas that are not used in the...