C H A P T E R  5
 Working With Matrices

Most matrices can be stored in ways that save both storage space and computation time. Sun Performance Library uses the following storage schemes:

• Banded storage
• Packed storage

The Sun Performance Library processes matrices that are in one of four forms:

• General
• Triangular
• Symmetric
• Tridiagonal

Storage schemes and matrix types are described in the following sections.

## Matrix Storage Schemes

Some Sun Performance Library routines that work with arrays stored normally have corresponding routines that take advantage of these special storage forms. For example, DGBMV will form the product of a general matrix in banded storage and a vector, and DTPMV will form the product of a triangular matrix in packed storage and a vector.

### Banded Storage

A banded matrix is stored so the jth column of the matrix corresponds to the jth column of the Fortran array.

The following code copies a banded general matrix in a general array into banded storage mode.

 ``` C Copy the matrix A from the array AG to the array AB. The ``` ``` C matrix is stored in general storage mode in AG and it will ``` ``` C be stored in banded storage mode in AB. The code to copy ``` ``` C from general to banded storage mode is taken from the ``` ``` C comment block in the original DGBFA by Cleve Moler. ``` ``` C ``` ``` NSUB = 1 ``` ``` NSUPER = 2 ``` ``` NDIAG = NSUB + 1 + NSUPER ``` ``` DO ICOL = 1, N ``` ``` I1 = MAX0 (1, ICOL - NSUPER) ``` ``` I2 = MIN0 (N, ICOL + NSUB) ``` ``` DO IROW = I1, I2 ``` ``` IROWB = IROW - ICOL + NDIAG ``` ``` AB(IROWB,ICOL) = AG(IROW,ICOL) ``` ``` END DO ``` ``` END DO ```

Note that this method of storing banded matrices is compatible with the storage method used by LAPACK, BLAS, and LINPACK, but is inconsistent with the method used by EISPACK.

### Packed Storage

A packed vector is an alternate representation for a triangular, symmetric, or Hermitian matrix. An array is packed into a vector by storing the elements sequentially column by column into the vector. Space for the diagonal elements is always reserved, even if the values of the diagonal elements are known, such as in a unit diagonal matrix.

An upper triangular matrix or a symmetric matrix whose upper triangle is stored in general storage in the array A, can be transferred to packed storage in the array AP as shown below. This code comes from the comment block of the LAPACK routine DTPTRI.

 ``` JC = 1 ``` ``` DO J = 1, N ``` ``` DO I = 1, J ``` ``` AP(JC+I-1) = A(I,J) ``` ``` END DO ``` ``` JC = JC + J ``` ``` END DO ```

Similarly, a lower triangular matrix or a symmetric matrix whose lower triangle is stored in general storage in the array A, can be transferred to packed storage in the array AP as shown below:

 ``` JC = 1 ``` ``` DO J = 1, N ``` ``` DO I = J, N ``` ``` AP(JC+I-1) = A(I,J) ``` ``` END DO ``` ``` JC = JC + N - J + 1 ``` ``` END DO ```

## Matrix Types

The general matrix form is the most common matrix, and most operations performed by the Sun Performance Library can be done on general arrays. In many cases, there are routines that will work with the other forms of the arrays. For example, DGEMM will form the product of two general matrices and DTRMM will form the product of a triangular and a general matrix.

### General Matrices

A general matrix is stored so that there is a one-to-one correspondence between the elements of the matrix and the elements of the array. Element Aij of a matrix A is stored in element A(I,J) of the corresponding array A. The general form is the most common form. A general matrix, because it is dense, has no special storage scheme. In a general banded matrix, however, the diagonal of the matrix is stored in the row below the upper diagonals.

For example, as shown below, the general banded matrix can be represented with banded storage. Elements shown with the symbol × are never accessed by routines that process banded arrays.

 General Banded Matrix General Banded Array in Banded Storage

### Triangular Matrices

A triangular matrix is stored so that there is a one-to-one correspondence between the nonzero elements of the matrix and the elements of the array, but the elements of the array corresponding to the zero elements of the matrix are never accessed by routines that process triangular arrays.

A triangular matrix can be stored using packed storage.

 Triangular Matrix Triangular Array in Packed Storage

A triangular banded matrix can be stored using banded storage as shown below. Elements shown with the symbol × are never accessed by routines that process banded arrays.

 Triangular Banded Matrix Triangular Banded Array in Banded Storage

### Symmetric Matrices

A symmetric matrix is similar to a triangular matrix in that the data in either the upper or lower triangle corresponds to the elements of the array. The contents of the other elements in the array are assumed and those array elements are never accessed by routines that process symmetric or Hermitian arrays.

A symmetric matrix can be stored using packed storage.

 Symmetric Matrix Symmetric Array in Packed Storage

A symmetric banded matrix can be stored using banded storage as shown below. Elements shown with the symbol × are never accessed by routines that process banded arrays.

 Symmetric Banded Matrix Symmetric Banded Array in Banded Storage

### Tridiagonal Matrices

A tridiagonal matrix has elements only on the main diagonal, the first superdiagonal, and the first subdiagonal. It is stored using three 1-dimensional arrays.

 Tridiagonal Matrix Tridiagonal Array in Tridiagonal Storage

## Sparse Matrices

The Sun Performance Library sparse solver package is a collection of routines that efficiently factor and solve sparse linear systems of equations. Use the sparse solver package to:

• Solve symmetric, structurally symmetric, and unsymmetric coefficient matrices
• Specify a choice of ordering methods, including user-specified orderings

The sparse solver package contains interfaces for FORTRAN 77. Fortran 95 and C interfaces are not currently provided. To use the sparse solver routines from Fortran 95, use the FORTRAN 77 interfaces. To use the sparse solver routines with C, append an underscore to the routine name (dgssin_(), dgssor_(), and so on), pass arguments by reference, and use 1-based array indexing.

### Sparse Solver Matrix Data Formats

Sparse matrices are usually represented in formats that minimize storage requirements. By taking advantage of the sparsity and not storing zeros, considerable storage space can be saved. The storage format used by the general sparse solver is the compressed sparse column (CSC) format (also called the Harwell-Boeing format).

The CSC format represents a sparse matrix with two integer arrays and one floating point array. The integer arrays (colptr and rowind) specify the location of the nonzeros of the sparse matrix, and the floating point array (values) is used for the nonzero values.

The column pointer (colptr) array consists of n+1 elements where colptr(i) points to the beginning of the ith column, and colptr(i + 1) - 1 points to the end of the ith column. The row indices (rowind) array contains the row indices of the nonzero values. The values arrays contains the corresponding nonzero numerical values.

The following matrix data formats exist for a sparse matrix of neqns equations and nnz nonzeros:

• Symmetric
• Structurally symmetric
• Unsymmetric

The most efficient data representation often depends on the specific problem. The following sections show examples of sparse matrix data formats.

#### Symmetric Sparse Matrices

A symmetric sparse matrix is a matrix where a(i, j) = a(j, i) for all i and j. Because of this symmetry, only the lower triangular values need to be passed to the solver routines. The upper triangle can be determined from the lower triangle.

An example of a symmetric matrix is shown below. This example is derived from A. George and J. W-H. Liu. "Computer Solution of Large Sparse Positive Definite Systems."

To represent A in CSC format:

• colptr: 1, 6, 7, 8, 9, 10
• rowind: 1, 2, 3, 4, 5, 2, 3, 4, 5
• values: 4.0, 1.0, 2.0, 0.5, 2.0, 0.5, 3.0, 0.625, 16.0

#### Structurally Symmetric Sparse Matrices

A structurally symmetric sparse matrix has nonzero values with the property that if a(i, j) 0, then a(j, i) 0 for all i and j. When solving a structurally symmetric system, the entire matrix must be passed to the solver routines.

An example of a structurally symmetric matrix is shown below.

To represent A in CSC format:

• colptr: 1, 3, 6, 7, 9
• rowind: 1, 2, 1, 2, 4, 3, 2, 4
• values: 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0

#### Unsymmetric Sparse Matrices

An unsymmetric sparse matrix does not have a(i, j) = a(j, i) for all i and j. The structure of the matrix does not have an apparent pattern. When solving an unsymmetric system, the entire matrix must be passed to the solver routines. An example of an unsymmetric matrix is shown below.

To represent A in CSC format:

• colptr: 1, 6, 7, 8, 9, 11
• rowind: 1, 2, 3, 4, 5, 2, 3, 4, 2, 5
• values: 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0

### Sun Performance Library Sparse BLAS

The Sun Performance Library sparse BLAS package is based on the following two packages:

• Netlib Sparse BLAS package, by Dodson, Grimes, and Lewis consists of sparse extensions to the Basic Linear Algebra Subroutines that operate on sparse vectors.
• NIST (National Institute of Standards and Technology) Fortran Sparse BLAS Library consists of routines that perform matrix products and solution of triangular systems for sparse matrices in a variety of storage formats.

Refer to the following sources for additional sparse BLAS information.

### Naming Conventions

The Netlib Sparse BLAS and NIST Fortran Sparse BLAS Library routines each use their own naming conventions, as described in the following two sections.

#### Netlib Sparse BLAS

Each Netlib Sparse BLAS routine has a name of the form Prefix-Root-Suffix where the:

• Prefix represents the data type.
• Root represents the operation.
• Suffix represents whether or not the routine is a direct extension of an existing dense BLAS routine.

TABLE 5-1 lists the naming conventions for the Netlib Sparse BLAS vector routines.

TABLE 5-1 Netlib Sparse BLAS Naming Conventions

Operation

Root of Name

Prefix and Suffix

Dot product

-DOT-

S-I D-I C-UI Z-UI C-CI Z-CI

Scalar times a vector added to a vector

-AXPY-

S-I D-I C-I Z-I

Apply Givens rotation

-ROT-

S-I D-I

Gather x into y

-GTHR-

S- D- C- Z- S-Z D-Z C-Z Z-Z

Scatter x into y

-SCTR-

S- D- C- Z-

The prefix can be one of the following data types:

• S: SINGLE
• D: DOUBLE
• C: COMPLEX
• Z: COMPLEX*16 or DOUBLE COMPLEX

The I, CI, and UI suffixes denote sparse BLAS routines that are direct extensions to dense BLAS routines.

#### NIST Fortran Sparse BLAS

Each NIST Fortran Sparse BLAS routine has a six-character name of the form XYYYZZ where:

• X represents the data type.
• YYY represents the sparse storage format.
• ZZ represents the operation.

TABLE 5-2 shows the values for X, Y, and Z.

TABLE 5-2 NIST Fortran Sparse BLAS Routine Naming Conventions

X: Data Type

X

S: single precision

D: double precision

C: complex

Z: double complex

YYY: Sparse Storage Format

YYY

Single entry formats:

COO: coordinate

CSC: compressed sparse column

CSR: compressed sparse row

DIA: diagonal

ELL: ellpack

SKY: skyline

Block entry formats:

BCO: block coordinate

BSC: block compressed sparse column

BSR: block compressed sparse row

BDI: block diagonal

BEL: block ellpack

VBR: block compressed sparse row

ZZ: Operation

ZZ

MM: matrix-matrix product

SM: solution of triangular system (supported for all formats except COO)

RP: right permutation (for JAD format only)

### Sparse Solver Routines

The Sun Performance Library sparse solver package contains the routines listed in TABLE 5-3.

 Routine Function DGSSFS() One call interface to sparse solver DGSSIN() Sparse solver initialization DGSSOR() Fill reducing ordering and symbolic factorization DGSSFA() Matrix value input and numeric factorization DGSSSL() Triangular solve Utility Routine Function DGSSUO() Sets user-specified ordering permutation. DGSSRP() Returns permutation used by solver. DGSSCO() Returns condition number estimate of coefficient matrix. DGSSDA() De-allocates sparse solver. DGSSPS() Prints solver statistics.

Use the regular interface to solve multiple matrices with the same structure, but different numerical values, as shown below:

 ```call dgssin() ! {initialization, input coefficient matrix ``` ``` ! structure} ``` ```call dgssor() ! {fill-reducing ordering, symbolic factorization} ``` ```do m = 1, number_of_structurally_identical_matrices ``` ``` call dgssfa() ! {input coefficient matrix values, numeric ``` ``` ! factorization} ``` ``` do r = 1, number_of_right_hand_sides ``` ``` call dgsssl() ! {triangular solve} ``` ``` enddo ``` ```enddo ```

The one-call interface is not as flexible as the regular interface, but it covers the most common case of factoring a single matrix and solving some number right-hand sides. Additional calls to dgsssl() are allowed to solve for additional right-hand sides, as shown in the following example.

 ```call dgssfs() ! {initialization, input coefficient matrix ``` ``` ! structure} ``` ``` ! {fill-reducing ordering, symbolic factorization} ``` ``` ! {input coefficient matrix values, numeric ``` ``` ! factorization} ``` ``` ! {triangular solve} ``` ```do r = 1, number_of_right_hand_sides ``` ``` call dgsssl() ! {triangular solve} ``` ```enddo ```

### Routine Calling Order

To solve problems with the sparse solver package, use the sparse solver routines in the order shown in TABLE 5-4.

 One Call Interface: For solving single matrix Start DGSSFS() Initialize, order, factor, solve DGSSSL() Additional solves (optional): repeat dgsssl() as needed DGSSDA() Deallocate working storage Finish End of One-Call Interface Regular Interface: For solving multiple matrices with the same structure Start DGSSIN() Initialize DGSSOR() Order DGSSFA() Factor DGSSSL() Solve: repeat dgssfa() or dgsssl() as needed DGSSDA() Deallocate working storage Finish End of Regular Interface

### Sparse Solver Examples

CODE EXAMPLE 5-1 shows solving a symmetric system using the one-call interface, and CODE EXAMPLE 5-2 shows solving a symmetric system using the regular interface.

 ``` ``` ```my_system% cat example_1call.f ``` ``` program example_1call ``` ```c ``` ```c This program is an example driver that calls the sparse solver. ``` ```c It factors and solves a symmetric system, by calling the ``` ```c one-call interface. ``` ```c ``` ``` implicit none ``` ```  ``` ``` integer neqns, ier, msglvl, outunt, ldrhs, nrhs ``` ``` character mtxtyp*2, pivot*1, ordmthd*3 ``` ``` double precision handle(150) ``` ``` integer colstr(6), rowind(9) ``` ``` double precision values(9), rhs(5), xexpct(5) ``` ``` integer i ``` ```c ``` ```c Sparse matrix structure and value arrays. From George and Liu, ``` ```c page 3. ``` ```c Ax = b, (solve for x) where: ``` ```c ``` ```c 4.0 1.0 2.0 0.5 2.0 2.0 7.0 ``` ```c 1.0 0.5 0.0 0.0 0.0 2.0 3.0 ``` ```c A = 2.0 0.0 3.0 0.0 0.0 x = 1.0 b = 7.0 ``` ```c 0.5 0.0 0.0 0.625 0.0 -8.0 -4.0 ``` ```c 2.0 0.0 0.0 0.0 16.0 -0.5 -4.0 ``` ```c ``` ``` data colstr / 1, 6, 7, 8, 9, 10 / ``` ``` data rowind / 1, 2, 3, 4, 5, 2, 3, 4, 5 / ``` ``` data values / 4.0d0, 1.0d0, 2.0d0, 0.5d0, 2.0d0, 0.5d0, 3.0d0, ``` ``` & 0.625d0, 16.0d0 / ``` ``` data rhs / 7.0d0, 3.0d0, 7.0d0, -4.0d0, -4.0d0 / ``` ``` data xexpct / 2.0d0, 2.0d0, 1.0d0, -8.0d0, -0.5d0 / ``` ```c ``` ```c set calling parameters ``` ```c ``` ``` mtxtyp= 'ss' ``` ``` pivot = 'n' ``` ``` neqns = 5 ``` ``` nrhs = 1 ``` ```  ``` ``` ldrhs = 5 ``` ``` outunt = 6 ``` ``` msglvl = 0 ``` ``` ordmthd = 'mmd' ``` ```c ``` ```c call single call interface ``` ```c ``` ``` call dgssfs ( mtxtyp, pivot, neqns , colstr, rowind, ``` ``` & values, nrhs , rhs, ldrhs , ordmthd, ``` ``` & outunt, msglvl, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c deallocate sparse solver storage ``` ```c ``` ``` call dgssda ( handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c print values of sol ``` ```c ``` ``` write(6,200) 'i', 'rhs(i)', 'expected rhs(i)', 'error' ``` ``` do i = 1, neqns ``` ``` write(6,300) i, rhs(i), xexpct(i), (rhs(i)-xexpct(i)) ``` ``` enddo ``` ``` stop ``` ``` 110 continue ``` ```c ``` ```c call to sparse solver returns an error ``` ```c ``` ``` write ( 6 , 400 ) ``` ``` & ' example: FAILED sparse solver error number = ', ier ``` ``` stop ``` ```  ``` ``` 200 format(a5,3a20) ``` ```  ``` ``` 300 format(i5,3d20.12) ! i/sol/xexpct values ``` ```  ``` ``` 400 format(a60,i20) ! fail message, sparse solver error number ``` ```  ``` ``` end ``` ```  ``` ```my_system% f95 -dalign example_1call.f -xlic_lib=sunperf ``` ```my_sytem% a.out ``` ``` i rhs(i) expected rhs(i) error ``` ``` 1 0.200000000000D+01 0.200000000000D+01 -0.528466159722D-13 ``` ``` 2 0.200000000000D+01 0.200000000000D+01 0.105249142734D-12 ``` ``` 3 0.100000000000D+01 0.100000000000D+01 0.350830475782D-13 ``` ``` 4 -0.800000000000D+01 -0.800000000000D+01 0.426325641456D-13 ``` ``` 5 -0.500000000000D+00 -0.500000000000D+00 0.660582699652D-14 ``` ```  ```

 ``` ``` ```my_system% cat example_ss.f ``` ``` program example_ss ``` ```c ``` ```c This program is an example driver that calls the sparse solver. ``` ```c It factors and solves a symmetric system. ``` ```  ``` ``` implicit none ``` ```  ``` ``` integer neqns, ier, msglvl, outunt, ldrhs, nrhs ``` ``` character mtxtyp*2, pivot*1, ordmthd*3 ``` ``` double precision handle(150) ``` ``` integer colstr(6), rowind(9) ``` ``` double precision values(9), rhs(5), xexpct(5) ``` ``` integer i ``` ```c ``` ```c Sparse matrix structure and value arrays. From George and Liu, ``` ```c page 3. ``` ```c Ax = b, (solve for x) where: ``` ```c ``` ```c 4.0 1.0 2.0 0.5 2.0 2.0 7.0 ``` ```c 1.0 0.5 0.0 0.0 0.0 2.0 3.0 ``` ```c A = 2.0 0.0 3.0 0.0 0.0 x = 1.0 b = 7.0 ``` ```c 0.5 0.0 0.0 0.625 0.0 -8.0 -4.0 ``` ```c 2.0 0.0 0.0 0.0 16.0 -0.5 -4.0 ``` ```c ``` ``` data colstr / 1, 6, 7, 8, 9, 10 / ``` ``` data rowind / 1, 2, 3, 4, 5, 2, 3, 4, 5 / ``` ``` data values / 4.0d0, 1.0d0, 2.0d0, 0.5d0, 2.0d0, 0.5d0, ``` ``` & 3.0d0, 0.625d0, 16.0d0 / ``` ``` data rhs / 7.0d0, 3.0d0, 7.0d0, -4.0d0, -4.0d0 / ``` ``` data xexpct / 2.0d0, 2.0d0, 1.0d0, -8.0d0, -0.5d0 / ``` ```  ``` ```c ``` ```c initialize solver ``` ```c ``` ``` mtxtyp= 'ss' ``` ``` pivot = 'n' ``` ``` neqns = 5 ``` ``` outunt = 6 ``` ``` msglvl = 0 ``` ```c ``` ```c call regular interface ``` ```c ``` ``` call dgssin ( mtxtyp, pivot, neqns , colstr, rowind, ``` ``` & outunt, msglvl, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c ordering and symbolic factorization ``` ```c ``` ``` ordmthd = 'mmd' ``` ``` call dgssor ( ordmthd, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c numeric factorization ``` ```c ``` ``` call dgssfa ( neqns, colstr, rowind, values, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c solution ``` ```c ``` ``` nrhs = 1 ``` ``` ldrhs = 5 ``` ``` call dgsssl ( nrhs, rhs, ldrhs, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c deallocate sparse solver storage ``` ```c ``` ``` call dgssda ( handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c print values of sol ``` ```c ``` ``` write(6,200) 'i', 'rhs(i)', 'expected rhs(i)', 'error' ``` ``` do i = 1, neqns ``` ``` write(6,300) i, rhs(i), xexpct(i), (rhs(i)-xexpct(i)) ``` ``` enddo ``` ``` stop ``` ```  ``` ``` 110 continue ``` ```c ``` ```c call to sparse solver returns an error ``` ```c ``` ``` write ( 6 , 400 ) ``` ``` & ' example: FAILED sparse solver error number = ', ier ``` ``` stop ``` ```  ``` ``` 200 format(a5,3a20) ``` ```  ``` ``` 300 format(i5,3d20.12) ! i/sol/xexpct values ``` ```  ``` ``` 400 format(a60,i20) ! fail message, sparse solver error number ``` ``` ``` ``` end ``` ```my_system% f95 -dalign example_ss.f -xlic_lib=sunperf ``` ```my_sytem% a.out ``` ``` i rhs(i) expected rhs(i) error ``` ``` 1 0.200000000000D+01 0.200000000000D+01 -0.528466159722D-13 ``` ``` 2 0.200000000000D+01 0.200000000000D+01 0.105249142734D-12 ``` ``` 3 0.100000000000D+01 0.100000000000D+01 0.350830475782D-13 ``` ``` 4 -0.800000000000D+01 -0.800000000000D+01 0.426325641456D-13 ``` ``` 5 -0.500000000000D+00 -0.500000000000D+00 0.660582699652D-14 ``` ```  ```

 ``` ``` ```my_system% cat example_su.f ``` ``` program example_su ``` ```c ``` ```c This program is an example driver that calls the sparse solver. ``` ```c It factors and solves a structurally symmetric system ``` ```c (w/unsymmetric values). ``` ```c ``` ``` implicit none ``` ```  ``` ``` integer neqns, ier, msglvl, outunt, ldrhs, nrhs ``` ``` character mtxtyp*2, pivot*1, ordmthd*3 ``` ``` double precision handle(150) ``` ``` integer colstr(5), rowind(8) ``` ``` double precision values(8), rhs(4), xexpct(4) ``` ``` integer i ``` ```  ``` ```c ``` ```c Sparse matrix structure and value arrays. Coefficient matrix ``` ```c has a symmetric structure and unsymmetric values. ``` ```c Ax = b, (solve for x) where: ``` ```c ``` ```c 1.0 3.0 0.0 0.0 1.0 7.0 ``` ```c 2.0 4.0 0.0 7.0 2.0 38.0 ``` ```c A = 0.0 0.0 6.0 0.0 x = 3.0 b = 18.0 ``` ```c 0.0 5.0 0.0 8.0 4.0 42.0 ``` ```c ``` ``` data colstr / 1, 3, 6, 7, 9 / ``` ``` data rowind / 1, 2, 1, 2, 4, 3, 2, 4 / ``` ``` data values / 1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0, 7.0d0, ``` ``` & 8.0d0 / ``` ``` data rhs / 7.0d0, 38.0d0, 18.0d0, 42.0d0 / ``` ``` data xexpct / 1.0d0, 2.0d0, 3.0d0, 4.0d0 / ``` ```c ``` ```c initialize solver ``` ```c ``` ``` mtxtyp= 'su' ``` ``` pivot = 'n' ``` ``` neqns = 4 ``` ``` outunt = 6 ``` ``` msglvl = 0 ``` ```c ``` ```c call regular interface ``` ```c ``` ``` call dgssin ( mtxtyp, pivot, neqns , colstr, rowind, ``` ``` & outunt, msglvl, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c ordering and symbolic factorization ``` ```c ``` ``` ordmthd = 'mmd' ``` ``` call dgssor ( ordmthd, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c numeric factorization ``` ```c ``` ``` call dgssfa ( neqns, colstr, rowind, values, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```  ``` ```c ``` ```c solution ``` ```c ``` ``` nrhs = 1 ``` ``` ldrhs = 4 ``` ``` call dgsssl ( nrhs, rhs, ldrhs, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c deallocate sparse solver storage ``` ```c ``` ``` call dgssda ( handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c print values of sol ``` ```c ``` ``` write(6,200) 'i', 'rhs(i)', 'expected rhs(i)', 'error' ``` ``` do i = 1, neqns ``` ``` write(6,300) i, rhs(i), xexpct(i), (rhs(i)-xexpct(i)) ``` ``` enddo ``` ``` stop ``` ``` 110 continue ``` ```c ``` ```c call to sparse solver returns an error ``` ```c ``` ``` write ( 6 , 400 ) ``` ``` & ' example: FAILED sparse solver error number = ', ier ``` ``` stop ``` ```  ``` ``` 200 format(a5,3a20) ``` ```  ``` ``` 300 format(i5,3d20.12) ! i/sol/xexpct values ``` ```  ``` ``` 400 format(a60,i20) ! fail message, sparse solver error number ``` ```  ``` ``` end ``` ```my_system% f95 -dalign example_su.f -xlic_lib=sunperf ``` ```my_system% a.out ``` ``` i rhs(i) expected rhs(i) error ``` ``` 1 0.100000000000D+01 0.100000000000D+01 0.000000000000D+00 ``` ``` 2 0.200000000000D+01 0.200000000000D+01 0.000000000000D+00 ``` ``` 3 0.300000000000D+01 0.300000000000D+01 0.000000000000D+00 ``` ``` 4 0.400000000000D+01 0.400000000000D+01 0.000000000000D+00 ``` ```  ```

 ``` ``` ```my_system% cat example_uu.f ``` ``` program example_uu ``` ```c ``` ```c This program is an example driver that calls the sparse solver. ``` ```c It factors and solves an unsymmetric system. ``` ```c ``` ``` implicit none ``` ```  ``` ``` integer neqns, ier, msglvl, outunt, ldrhs, nrhs ``` ``` character mtxtyp*2, pivot*1, ordmthd*3 ``` ``` double precision handle(150) ``` ``` integer colstr(6), rowind(10) ``` ``` double precision values(10), rhs(5), xexpct(5) ``` ``` integer i ``` ```c ``` ```c Sparse matrix structure and value arrays. Unsummetric matrix A. ``` ```c Ax = b, (solve for x) where: ``` ```c ``` ```c 1.0 0.0 0.0 0.0 0.0 1.0 1.0 ``` ```c 2.0 6.0 0.0 0.0 9.0 2.0 59.0 ``` ```c A = 3.0 0.0 7.0 0.0 0.0 x = 3.0 b = 24.0 ``` ```c 4.0 0.0 0.0 8.0 0.0 4.0 36.0 ``` ```c 5.0 0.0 0.0 0.0 10.0 5.0 55.0 ``` ```c ``` ``` data colstr / 1, 6, 7, 8, 9, 11 / ``` ``` data rowind / 1, 2, 3, 4, 5, 2, 3, 4, 2, 5 / ``` ``` data values / 1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0, 7.0d0, ``` ``` & 8.0d0, 9.0d0, 10.0d0 / ``` ``` data rhs / 1.0d0, 59.0d0, 24.0d0, 36.0d0, 55.0d0 / ``` ``` data xexpct / 1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0 / ``` ```c ``` ```c initialize solver ``` ```c ``` ``` mtxtyp= 'uu' ``` ``` pivot = 'n' ``` ``` neqns = 5 ``` ``` outunt = 6 ``` ``` msglvl = 3 ``` ``` call dgssin ( mtxtyp, pivot, neqns , colstr, rowind, ``` ``` & outunt, msglvl, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```  ``` ```c ``` ```c ordering and symbolic factorization ``` ```c ``` ``` ordmthd = 'mmd' ``` ``` call dgssor ( ordmthd, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c numeric factorization ``` ```c ``` ``` call dgssfa ( neqns, colstr, rowind, values, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c solution ``` ```c ``` ``` nrhs = 1 ``` ``` ldrhs = 5 ``` ``` call dgsssl ( nrhs, rhs, ldrhs, handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c deallocate sparse solver storage ``` ```c ``` ``` call dgssda ( handle, ier ) ``` ``` if ( ier .ne. 0 ) goto 110 ``` ```c ``` ```c print values of sol ``` ```c ``` ``` write(6,200) 'i', 'rhs(i)', 'expected rhs(i)', 'error' ``` ``` do i = 1, neqns ``` ``` write(6,300) i, rhs(i), xexpct(i), (rhs(i)-xexpct(i)) ``` ``` enddo ``` ``` stop ``` ``` 110 continue ``` ```c ``` ```c call to sparse solver returns an error ``` ```c ``` ``` write ( 6 , 400 ) ``` ``` & ' example: FAILED sparse solver error number = ', ier ``` ``` stop ``` ```  ``` ``` 200 format(a5,3a20) ``` ```  ``` ``` 300 format(i5,3d20.12) ! i/sol/xexpct values ``` ```  ``` ``` 400 format(a60,i20) ! fail message, sparse solver error number ``` ``` end ``` ```  ``` ```my_system% f95 -dalign example_uu.f -xlic_lib=sunperf ``` ```my_system% a.out ``` ``` i rhs(i) expected rhs(i) error ``` ``` 1 0.100000000000D+01 0.100000000000D+01 0.000000000000D+00 ``` ``` 2 0.200000000000D+01 0.200000000000D+01 0.000000000000D+00 ``` ``` 3 0.300000000000D+01 0.300000000000D+01 0.000000000000D+00 ``` ``` 4 0.400000000000D+01 0.400000000000D+01 0.000000000000D+00 ``` ``` 5 0.500000000000D+01 0.500000000000D+01 0.000000000000D+00 ``` ```  ```

## References

The following books and papers provide additional information for the sparse BLAS and sparse solver routines.

• Dodson, D.S, R.G. Grimes, and J.G. Lewis. "Sparse Extensions to the Fortran Basic Linear Algebra Subprograms." ACM Transactions on Mathematical Software, June 1991, Vol 17, No. 2.
• A. George and J. W-H. Liu. "Computer Solution of Large Sparse Positive Definite Systems." Prentice-Hall Inc., Englewood Cliffs, New Jersey, 1981.
• E. Ng and B. W. Peyton. "Block Sparse Cholesky Algorithms on Advanced Uniprocessor Computers." SIAM M. Sci Comput., 14:1034-1056, 1993.
• Ian S. Duff, Roger G. Grimes and John G. Lewis, "User's Guide for the Harwell-Boeing Sparse Matrix Collection (Release I)," Technical Report TR/PA/92/86, CERFACS, Lyon, France, October 1992.