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