Sparse Matrix Computations

Introduction

This section is concerned with methods for computing with sparse matrices. Our primary goal is to give the appearance of simplicity and allow ease-of-use in dealing with these calculations. The underlying principle in our design is to use Fortran 2003 standard support for derived types with initialized and allocatable components. To perform data storage and conversions we use overloaded assignment to hide complexity. The operations currently supported are:

defining entries of the matrices,
adding sparse matrices,
forming products of sparse matrices and dense vectors or matrices,
solving linear systems of algebraic equations
condition number computation
conversion of sparse matrices or dense arrays to the converse
storage management operations

The definition of the sparse matrices starts with a triplet consisting of the row and column indices and a value at that entry. By setting a flag in the derived type SLU_Options, repeated values may be accumulated to yield a value that is the sum of all triplets for that matrix entry. A diagram for constructing a single precision sparse 10000 × 10000 matrix, H, is illustrated with the pseudocode fragment:

Use linear_operators
Integer I, J; Real(Kind(1.e0)) value, x(10000)
Type(s_sparse) A
Type(s_hbc_sparse) H
                                                

Define non-zero values of A with repeated overloaded assignments
A = s_entry(I, J, value).

Convert to computational Harwell-Boeing form with the overloaded assignment H = A.

Compute with sparse matrix H, e. g., x = H .ix. x.

A basic feature is that there are four sparse matrix derived types, Types (s_hbc_sparse), (d_hbc_sparse), (c_hbc_sparse), and (z_hbc_sparse). These respectively handle single, double, complex and double-complex data. The defined operators work with a sparse matrix and a corresponding dense array of the same precision and data type. There is no mixing of data types such as a sparse double precision matrix multiplied by a single precision vector. To accommodate that case an intermediate double precision quantity will be created that ascends the single precision vector to a double precision vector. The table below shows the operations that are valid with sparse matrix types.

Mathematical Operation	Operation Notation	Input Terms	Output Terms
y = H-1x	y = H .ix. x	Hn×n sparse, x(1:k), k ≥ n	y(1:n)
y = xTH-1 ≡ H-Tx	y = x .xi. H	Hn×n sparse, x(1:k), k ≥ n	y(1:n)
Y = H-1Xn×r	Y= H .ix. X	Hn×n sparse, X(1:k,1:r), k ≥ n	Y(1:n,1:r)
Y = Xr×nH-1 ≡ (H-T XT)T	Y = X .xi. H	Hn×n sparse, X(1:r,1:k), k ≥ n	Y(1:r,1:n)
y = Hx	y = H .x. x	Hm×n sparse, x(1:k), k ≥ n	y(1:m)
y = xTH ≡ HTx	y = x .x. H	Hm×n sparse, x(1:k), k ≥ m	y(1:n)
Y = HXn×r	Y = H .x. X	Hm×n sparse, X(1:k,1:r), k ≥ n	Y(1:m,1:r)
Y = Xr×mH	Y = X .x. H	Hm×n sparse, X(1:r,1:k), k ≥ m	Y(1:r,1:n)
K = HT	K = .t. H	Hm×n sparse	Kn×m sparse
	K = .h. H	Hm×n sparse, complex	Kn×m sparse
y = HTx	y = H .tx. x	Hm×n sparse, x(1:k), k ≥ m	y(1:n)
Y = HTXm×r	Y = H .tx. X	Hm×n sparse, X(1:k,1:r), k ≥ m	Y(1:n,1:r)
y = xTH	Y = x .tx. H	Hm×n sparse, x(1:k), k ≥ m	y(1:n)
Y = XTr×mH	Y = X .tx. H	Hm×n sparse, X(1:k,1:r), k ≥ m	Y(1:r,1:n)
y = HxT	y = H .xt. x	Hm×n sparse, x(1:k), k ≥ n	y(1:m)
Y = HXTn×r	Y = H .xt. X	Hm×n sparse, x(1:k,1:r), k ≥ n	Y(1:m,1:r)
y = xHT	y = x .xt. H	Hm×n sparse, x(1:k), k ≥ n	y(1:m)
Y = Xr×nHT	Y = X .xt. H	Hm×n sparse, x(1:r,1:k), k ≥ n	Y(1:r,1:m)
	y = H .hx. x	Hm×n sparse1, x(1:k), k ≥ m	y(1:n)
	Y = H .hx. X	Hm×n sparse, X(1:k,1:r), k ≥ m	Y(1:n,1:r)
	Y = x .hx. H	Hm×n sparse, x(1:k), k ≥ m	y(1:n)
	Y = X .hx. H	Hm×n sparse, X(1:k,1:r), k ≥ m	Y(1:r,1:n)
	y = H .xh. x	Hm×n sparse, x(1:k), k ≥ n	y(1:m)
	Y = H .xh. X	Hm×n sparse, x(1:k,1:r), k ≥ n	Y(1:m,1:r)
	y = x .xh H	Hm×n sparse, x(1:k), k ≥ n	y(1:m)
	Y = X .xh. H	Hm×n sparse, x(1:r,1:k), k ≥ n	Y(1:r,1:m)

Note: 1 The operators .hx. and .xh. apply to sparse complex matrices only. For real matrices use the .tx. and .xt. operators.

Derived Type Definitions

A derived type is used for the entries (triplets or coordinate format) of a sparse matrix, which consists of row and column coordinates and a corresponding value:

    type s_entry
        integer irow
        integer jcol
        real(kind(1.e0)) value
    end type
                                                

Additionally, type (d_entry), type (c_entry), and type (z_entry) are defined similarly. These support double precision, complex and complex-double precision accuracy and types.

Thus for a sparse matrix A, the entry at the intersection of row irow and column jcol is the scalar value. We define a sparse matrix representation in terms of a collection of triplets. This is a convenient way for a user to define a sparse matrix. This representation is used to define the matrix entries in a user’s program using overloaded assignment. There is no implied order on the collection of triplets that define this sparse matrix. Our experience shows that for writing application code the technique of using triplets to define the matrix entries is convenient and provides a workable transition from mathematical definitions of the entries to computer code. Also note that there is generally no need for the programmer to allocate the components of a matrix of type s_sparse when using the overloaded assignment: s_sparse = s_entry. The software handles this detail by reallocating and expanding those components of the s_sparse matrix as required. (For this task we use the Fortran 2003 intrinsic subroutine move_alloc(), when it is available. This routine provides an efficient way to perform a reallocation.) The amount reallocated is controlled by an expansion factor that is a component of the derived type SLU_options.

     type s_sparse
         integer :: mrows = 0
         integer :: ncols = 0
         integer :: numnz = 0
         integer, allocatable, dimension(:) :: irow
         integer, allocatable, dimension(:) :: jcol
         real(kind(1.e0)), allocatable, dimension(:) :: value
         type (SLU_options) options
     end type
                                                

When performing matrix computations we use the Harwell-Boeing column-oriented derived type. The row indices, for each column, are unique and increasing. The values in the colptr(1:ncols) component mark the start of the row indices and corresponding matrix entries for that column. The value colptr(ncols+1)-1 will equal the value numnz after the matrix is defined with non-zero entries. The row indices for each column are in array irow(:). They are unique and sorted into increasing order.

     type s_hbc_sparse
         integer :: mrows = 0
         integer :: ncols = 0
         integer :: numnz = 0
         integer, allocatable, dimension(:) :: irow
         integer, allocatable, dimension(:) :: colptr
         real(kind(1.e0)), allocatable, dimension(:) :: value
         type(SLU_options) options
     end type
                                                

Additionally we support types (d_hbc_sparse), type (c_hbc_sparse), and type (z_hbc_sparse). These will have analogous support for the operations defined with type (s_hbc_sparse) and others that follow. From now on we only mention type (s_hbc_sparse).

All components of the type (s_sparse) object are self-explanatory except for the one named type(SLU_options). This component contains various parameters for managing the data structure, and for computing matrix products and linear system solutions. Normally these components do not need to be changed from their default values.

The derived type SLU_Options carries extra required information. That data needed for SuperLU2 is labeled with a comment. The remaining data is needed by IMSL codes that call on SuperLU. Of particular importance is the Sequence attribute statement. This prevents the Fortran compiler from rearranging the order of the components. Maintaining this order is required since the derived type SLU_Options is passed to a IMSL C code that uses the information as a C structure. The Sequence statement orders the Fortran-defined data so that it matches the C code structure.

Note: SuperLU is used to support the defined operations .ix. and .xi., and the condition number function, cond(). SuperLU is well-tested. Distributed and threaded versions are available but these are not used here in our software at present. SuperLU was developed by James W. Demmel, Stanley C. Eisenstat, John R. Gilbert, Xiaoye S. Li, and Joseph W. H. Liu. Note that the authors do not support the package in the context used in the IMSL Libraries.

Type SLU_options
         Sequence
         Integer :: unique = 1 ! Each new entry is unique –IMSL
         Integer :: Accumulate = 0
                            ! Accumulate or assemble duplicated entries in
                            ! a ?_sparse matrix.  This flag is checked
                            ! when executing an overloaded assignment
                            ! with a Harwell-Boeing = ?_sparse matrix.
                            ! The default is not to accumulate (0)
                            ! Assign the value 1 to accumulate.   
         Integer :: handle(2) = 0
                            ! Each HBC matrix requiring an LU
                            ! decomposition will have allocated
                            ! arrays whose start is pointed to by
                            ! this value.  In cases where the OS
                            ! uses 64 bit addressing 8 bytes are used.
         Integer :: Info = - 1
                            ! Flag returned after LU factorization (SuperLU) 
         Integer :: Fact = 0 !DOFACT - SuperLU
         Integer :: Equil = 1 !YES
         Integer :: ColPerm = 3 !COLAMD
         Integer :: Trans = 0 !NOTRANS
         Integer :: IterRefine = 1 !REFINE
         Integer :: PrintStat = 0  !NO
         Integer :: SymmetricMode = 0 !NO
         Integer :: PivotGrowth = 0 !NO
         Integer :: ConditionNumber = 0 !NO
         Integer :: RowPerm = 0 !NO
         Integer :: ReplaceTinyPivot = 0 !NO
         Integer :: SolveInitialized = 0 !NO
         Integer :: RefineInitialized = 0 !NO
         Real (Kind(1.d0)) :: DiagPivotThresh = 1.d0 ! SuperLU
         Real (Kind(1.d0)) :: expansion_factor = 1.2 ! VNI –
! The factor to use when expanding storage.  Any value > 1.
! can be used such that the integer part of this factor times
! any integer > 9 provides at least a value of 1 increase.
         Integer :: Cond_Iteration_Max = 30
! Maximum number of Lanczos and inverse iterations with sparse COND().
          Integer Alignment_Dummy
     End Type
                                                

Overloaded Assignments

A natural way to define a sparse matrix is in terms of its triplets. The basic tool used here to define all the non-zero entries is overloaded assignment. Fortran 90, and further updates to the standard, supports a hidden subroutine call, packaged in a module, when an assignment is executed between differing derived types. Thus if a Fortran program has a declaration type(s_sparse) A, then the overloaded assignment statement

A = s_entry(I, J, value)

has the effect of calling subroutines that result in joining the matrix entry value at the intersection of row I and column J. The components of A are managed to hold any number of values. The number of rows, columns and non-zero values are updated as new triplets are assigned. Also the arrays that hold the triplets are re-allocated and expanded, as required, to hold newly assigned triplets.

The code snippet for this operation, and others that follow, will require use of the module linear_operators. If new space is required in the assignment, a reallocation of the components of the matrix A will occur. The user does not have to manage the details.

Use linear_operators
Type(s_sparse) A
…
 
{For all entries in A, A = s_entry(I, J, Value)}
                                                

Sparse = Collection of Triplets

The Harwell-Boeing sparse matrix data types are used for computations. An assignment, H = A, implies deallocating any allocated components of H, allocating new storage, and sorting the collection of triplets provided as input in the sparse matrix A. If the accumulation flag is set in H%options%accumulate, the duplicate row indices in a column are reduced to a single entry and the corresponding values are added to yield a final value. The assignment H = 0 deallocates the allocated components and returns H to its initialized state, except for any changes to the component SLU_options. A similar comment holds for the assignment, A = 0.

Use linear_operators
Type(s_sparse) A
Type(s_hbc_sparse) H
…
{For all nonzero matrix entries, A = s_entry(I, J, Value)}
 
H = A
A = 0   ! Clear and deallocate components of A
…
H = 0   ! Clear and deallocate components of H
                                                

Sparse = Dense

The non-zero entries of the dense array are converted to a Harwell-Boeing sparse matrix. As a first step any allocated components are cleared and then allocated as needed to hold the non-zero values of the dense array. The specific dimensions of array D are arbitrary.

Use linear_operators
Type(s_hbc_sparse) H
Integer, parameter :: M=1000, N=1000
Real (kind(1.e0)) D(M,N) 
{Define entries of D}
H = D
                                                

Dense = Sparse

For some applications it is convenient to expand a sparse matrix into a dense matrix. The specific dimensions of array D are arbitrary.

Use linear_operators
Type(s_hbc_sparse) H
Integer, parameter :: M=1000, N=1000
Real (kind(1.e0)) D(M,N) 
{Define entries of H}
D = H
                                                

Scalar = s_hbc_entry(Sparse, I, J)

This assignment gets the value at the intersection of row I and column J of the Harwell-Boeing sparse matrix. There must be type agreement with the function and sparse matrix type. Use a prefix of d_, c_, or z_ for double, complex, or double complex values.

Use inear_operators
Type(s_hbc_sparse) H
Real (kind(1.e0)) value
{Define entries of H, I and J}
value = s_hbc_entry(H, I, J) 
                                                