| 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+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
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: 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.
TABLE 5-3 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 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.
TABLE 5-4 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 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.
CODE EXAMPLE 5-1 Solving a Symmetric System-One-Call 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
|
CODE EXAMPLE 5-2 Solving a Symmetric System-Regular 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.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
|
CODE EXAMPLE 5-3 Solving a Structurally Symmetric System With Unsymmetric Values-Regular 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 5-4 Solving an Unsymmetric System-Regular 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. 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.
| Sun Performance Library User's Guide
|
819-0498-10
|
|
Copyright © 2004, Sun Microsystems, Inc. All Rights Reserved.