Previous Next Contents Generated Index Doc Set Home



Reciprocal Condition Number of a Cholesky-Factored Symmetric Positive Definite Matrix in Banded Storage

The subroutines described in this section estimate the reciprocal condition number of a real symmetric (or Hermitian) positive definite matrix A in banded storage, which has been Cholesky-factored by xPBTRF.

Calling Sequence

CALL DPBCON 
(UPLO, N, NDIAG, DA, LDA, DANORM, DRCOND, DWORK, 
IWORK2, INFO)
CALL SPBCON 
(UPLO, N, NDIAG, SA, LDA, SANORM, SRCOND, SWORK, 
IWORK2, INFO)
CALL ZPBCON 
(UPLO, N, NDIAG, ZA, LDA, DANORM, DRCOND, ZWORK, 
DWORK2, INFO)
CALL CPBCON 
(UPLO, N, NDIAG, CA, LDA, SANORM, SRCOND, CWORK, 
SWORK2, INFO)






void dpbcon 
(char uplo, int n, int ndiag, double *da, int lda, 
double danorm, double *drcond, int *info)
void spbcon 
(char uplo, int n, int ndiag, float *sa, int lda, float 
sanorm, float *srcond, int *info)
void zpbcon 
(char uplo, int n, int ndiag, doublecomplex *za, int 
lda, double danorm, double *drcond, int *info)
void cpbcon 
(char uplo, int n, int ndiag, complex *ca, int lda, 
float sanorm, float *srcond, int *info)

Arguments

UPLO

Indicates whether xA contains the upper or lower triangle of the matrix. The legal values for UPLO are listed below. Any values not listed below are illegal.

'U' or 'u'

xA contains the upper triangle.

'L' or 'l'

xA contains the lower triangle.

N

Order of the matrix A. N 0.

NDIAG

Number of superdiagonals or subdiagonals of the matrix A. N-1 NDIAG 0 but if N = 0 then NDIAG = 0.

If UPLO = 'U' or 'u', NDIAG is the number of superdiagonals.

If UPLO = 'L' or 'l', NDIAG is the number of subdiagonals.

xA

Cholesky factorization of the matrix A, as computed by xPBTRF.

LDA

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

xANORM

The 1-norm or -norm of matrix A.

xRCOND

On exit, the estimated reciprocal of the condition number of the matrix A where the reciprocal condition number of A is defined to be 1 / (ANORM × ||A-1||). The reciprocal of the condition number is estimated instead of the condition number itself to avoid overflow or division by zero. If RCOND is less than machine precision (in particular, if RCOND = 0) then A is singular to working precision.

xWORK

Scratch array with a dimension of 3 × N for real subroutines or 2 × N for complex subroutines.

xWORK2

Scratch array with a dimension of N.

INFO

On exit:

INFO = 0

Subroutine completed normally.

INFO < 0

The ith argument, where i = |INFO|, had an illegal value.

Sample Program




      PROGRAM TEST
      IMPLICIT NONE
C
      INTEGER           LDA, LDWORK, N, NDIAG
      PARAMETER        (N = 4)
      PARAMETER        (NDIAG = 1)
      PARAMETER        (LDA = NDIAG + 1)
      PARAMETER        (LDWORK = 3 * N)
C
      DOUBLE PRECISION  A(LDA,N), ANORM, RCOND, WORK(LDWORK)
      INTEGER           ICOL, INFO, IROW, IWORK(N)
C
      EXTERNAL          DPBCON, DPBTRF
      INTRINSIC         ABS
C
C     Initialize the array A to store in symmetric banded form
C     the 4x4 symmetric positive definite matrix A with one
C     subdiagonal and one superdiagonal shown below.
C
C          2  -1   0   0
C     A = -1   2  -1   0
C          0  -1   2  -1
C          0   0  -1   2
C
      DATA A /    8D8, 2.0D0, -1.0D0, 2.0D0,
     $         -1.0D0, 2.0D0, -1.0D0, 2.0D0 /
C
C     Print the initial values of the arrays.
C
      PRINT 1000
      PRINT 1010, A(2,1), A(1,2),  0.0D0, 0.0D0
      PRINT 1010, A(1,2), A(2,2), A(1,3), 0.0D0
      PRINT 1010,  0.0D0, A(1,3), A(2,3), A(1,4)
      PRINT 1010,  0.0D0,  0.0D0, A(1,4), A(2,4)
      PRINT 1020
      PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, LDA)
C
C     Cholesky factor A.
C
      CALL DPBTRF ('UPPER TRIANGLE OF A STORED', N, NDIAG, A, LDA,
     $             INFO)
      IF (INFO .NE. 0) THEN
        PRINT 1030, INFO
        STOP 1
      END IF
C
C     Compute and print the condition number of A.
C
      ANORM = 3.0D0
      CALL DPBCON ('UPPER TRIANGLE FACTOR STORED', N, NDIAG, A,
     $             LDA, ANORM, RCOND, WORK, IWORK, INFO)
      IF (INFO .NE. 0) THEN
        PRINT 1040, ABS(INFO)
        STOP 2
      END IF
      PRINT 1050, 1.0D0 / RCOND
C
 1000 FORMAT (1X, 'A in full form:')
 1010 FORMAT (4(3X, F6.3))
 1020 FORMAT (/1X, 'A in banded form:  (* in unused elements)')
 1030 FORMAT (1X, 'Error factoring A, INFO = ', I5)
 1040 FORMAT (1X, 'Illegal argument to DPBCON, argument #', I2)
 1050 FORMAT (/1X, 'Estimated condition number of A = ', F8.4)
C
      END
 

Sample Output

 
 A in full form:
    2.000   -1.000    0.000    0.000
   -1.000    2.000   -1.000    0.000
    0.000   -1.000    2.000   -1.000
    0.000    0.000   -1.000    2.000



 A in banded form:  (* in unused elements)
   ******   -1.000   -1.000   -1.000
    2.000    2.000    2.000    2.000



 Estimated condition number of A =   9.0000






Previous Next Contents Generated Index Doc Set Home