1      SUBROUTINE DPPCON( UPLO, N, AP, ANORM, RCOND, WORK, IWORK, INFO )
2*
3*  -- LAPACK routine (version 3.0) --
4*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5*     Courant Institute, Argonne National Lab, and Rice University
6*     March 31, 1993
7*
8*     .. Scalar Arguments ..
9      CHARACTER          UPLO
10      INTEGER            INFO, N
11      DOUBLE PRECISION   ANORM, RCOND
12*     ..
13*     .. Array Arguments ..
14      INTEGER            IWORK( * )
15      DOUBLE PRECISION   AP( * ), WORK( * )
16*     ..
17*
18*  Purpose
19*  =======
20*
21*  DPPCON estimates the reciprocal of the condition number (in the
22*  1-norm) of a real symmetric positive definite packed matrix using
23*  the Cholesky factorization A = U**T*U or A = L*L**T computed by
24*  DPPTRF.
25*
26*  An estimate is obtained for norm(inv(A)), and the reciprocal of the
27*  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
28*
29*  Arguments
30*  =========
31*
32*  UPLO    (input) CHARACTER*1
33*          = 'U':  Upper triangle of A is stored;
34*          = 'L':  Lower triangle of A is stored.
35*
36*  N       (input) INTEGER
37*          The order of the matrix A.  N >= 0.
38*
39*  AP      (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
40*          The triangular factor U or L from the Cholesky factorization
41*          A = U**T*U or A = L*L**T, packed columnwise in a linear
42*          array.  The j-th column of U or L is stored in the array AP
43*          as follows:
44*          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
45*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
46*
47*  ANORM   (input) DOUBLE PRECISION
48*          The 1-norm (or infinity-norm) of the symmetric matrix A.
49*
50*  RCOND   (output) DOUBLE PRECISION
51*          The reciprocal of the condition number of the matrix A,
52*          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
53*          estimate of the 1-norm of inv(A) computed in this routine.
54*
55*  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
56*
57*  IWORK   (workspace) INTEGER array, dimension (N)
58*
59*  INFO    (output) INTEGER
60*          = 0:  successful exit
61*          < 0:  if INFO = -i, the i-th argument had an illegal value
62*
63*  =====================================================================
64*
65*     .. Parameters ..
66      DOUBLE PRECISION   ONE, ZERO
67      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
68*     ..
69*     .. Local Scalars ..
70      LOGICAL            UPPER
71      CHARACTER          NORMIN
72      INTEGER            IX, KASE
73      DOUBLE PRECISION   AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
74*     ..
75*     .. External Functions ..
76      LOGICAL            LSAME
77      INTEGER            IDAMAX
78      DOUBLE PRECISION   DLAMCH
79      EXTERNAL           LSAME, IDAMAX, DLAMCH
80*     ..
81*     .. External Subroutines ..
82      EXTERNAL           DLACON, DLATPS, DRSCL, XERBLA
83*     ..
84*     .. Intrinsic Functions ..
85      INTRINSIC          ABS
86*     ..
87*     .. Executable Statements ..
88*
89*     Test the input parameters.
90*
91      INFO = 0
92      UPPER = LSAME( UPLO, 'U' )
93      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
94         INFO = -1
95      ELSE IF( N.LT.0 ) THEN
96         INFO = -2
97      ELSE IF( ANORM.LT.ZERO ) THEN
98         INFO = -4
99      END IF
100      IF( INFO.NE.0 ) THEN
101         CALL XERBLA( 'DPPCON', -INFO )
102         RETURN
103      END IF
104*
105*     Quick return if possible
106*
107      RCOND = ZERO
108      IF( N.EQ.0 ) THEN
109         RCOND = ONE
110         RETURN
111      ELSE IF( ANORM.EQ.ZERO ) THEN
112         RETURN
113      END IF
114*
115      SMLNUM = DLAMCH( 'Safe minimum' )
116*
117*     Estimate the 1-norm of the inverse.
118*
119      KASE = 0
120      NORMIN = 'N'
121   10 CONTINUE
122      CALL DLACON( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE )
123      IF( KASE.NE.0 ) THEN
124         IF( UPPER ) THEN
125*
126*           Multiply by inv(U').
127*
128            CALL DLATPS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
129     $                   AP, WORK, SCALEL, WORK( 2*N+1 ), INFO )
130            NORMIN = 'Y'
131*
132*           Multiply by inv(U).
133*
134            CALL DLATPS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
135     $                   AP, WORK, SCALEU, WORK( 2*N+1 ), INFO )
136         ELSE
137*
138*           Multiply by inv(L).
139*
140            CALL DLATPS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
141     $                   AP, WORK, SCALEL, WORK( 2*N+1 ), INFO )
142            NORMIN = 'Y'
143*
144*           Multiply by inv(L').
145*
146            CALL DLATPS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N,
147     $                   AP, WORK, SCALEU, WORK( 2*N+1 ), INFO )
148         END IF
149*
150*        Multiply by 1/SCALE if doing so will not cause overflow.
151*
152         SCALE = SCALEL*SCALEU
153         IF( SCALE.NE.ONE ) THEN
154            IX = IDAMAX( N, WORK, 1 )
155            IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
156     $         GO TO 20
157            CALL DRSCL( N, SCALE, WORK, 1 )
158         END IF
159         GO TO 10
160      END IF
161*
162*     Compute the estimate of the reciprocal condition number.
163*
164      IF( AINVNM.NE.ZERO )
165     $   RCOND = ( ONE / AINVNM ) / ANORM
166*
167   20 CONTINUE
168      RETURN
169*
170*     End of DPPCON
171*
172      END
173