Sun Performance Library User's Guide


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+I1) = 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+I1) = 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 onetoone correspondence between the elements of the matrix and the elements of the array. Element A_{ij} 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 onetoone 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 1dimensional 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 userspecified 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 1based 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 HarwellBoeing 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. WH. 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 PrefixRootSuffix 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 51 lists the naming conventions for the Netlib Sparse BLAS vector routines.
TABLE 51 Netlib Sparse BLAS Naming Conventions
Operation

Root of Name

Prefix and Suffix

Dot product

DOT

SI DI CUI ZUI CCI ZCI

Scalar times a vector added to a vector

AXPY

SI DI CI ZI

Apply Givens rotation

ROT

SI DI

Gather x into y

GTHR

S D C Z SZ DZ CZ ZZ

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 sixcharacter name of the form XYYYZZ where:
 X represents the data type.
 YYY represents the sparse storage format.
 ZZ represents the operation.
TABLE 52 shows the values for X, Y, and Z.
TABLE 52 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
JAD: jagged diagonal
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: matrixmatrix 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 53.
TABLE 53 Sparse Solver Routines
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 userspecified ordering permutation.

DGSSRP()

Returns permutation used by solver.

DGSSCO()

Returns condition number estimate of coefficient matrix.

DGSSDA()

Deallocates 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() ! {fillreducing 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 onecall 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 righthand sides. Additional calls to dgsssl() are allowed to solve for additional righthand sides, as shown in the following example.
call dgssfs() ! {initialization, input coefficient matrix
! structure}
! {fillreducing 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 54.
TABLE 54 Sparse Solver Routine Calling Order
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 OneCall 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 51 shows solving a symmetric system using the onecall interface, and CODE EXAMPLE 52 shows solving a symmetric system using the regular interface.
CODE EXAMPLE 51 Solving a Symmetric SystemOneCall 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 onecall 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.528466159722D13
2 0.200000000000D+01 0.200000000000D+01 0.105249142734D12
3 0.100000000000D+01 0.100000000000D+01 0.350830475782D13
4 0.800000000000D+01 0.800000000000D+01 0.426325641456D13
5 0.500000000000D+00 0.500000000000D+00 0.660582699652D14

CODE EXAMPLE 52 Solving a Symmetric SystemRegular Interface
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.528466159722D13
2 0.200000000000D+01 0.200000000000D+01 0.105249142734D12
3 0.100000000000D+01 0.100000000000D+01 0.350830475782D13
4 0.800000000000D+01 0.800000000000D+01 0.426325641456D13
5 0.500000000000D+00 0.500000000000D+00 0.660582699652D14

CODE EXAMPLE 53 Solving a Structurally Symmetric System With Unsymmetric ValuesRegular Interface
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

CODE EXAMPLE 54 Solving an Unsymmetric SystemRegular Interface
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. WH. Liu. "Computer Solution of Large Sparse Positive Definite Systems." PrenticeHall 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:10341056, 1993.
 Ian S. Duff, Roger G. Grimes and John G. Lewis, "User's Guide for the HarwellBoeing Sparse Matrix Collection (Release I)," Technical Report TR/PA/92/86, CERFACS, Lyon, France, October 1992.
Sun Performance Library User's Guide

819049810


Copyright © 2004, Sun Microsystems, Inc. All Rights Reserved.