JCL
Library for Reducing Bandwidth and Profile of Sparse Matrix

Feb. 11, 2008

Copyright © 2006-2008
By Equation Solution

LIST OF CONTENTS
Chapter 1. Introduction ................................................................................................. 1
1.1 Non-zero fill-ins and Non-zero coefficients ........................................................... 2
Chapter 2. Calling Syntax ............................................................................................. 3
2.1 A Fortran Example .................................................................................................. 4
2.2 A Bad Example ....................................................................................................... 7
Chapter 3. Renumbering Before the Determination of Sparse Matrix ......................... 10
Chapter 4. Renumbering Asymmetric Sparse Matrix ................................................... 14
Chapter 5. Data Storage Schemes .................................................................................. 15
5.1 Constant Bandwidth Sparse Symmetric Matrix ...................................................... 15
5.2 Variable Bandwidth Sparse Symmetric Matrix ...................................................... 16-i-Chapter 1. Introduction
JCL is a Fortran library for reducing bandwidth and profile of sparse matrix. JCL isalso callable by C language. JCL applies the algorithm, “Algorithms for reducing thebandwidth and profile of a sparse matrix”, authored by J.-C. Luo, published incomputers & structures, vol. 44, no. 3, 1992. The algorithm has been proved to be themost reliable.
Most scientific and engineering problems are formulated into a symmetric system oflinear/non-linear equations, e. g., [ A]{X }={B} . Because of piecewise approximation,
the symmetric matrix [ A] may lead to a sparse figure, e. g.,[ ]°°° °° °°°°°°°°°°°°°°°

(1.1)

°°

where ° represents a non-zero fill-in, and other blanks represent zero coefficient. An
arbitrary symmetric sparse matrix, e. g., the above example, cannot take advantage ofsparse solvers. JCL renumbers the unknowns, e. g., rearrangement of rows and columns,
such that non-zero fill-ins are around the diagonal as:[ ]° ° °° ° °° ° ° °° ° °° ° °° °(1.2)

° °° °The sparse matrix, as shown in Eq. (1.2), then can be solved by constant bandwidth
solvers or variable bandwidth solvers. Variable bandwidth solver is also called skyline-1-solver.
The cost to solve [A]{X}={B} depends on bandwidth (or half bandwidth) andprofile. Half bandwidth is defined as the maximal ∣i− j∣ where i and j are subscripts
nof non-zero fill-in Aij or A ji , and profile is defined as ∑i=1 max of ∣i− j∣ where i andj are subscripts of non-zero fill-in A ji . Renumbering can significantly reduce solution

costs.
JCL provides an efficient and reliable way to reduce bandwidth and profile of asparse matrix. JCL is thread-safe, and can be called simultaneously in a program.

1.1 Non-zero fill-ins and Non-zero coefficientsNon-zero fill-in is different from non-zero coefficient. For the following example,[1.0 2.0
3.0 4.0 5.0[ A]=
5.0 6.0 7.0
7.0 8.0 9.0
9.0 1.0]

The matrix has 13 non-zero coefficients. The non-zero coefficients are asymmetricbecause ofA12≠ A21 . Non-zero fill-in is an entry with non-zero coefficient. The non-zero fill-ins of matrix [A] are represented as follows:[ ]° °° ° °° ° °° ° °° °In the above example, non-zero coefficients are asymmetric, but the non-zero fill-ins aresymmetric. Once non-zero fill-ins are symmetric, JCL is applicable.

-2-Chapter 2. Calling Syntax
JCL is programmed in Fortran. The entry subroutine to JCL is defined as follows:jcl(nas,n,adj,output,work,limit,flag)
where
1.

nas: 4-byte integer array, dimension(n), input the number of vertices adjacent
to each vertex. Concept of adjacent vertices comes from graph theory bywhich the renumbering algorithms are developed. In finite element, adjacentvertices are similar to connected degrees of freedom. In algebra, adjacentvertices are similar to coupled unknowns.

2.

n: a 4-byte integer, input the number of unknowns, or system size, or order.

3.

adj: 4-byte integer array, dimension(sum of values in array nas), input thevertices adjacent to each vertex. E. g., vertices adjacent to the first vertex, andthen vertices adjacent to the second vertex, and so on. The required length ofarray adj is equal to the sum of each value in array nas.

4. output: 4-byte integer array, dimension(n), returns the new number ofvertices.
5.

work: 4-byte integer array, dimension(limit), is a working array. This arrayshould have a sufficient space.

6.

limit: a 4-byte integer, input the dimension of array work. The dimension
cannot be determined in advance because of the nature of algorithm. Anestimation is as:
13*n+(sum of values in array nas)

JCL is a memory eater, but is not a CPU time eater. Computer memory is anot an issue on modern computer. The estimation can fit for most problems.
7.

flag: 4-byte integer, returns the status=1 , successful entry. The array output has new number ofvertices.
=0 , the limit of array work is insufficient. User can re-run JCL

with a big dimension of work.
-3-

= −i ,incorrect adjacent vertices, where i is an integer value.
User should check vertices adjacent to the vertex “i” andvertices adjacent to adjacent vertices of vertex “i”.

2.1 A Fortran Example
Let us consider the8×8 example as shown in Eq. (1.1). The information ofadjacent vertices can be summarized as:

VertexAdjacent Vertices

Number ofAdjacent Vertices

15123,6,83

32,52

47

1

51,32

62,8274

182,62

The following is a Fortran example showing how to call JCL.
integer (4), parameter :: n = 8
integer (4) :: nas(n)integer (4) :: output(n)integer (4) :: adj(14) !! why 14?
!! see program segment of!! “number of adjacent vertices”
integer (4), parameter :: limit = 13*n+14!! JCL is memory eater
integer (4) :: work(limit)
integer (4) :: flag!
!
!
!

number of adjacent vertices
by the 3rd column of above table-4-nas(1) = 1
nas(2) = 3
nas(3) = 2
nas(4) = 1
nas(5) = 2
nas(6) = 2
nas(7) = 1
nas(8) = 2!
!
!
!

!! sum of values = 14!! adj is declared as dimension(14)adjacent vertices
by the 2nd column of the above table!! vertex adjacent to vertex 1
adj(1) = 5!! vertices adjacent to vertex 2adj(2) = 3adj(3) = 6adj(4) = 8!! vertices adjacent to vertex 3adj(5) = 2adj(6) = 5!! vertex adjacent to vertex 4adj(7) = 7!! vertices adjacent to vertex 5adj(8) = 1adj(9) = 3!! vertices adjacent to vertex 6adj(10) = 2adj(11) = 8!! vertex adjacent to vertex 7adj(12) = 4!! vertices adjacent to vertex 8adj(13) = 2adj(14) = 6!
!
!

renumber by JCL!
!
!

output the new number

call jcl(nas,n,adj,output,work,limit,flag)

write(*,*) "flag",flagif(flag.eq.1) write(*,*) output

-5-!
!
!

end of programend

The output of new order is:

63

4852

71 . Thatmeans the vertex 1 should be numbered as 6, and vertex 2 should be numbered as 3, andso on. By the new order, user can copy rows of the original matrix in Eq. (1.1) to a newmatrix. For example, user can copy row 1 to row 6 of new matrix, and similarly copy row
2 to row 3, and so on. The row arrangement gets the following sparse figure:[ ]°°° °° °°°°°°°°°°°°°°°

(1.3)

°°

Then, by the new order (6 3 4 8 5 2 7 1), user can copy columns in Eq. (1.3) to a newmatrix. For example, copy column 1 to column 6 of new matrix, and similarly copy
column 2 to column 3, and so on. Arrangement of columns of the matrix in Eq. (1.3) getsthe perfect sparse figure as shown in Eq. (1.2).
The above example shows that bandwidth and profile of a sparse matrix are reduced
by rearranging rows and columns of matrix [A]. Renumbering after the matrix [A] isdetermined is less efficient and less popular. In most scientific and engineering
applications, vertices are renumbered before a sparse matrix is determined. The next
Chapter will illustrate how to renumber vertices without matrix [A].

2.2 A Bad Example
Let us use the previous

8×8 example to illustrate how to find error in input data.

-6-The information of adjacent vertices is as follows:

VertexAdjacent Vertices

Number ofAdjacent Vertices

151

2

3,6,83

3

2,5247151,32

62,8274

182,62The vertex 6 has adjacent vertices 2 and 8. Here on purpose input incorrect adjacentvertices 3 and 8. The example program is as follow:
integer (4), parameter :: n = 8
integer (4) :: nas(n)integer (4) :: output(n)integer (4) :: adj(14)integer (4), parameter :: limit = 13*n+14
integer (4) :: work(limit)
integer (4) :: flag!
!
!

number of adjacent vertices!
!
!

adjacent vertices

nas(1) = 1
nas(2) = 3
nas(3) = 2nas(4) = 1nas(5) = 2
nas(6) = 2nas(7) = 1nas(8) = 2!! vertex adjacent to vertex 1adj(1) = 5!! vertices adjacent to vertex 2-7-adj(2) = 3adj(3) = 6adj(4) = 8!! vertices adjacent to vertex 3adj(5) = 2adj(6) = 5!! vertex adjacent to vertex 4adj(7) = 7!! vertices adjacent to vertex 5adj(8) = 1adj(9) = 3!! vertices adjacent to vertex 6adj(10) = 3 !! (error: It should be 2)
adj(11) = 8!! vertex adjacent to vertex 7
adj(12) = 4!! vertices adjacent to vertex 8adj(13) = 2
adj(14) = 6!
!
!

renumber by JCL!
!
!

output the new number

call jcl(nas,n,adj,output,work,limit,flag)

write(*,*) "flag",flag
if(flag.eq.1) write(*,*) output!
!
!

end of programend

The output flag is -2. That gives a hint that an error may happen to the vertices adjacent
to vertex 2 or happens to the vertices adjacent to vertices adjacent to vertex 2. User needs
to check the following:
1) First, check vertices adjacent to vertex 2. Vertices adjacent to vertex 2 includes
3, 6, and 8. Data input to JCL is 3, 6, and 8. The input is correct. Then, checkvertices adjacent to 3, 6, and 8.
2) Vertices adjacent to vertex 3 are 2 and 5. Data input to JCL is 2 and 5. It iscorrect.
-8-3) Vertices adjacent to vertex 6 are 2 and 8. However, Data input to JCL is 3 and 8.
We find the error.
The return flag gives a hint to find an input error.

-9-Chapter 3. Renumbering Before the Determination ofSparse Matrix
This Chapter illustrates how to renumber vertices before the determination of matrix[A]. Renumbering is important in scientific and engineering computing. Most scientificand engineering problems can be formulated into a system of equations, i.e.,[ A]{X }={B} where matrix[ A] is symmetric and sparse. The unknowns can berenumbered before the determination of matrix[ A] . For example, let us consider theengineering problem, as shown in Figure 1.

It is a cable. The example solves horizontal and vertical displacements at five joints,
each of which has horizontal and vertical displacements. Each displacement is equivalent
to a vertex in renumbering algorithms. This example has a total of 10 vertices, and are
initially labeled from 1 to 10 in Figure 1. We won't use the initial number to solve
displacements, but use the initial number to define adjacent vertices. The adjacentvertices are required by JCL. By Figure 1, we can define the adjacent vertices as:
Vertex Adjacent Vertices

Number ofAdjacent Vertices

1

2,6,7

3

2

1,3,6,7,8

5

3

2,4,7,8,9

5- 10 -

Vertex Adjacent Vertices

Number ofAdjacent Vertices

4

3,5,8,9,10

5

5

4,9,10

3

6

1,2,7

3

7

1,2,3,6,8

5

8

2,3,4,7,95

9

3,4,5,8,10

5

10

4,5,9

3

A vertex is equivalent to a degree of freedom (hereinafter “dof”) in finite element
method, and adjacent vertices are equivalent to coupled degrees of freedom. If user has abackground with finite element method, it is easy to verify the coupled dof with respect toeach dof can be summarized as:
dofcoupled dof

1

2,6,7

2

1,3,6,7,8

3

2,4,7,8,9

4

3,5,8,9,10

5

4,9,10

6

1,2,7

7

1,2,3,6,8

8

2,3,4,7,9

9

3,4,5,8,10

10

4,5,9

Adjacent vertices and coupled dof identify connectivity. This manual uses the term“adjacent vertices” for illustration. Applying the above adjacent vertices to JCL, the newnumber of vertices is as:108642975

3

1E. g., replace “number 1” in Figure 1 with “number 10”, and replace “number 2” with- 11 -

“number 8”, and so on. The new number is labeled as:

We apply the new number labeled in Figure 2 to solve displacements, e.g., by finite
element method. The connectivity in Figure 2 shows the following:1) At vertex 1, vertices 1, 2, 3, and 4 are connected (or coupled). That means, on thefirst row of [A], columns 1, 2, 3, and 4 have non-zero fill-ins.
2) At vertex 2, vertices 1, 2, 3, and 4 are connected. On the second row of [A], columns
1, 2, 3, and 4 have non-zero fill-ins.
3) At vertex 3, vertices 1, 2, 3, 4, 5, and 6 are connected. One the third row of [A],columns 1, 2, 3, 4, 5, and 6 have non-zero fill-ins.
4) .... and so on.

[The new number of vertices leads the following matrix [ A] :A11 A12A21 A22A31 A32
A41 A42A13A23A33A43A53A63

A14A24A34A44A54A64

A35A45A55
A65
A75
A85

A36
A46
A56
A66
A76
A86

- 12 -A57
A67
A77
A87
A97
A10,7A58
A68
A78
A88A98
A10,8A79
A89
A99
A10,9A7,10A8,10A9,10A10,10]User can use the new numbers to calculate matrix [A] without a rearrangement of rows
and columns. The procedure to renumber vertices can be summarized as follows:
1) Arbitrarily number the vertices2) Define adjacent vertices by the initial number
3) Apply JCL to get new number.

- 13 -Chapter 4. Renumbering Asymmetric Sparse Matrix
When renumbering a matrix with asymmetric non-zero fill-ins, it is necessary to
expand the asymmetric figure into a symmetric figure. For example, the following has anasymmetric non-zero fill-ins:[ ]° °°°° ° °° ° °°°The asymmetric figure needs to be expanded into a symmetric figure by copying thetranspose of upper triangular part into the lower triangular part, and copying the transposeof lower triangular part into the upper triangular part. The expanded symmetric figure
becomes:[ ]° °° ° °° °°°°° °°° °° °Then, JCL renumbers the expanded sparse matrix that has a symmetric non-zero fill-ins.

- 14 -CHAPTER 5. Data Storage Schemes
Sparse matrix can be efficiently handled by a proper data storage schemes. Among
computer languages, Fortran is most suitable for sparse matrix operations. This chapter
will introduces data storage schemes for constant-bandwidth and variable-bandwidth
sparse symmetric matrices. For sparse asymmetric matrices, please refer to LAIPE
manual.

5.1 Constant Bandwidth Sparse Symmetric Matrix
The data storage scheme is applicable in subroutine (or function). Sparse matrix ispassed as a dummy argument. Since it is symmetric, we handle the lower triangular part.
For example,[14 25sym.
2 29 88
9 34 893 23 45
11 722329]

Data storage scheme for constant bandwidth sparse symmetric matrix is declared in a
Fortran program, for example,Integer :: mReal :: a(m,1)
where m is the half bandwidth and a is the sparse matrix. The total length of matrix a is(n-1)*m+n. Then, the coefficient

Aij is programmed in fortran as a(i,j). Fortranprovides a very convenient way to handle sparse matrix.

- 15 -

In the above example, m=2 and n=7, the total required length for matrix a is (71)*2+7=19. A subroutine to input the above example is as follow.
Subroutine input(a,m)
integer :: mreal :: a(m,1)
a(1,1) = 1.0a(2,1) = 4.0
a(3,1) = 2.0a(2,2) = 25.0a(3,2) = 29.0a(4,2) = 9.0
a(3,3) = 88.0
a(4,3) = 34.0a(5,3) = 3.0
a(4,4) = 89.0
a(5,4) = 23.0
a(6,4) = 11.0a(5,5) = 45.0a(6,5) = 7.0
a(7,5) = 3.0
a(6,6) = 22.0
a(7,6) = 2.0
a(7,7) = 9.0
returnend
Matrix a requires a length of 19. In the main program (or caller) needs to declare as:real :: a(19)and also need to define m=2.

5.2 Variable Bandwidth Sparse Symmetric Matrix
We handle the upper triangular part for variable bandwidth sparse symmetric matrix.
For example,- 16 -

[14272
3 6 234 022
5 22 13sym.
6 437]

Variable-bandwidth sparse symmetric solver is also called skyline solver. Date storage
scheme includes an address label, label, whose definition is as:
Set label (1) = 1
For i = 2 to n, do the followinglabel(i) = label(i-1)+[number of non-zerofill-ins above thediagonal in thei-th column]
The address labels for the above example are 1, 2, 3, 4, 7, 8, and 11. The data storage
scheme declares the sparse matrix in Fortran program, for example, asreal :: a(1,1)
Then, coefficient

Aij is programmed in Fortran program as A(i,label(j)). Fortranprovides a very convenient way to handle sparse matrix. The total required length ofmatrix a is defined as lable(n)-1+n where n is matrix order. In the above example, n=7
and label(7)=11, the required length is 11-1+7 = 17. It can be verified there are 17 nonzero fill-ins. A fortran subroutine to input the above example is as:
subroutine input (a,label)real :: a(1,1)
integer :: label(1)a(1,label(1)) = 1.0a(1,label(2)) = 4.0
a(2,label(2)) = 2.0
a(2,label(3)) = 7.0a(3,label(3)) = 3.0a(3,label(4)) = 6.0
a(4,label(4)) = 4.0a(3,label(5)) = 2.0- 17 -a(4,label(5)) = 23.0a(5,label(5)) = 0.0a(4,label(5)) = 5.0a(5,label(6)) = 22.0a(6,label(6)) = 6.0a(4,label(7)) = 22.0a(5,label(7)) = 13.0a(6,label(7)) = 43.0a(7,label(7)) = 7.0
returnend- 18 -