Lower Triangle Sparse Matrix Algorithm

Material Information

Lower Triangle Sparse Matrix Algorithm
Li, Pei Li
Davis, Timothy ( Mentor )
Place of Publication:
Gainesville, Fla.
University of Florida
Publication Date:


serial ( sobekcm )

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.


This item has the following downloads:

Full Text

JOurn.31 or in.lnerr.3dua.3E- -ese-arch

..Oluin-e 4, issue 3 - 2No.en-i1er 2i1:112

Lower Triangle Sparse Matrix Algorithm

Pei Li Li


The project is on sparse matrix algorithm. Given the problem Lx = b, solve for x. L is a lower triangle sparse

matrix and b and x are sparse n-by-1 column vectors. The program is written in MATLAB (version 6)

programming environment. MATLAB stores sparse matrices in three column arrays to save memory space

and enables the program to solve for x more efficiently. The optimal algorithm requires a depth first traversal of

the graph of L. Perform depth first traversal of graph L to find a pattern of x, where the non-zero solution of x

would be. Since L and b are sparse, the solution x will also be sparse. Solving for x in this manner

saves computational time if the matrix is large.

The program is MEX-file written in ANSI C and has a MATLAB interface. The program is callable from

MATLAB command window. The MEX-file contains two routines, a computational routine that performs the

function and a gateway routine that interfaces to MATLAB. The main part of the computational routine

involves finding the pattern of x.


Matrix computations are necessary to solve many mathematical models that arise from practical problems.

As technology advances, in particular Internet technology, more data is produced and causes matrix problems

to increase in size. Research is actively conducted to discover better ways to store and interpret the data.

Many practical problems produce a system of linear algebraic equation. As the size of these problems A

increases, more efficient algorithms are needed to solve for the solution in a timely manner. A system of n

linear algebraic equations:

aIX1 + ai.X + 813X3 + .... + alnXn = bi
811XI + 822X2 + 823X3 + .... + xn = b2
a31x1 + a3�X2 + a33x3 + .... + a3nXn = b3

SnilXI + an2x2 + Hn1X3 + .... + nnXn = b,

is represented as Ax = b where aij ( 1 _ i _ n, 1 _ j _ n ) is the numerical value stored in row i and column j

of matrix A, and bi ( 1 _ i _ n ) is stored in column vector. The solution x result in a n-by-1 column vector. In

many practical problems, A contains a small number of non-zero entries and is sparse. Sparse matrices have few

non-zero entries relative to the size of n^2 entries.

The project is concerned with solving Ax=b when A takes the structure of a lower triangle sparse matrix. A

lower triangle matrix has all the non-zero entries on and below the diagonal. Replace the variable A with L to

get Lx=b where L is a sparse lower triangle matrix, and x and b are n-by-1 sparse column vectors. Let Lij denote

the entry in row i and column j of L. Lij * 0, ( i = j ) and Lij 4 0 when i > j. The optimal algorithm to solve Lx =

b was discovered in mid 80s (1). It requires a graph representation of L. The solution provided is a MEX-file written

in ANSI C in MATLAB (version 6) programming environment.


The program is written in MATLAB programming environment. MATLAB has built-in functions to perform

matrix operations. MATLAB is also a high level programming language that enables users to write m-files, which

are programs that behave like MATLAB built-in functions. MATLAB Application Program Interface enables MATLAB

to call functions written in C or Fortran. These programs are referred to as MEX-file. A MEX-file consists of

two routines, a computational routine and a gateway routine. The computational routine is the source code

that performs the operation of the program. The MEX-function calls the computational routine and provides

the interface to MATLAB. All MEX-functions must contain the heading "void mexFunction ( int nlhs, mxArray *plhs

[ ], int nrhs, const mxArray *prhs [ ] )". The variable nlhs refers to the number of left hand parameters passed to

the function, and plhs is a pointer to all the parameters passed to the function. The variable nrhs refers to

the number of right hand parameter passed back to MATLAB and prhs is the pointer to all of the parameter

passed back to MATLAB. Once a MEX-function is compiled, it behaves similar to MATLAB built-in functions. A

MEX-function is called by the function name.

MATLAB program stores all data types in the object MATLAB array. MATLAB array is declared as type mxArray in

C. The non-zero entries of sparse matrix are stored in three arrays. Pointer pr points to the array that contains

the numerical value of Lij and pointer ir points to the array that contains the corresponding row index. The third

array is referenced by pointer jc. The array stores 0 at jc[0], and the value of jc[i+l] is the value of jc[i] plus

the number of nonzero entries of column i.


A direct method to solve Lx = b takes time complexity 0 (n). The non-optimal algorithm below displays the steps

to solve for x.

1. x = b

2. for i = 1 to n

3. if xi O 0

4. xi = xi/ Li

xi+l:n = xi+l:n - Li+l:n, i * xi

Simple non-optimal algorithm (1)

The optimal algorithm reduces the time complexity to be proportional to the number of binary floating-

point operations. Hence the more non-zero entries L and b contain, the more time it takes to solve for x.

XX = { i xi 4 0 } denote the pattern of x

B = { i | bi * 0 } denote the pattern of b

Lj = { i I Lij 0 } denote the pattern of column j of L (1).

From the matrix of L, graph G (L) was constructed. Let G = (V, E) where V is a vertex or node and E is a direct

edge between two nodes. The nodes of G (L) are integers 1 to n. For each column i of L ( 1 �< i �< n ), connect node

i directly to Lj. Take for example i = 3; if L33, L34 and L35 are non-zero, a direct edge connects node 3 to node 4

and a direct edge connects node 3 to node 5.

Traverse the graph of L to find the pattern of x. XX contains where x has nonzero values. Change line 2 of

algorithm 1 to " for i eXX " provides the optimal algorithm.

Since L and b are sparse, observe a couple of properties that permit XX to be found. bi * 0 implies xi * 0, \I. xi *

0 implies xj * 0, Id and 'vJ e Li (1). Traverse the graph of L. Start at vertex Bo0 and process Bo. Then process all

the nodes reachable from B0. At each node i process all nodes reachable from node i. If there is no vertex

reachable from node i store i in XX and proceed to process the previous node. Once B0 is stored in XX, proceed

to process all the elements of B in a similar manner.

The main part of the computational routine is to find the pattern of x. Due to the storage convention of the graph

of L, each time a vertex i is processed the column i of L is also processed. As a result, a node can be processed

more than once. An array state[n] is created to ensure each node is only process once. State[ ] is initialize to null

at the beginning of the program. Once state [i] is processed, state [i] is changed to 1. Depth-first search method

is implemented with a stack. Since MATLAB stores all data in arrays, the program treats an array as a stack

using Last In First Out principle. After XX is found, the binary floating-point operations are carried to produce

the solution x. Since L and b are sparse, x will also be sparse.

The MEX-function takes in L and b and passes them to the computational routine. The solution is passed back to

MEX-routine and is displayed in MATLAB command window. MATLAB has a built-in function to solve Ax = b and

hence it is used to check the program's accuracy. In MATLAB the command ">> (LA(-1)) * b" provides the solution x.


The author would like to thank her advisor, Dr. Timothy Davis for his guidance, time and encouragement

toward successful completion of this project. She would also like to thank Qi Li Li for proofreading this report.


1. Davis, Tim. "COT 5404 &endash; Analysis and Design of Algorithm",

COT5404 lc Note_Set_l.pdf, University of Florida (as of 2 Aug 2001).

2. Davis, Tim. "Sparse Matrix Algorithms Research at the University of Florida",

sparse/, University of Florida (as of 7 July 2001)

3. Davis, Timothy A. "UMFPACK Version 3.0 User Guide",

UMFPACK3/userGuide.pdf, University of Florida (as-of 2 Aug 2001).

4. Gilberg, Richard and B.A. Forouzan. Data Structure A Pseudocode Approach with C, Boston, MA: PWS Publishing

Company (1998).

5. Lutkepohl, H. Handbook of Matrices, New York, NY: John Wiley & Son, Inc. (1996).

6. Ritchie, M. Dennis and Brian W. Kernighan. The C Programming Language, Englewood Cliffs, New Jersey:

Prentice Hall P T R Prentice Hall, Inc. (1988).

7. Waite, M. and S. Prata and D. Martin. C Primer Plus, Indianapolis, Indiana: Howard W. Sons &Co., Inc. (1984).


Back to the Journal of Undergraduate Research

College of Liberal Arts and Sciences I University Scholars Program I University of Florida |

� University of Florida, Gainesville, FL 32611; (352) 846-2032.

'Pic Fohi,,4ioni f6'r 'Thr Guwo N*ofswff