Sun Performance Library User's Guide Table of ContentsPrevious ChapterNext ChapterBook Index


C H A P T E R  4
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:

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

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 

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 

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 

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 

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 

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 

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."

symmetric matrix 

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) not equal 0, then a(j, i) not equal 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.

struturally symmetric matrix 

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.

unsymmetric matix 

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 4-1 lists the naming conventions for the Netlib Sparse BLAS vector routines.

TABLE 4-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 4-2 shows the values for X, Y, and Z.

TABLE 4-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 4-3.

TABLE 4-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 4-4.

TABLE 4-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 4-1 shows solving a symmetric system using the one-call interface, and CODE EXAMPLE 4-2 shows solving a symmetric system using the regular interface.

CODE EXAMPLE 4-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 4-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 4-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 4-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 816-2463-10 Table of ContentsPrevious ChapterNext ChapterBook Index


Copyright © 2002, Sun Microsystems, Inc. All rights reserved.