dpbco(l) — SunSoft Performance Library

NAME

dpbco - compute a Cholesky factorization and condition number of a symmetric positive definite matrix A in banded storage. If the condition number is not needed then xPBFA is slightly faster. It is typical to follow a call to xPBCO with a call to xPBSL to solve Ax = b or to xPBDI to compute the determinant of A.

SYNOPSIS

CALL DPBCO (DA, LDA, N, NDIAG, DRCOND, DWORK, INFO)

CALL SPBCO (SA, LDA, N, NDIAG, SRCOND, SWORK, INFO)

CALL ZPBCO (ZA, LDA, N, NDIAG, DRCOND, ZWORK, INFO)

CALL CPBCO (CA, LDA, N, NDIAG, SRCOND, CWORK, INFO)

ARGUMENTS

xAOn entry, the upper triangle of the matrix A.
On exit, a Cholesky factorization of the matrix A.

LDALeading dimension of the array A as specified in a dimension or
type statement. LDA >= NDIAG + 1.

NOrder of the matrix A. N ∗ 0.

NDIAGNumber of diagonals. N-1 >= NDIAG >= 0 but if N = 0 then NDIAG = 0.

xRCONDOn exit, an estimate of the reciprocal condition number of A.
0.0 <= RCOND <= 1.0. As the value of RCOND gets smaller, operations with A such as solving Ax = b may become less stable. If RCOND satisfies RCOND + 1.0 = 1.0 then A may be singular to working precision.

xWORKScratch array with a dimension of N.

INFOOn exit:
INFO = 0Subroutine completed normally.
INFO ∗ 0Returns a value k if the leading minor of order k is not positive definite.

SAMPLE PROGRAM

      PROGRAM TEST
      IMPLICIT NONE
C
      INTEGER           LDA, N, NDIAG
      PARAMETER        (N = 4)
      PARAMETER        (NDIAG = 1)
      PARAMETER        (LDA = NDIAG + 1)
C
      DOUBLE PRECISION A(LDA,N), B(N), RCOND, WORK(N)
      INTEGER           ICOL, INFO, IROW
C
      EXTERNAL          DPBCO, DPBSL
C
C     Initialize the array A to store in banded storage mode
C     the matrix A shown below. Initialize the array B to
C     store the vector B shown below.
C
C          2 -1   0   0        60
C     A = -1   2 -1   0    b = 60
C          0 -1   2 -1        60
C          0   0 -1   2        60
C
      DATA A / 8D8, 2.0D0, -1.0D0, 2.0D0, -1.0D0, 2.0D0, -1.0D0, 2.0D0 /
      DATA B / 4∗6.0D1 /
C
      PRINT 1000
      PRINT 1010, A(2,1), A(1,2)
      PRINT 1010, A(3,1), A(2,2), A(1,3)
      PRINT 1020,         A(3,2), A(2,3), A(1,4)
      PRINT 1030,                 A(3,3), A(2,4)
      PRINT 1040
      PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, LDA)
      PRINT 1050
      PRINT 1060, B
      CALL DPBCO (A, LDA, N, NDIAG, RCOND, WORK, INFO)
      IF (INFO .EQ. 0) THEN
        IF ((RCOND + 1.0D0) .EQ. 1.0D0) THEN
          PRINT 1100
        END IF
        CALL DPBSL (A, LDA, N, NDIAG, B)
        PRINT 1070, RCOND
        PRINT 1080
        PRINT 1060, B
      ELSE
        PRINT 1090
      END IF
C
1000 FORMAT (1X, ’A in full form:’)
1010 FORMAT (4(3X, F5.1))
1020 FORMAT (8X, 3(3X, F5.1))
1030 FORMAT (16X, 3(3X, F5.1))
1040 FORMAT (/1X, ’A in banded form: (∗ in unused entries)’)
1050 FORMAT (/1X, ’b:’)
1060 FORMAT (3X, F5.1)
1070 FORMAT (/1X, ’Reciprocal condition number of A:’, F5.1)
1080 FORMAT (/1X, ’A∗∗(-1) ∗ b:’)
1090 FORMAT (/1X, ’A is not positive definite.’)
1100 FORMAT (1X, ’A may be singular to working precision.’)
C
      END

SAMPLE OUTPUT

A in full form:
     2.0    -1.0
    -1.0     2.0    -1.0
            -1.0     2.0    -1.0
                    -1.0     2.0

A in banded form: (∗ in unused entries)
   ∗∗∗∗∗    -1.0    -1.0    -1.0
     2.0     2.0     2.0     2.0

b:
    60.0
    60.0
    60.0
    60.0

Reciprocal condition number of A: 0.1

A∗∗(-1) ∗ b:
   120.0
   180.0
   180.0
   120.0