cspco(3P)

NAME

cspco - compute the UDU factorization and condition number of a symmetric matrix A in packed storage. If the condition number is not needed then xSPFA is slightly faster. It is typical to follow a call to xSPCO with a call to xSPSL to solve Ax = b or to xSPDI to compute the determinant, inverse, and inertia of A.

SYNOPSIS

CALL DSPCO (DA, N, IPIVOT, DRCOND, DWORK)

CALL SSPCO (SA, N, IPIVOT, SRCOND, SWORK)

CALL ZSPCO (ZA, N, IPIVOT, DRCOND, ZWORK)

CALL CSPCO (CA, N, IPIVOT, SRCOND, CWORK)

void dspco(double ∗dap, long int n, long int ∗kpvt, double ∗rcond)

void sspco(float ∗sap, long int n, long int ∗kpvt, float ∗rcond)

void zspco(doublecomplex ∗zap, long int n,
long int ∗kpvt, double ∗rcond)

void cspco(complex ∗cap, long int n, long int ∗kpvt,
float ∗rcond)

ARGUMENTS

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

NOrder of the matrix A. N >= 0.

IPIVOTOn exit, a vector of pivot indices.

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.

SAMPLE PROGRAM

      PROGRAM TEST
      IMPLICIT NONE
C
      INTEGER           LENGTA, N
      PARAMETER        (N = 3)
      PARAMETER        (LENGTA = (N ∗ N + N) / 2)
C
      DOUBLE PRECISION A(LENGTA), B(N), RCOND, WORK(N)
      INTEGER           IPIVOT(N)
C
      EXTERNAL          DSPCO, DSPSL
C
C     Initialize the array A to store in packed symmetric format
C     the matrix A shown below. Initialize the array B to store
C     the vector b shown below.
C
C         1 0 4        30
C     A = 0 2 0    b = 4
C         4 0 1        15
C
      DATA A / 1.0D0, 0.0D0, 2.0D0, 4.0D0, 0.0D0, 1.0D0 /
      DATA B / 3.0D1, 4.0D0, 1.5D1 /
C
      PRINT 1000
      PRINT 1010, A(1), A(2), A(4)
      PRINT 1010, A(2), A(3), A(5)
      PRINT 1010, A(4), A(5), A(6)
      PRINT 1020
      PRINT 1030, B
      CALL DSPCO (A, N, IPIVOT, RCOND, WORK)
      IF ((RCOND + 1.0D0) .EQ. RCOND) THEN
        PRINT 1040
      END IF
      CALL DSPSL (A, N, IPIVOT, B)
      PRINT 1050
      PRINT 1030, B
C
1000 FORMAT (1X, ’A:’)
1010 FORMAT (3(3X, F4.1))
1020 FORMAT (/1X, ’b:’)
1030 FORMAT (3X, F4.1)
1040 FORMAT (1X, ’A may be singular to working precision.’)
1050 FORMAT (/1X, ’A∗∗(-1) ∗ b:’)
C
      END

SAMPLE OUTPUT

A:
    1.0    0.0    4.0
    0.0    2.0    0.0
    4.0    0.0    1.0

b:
   30.0
    4.0
   15.0

A∗∗(-1) ∗ b:
    2.0
    2.0
    7.0