001:       SUBROUTINE CPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
002:      $                   RWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH.
010: *
011: *     .. Scalar Arguments ..
012:       CHARACTER          UPLO
013:       INTEGER            INFO, KD, LDAB, N
014:       REAL               ANORM, RCOND
015: *     ..
016: *     .. Array Arguments ..
017:       REAL               RWORK( * )
018:       COMPLEX            AB( LDAB, * ), WORK( * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  CPBCON estimates the reciprocal of the condition number (in the
025: *  1-norm) of a complex Hermitian positive definite band matrix using
026: *  the Cholesky factorization A = U**H*U or A = L*L**H computed by
027: *  CPBTRF.
028: *
029: *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
030: *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
031: *
032: *  Arguments
033: *  =========
034: *
035: *  UPLO    (input) CHARACTER*1
036: *          = 'U':  Upper triangular factor stored in AB;
037: *          = 'L':  Lower triangular factor stored in AB.
038: *
039: *  N       (input) INTEGER
040: *          The order of the matrix A.  N >= 0.
041: *
042: *  KD      (input) INTEGER
043: *          The number of superdiagonals of the matrix A if UPLO = 'U',
044: *          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
045: *
046: *  AB      (input) COMPLEX array, dimension (LDAB,N)
047: *          The triangular factor U or L from the Cholesky factorization
048: *          A = U**H*U or A = L*L**H of the band matrix A, stored in the
049: *          first KD+1 rows of the array.  The j-th column of U or L is
050: *          stored in the j-th column of the array AB as follows:
051: *          if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
052: *          if UPLO ='L', AB(1+i-j,j)    = L(i,j) for j<=i<=min(n,j+kd).
053: *
054: *  LDAB    (input) INTEGER
055: *          The leading dimension of the array AB.  LDAB >= KD+1.
056: *
057: *  ANORM   (input) REAL
058: *          The 1-norm (or infinity-norm) of the Hermitian band matrix A.
059: *
060: *  RCOND   (output) REAL
061: *          The reciprocal of the condition number of the matrix A,
062: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
063: *          estimate of the 1-norm of inv(A) computed in this routine.
064: *
065: *  WORK    (workspace) COMPLEX array, dimension (2*N)
066: *
067: *  RWORK   (workspace) REAL array, dimension (N)
068: *
069: *  INFO    (output) INTEGER
070: *          = 0:  successful exit
071: *          < 0:  if INFO = -i, the i-th argument had an illegal value
072: *
073: *  =====================================================================
074: *
075: *     .. Parameters ..
076:       REAL               ONE, ZERO
077:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
078: *     ..
079: *     .. Local Scalars ..
080:       LOGICAL            UPPER
081:       CHARACTER          NORMIN
082:       INTEGER            IX, KASE
083:       REAL               AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
084:       COMPLEX            ZDUM
085: *     ..
086: *     .. Local Arrays ..
087:       INTEGER            ISAVE( 3 )
088: *     ..
089: *     .. External Functions ..
090:       LOGICAL            LSAME
091:       INTEGER            ICAMAX
092:       REAL               SLAMCH
093:       EXTERNAL           LSAME, ICAMAX, SLAMCH
094: *     ..
095: *     .. External Subroutines ..
096:       EXTERNAL           CLACN2, CLATBS, CSRSCL, XERBLA
097: *     ..
098: *     .. Intrinsic Functions ..
099:       INTRINSIC          ABS, AIMAG, REAL
100: *     ..
101: *     .. Statement Functions ..
102:       REAL               CABS1
103: *     ..
104: *     .. Statement Function definitions ..
105:       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
106: *     ..
107: *     .. Executable Statements ..
108: *
109: *     Test the input parameters.
110: *
111:       INFO = 0
112:       UPPER = LSAME( UPLO, 'U' )
113:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
114:          INFO = -1
115:       ELSE IF( N.LT.0 ) THEN
116:          INFO = -2
117:       ELSE IF( KD.LT.0 ) THEN
118:          INFO = -3
119:       ELSE IF( LDAB.LT.KD+1 ) THEN
120:          INFO = -5
121:       ELSE IF( ANORM.LT.ZERO ) THEN
122:          INFO = -6
123:       END IF
124:       IF( INFO.NE.0 ) THEN
125:          CALL XERBLA( 'CPBCON', -INFO )
126:          RETURN
127:       END IF
128: *
129: *     Quick return if possible
130: *
131:       RCOND = ZERO
132:       IF( N.EQ.0 ) THEN
133:          RCOND = ONE
134:          RETURN
135:       ELSE IF( ANORM.EQ.ZERO ) THEN
136:          RETURN
137:       END IF
138: *
139:       SMLNUM = SLAMCH( 'Safe minimum' )
140: *
141: *     Estimate the 1-norm of the inverse.
142: *
143:       KASE = 0
144:       NORMIN = 'N'
145:    10 CONTINUE
146:       CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
147:       IF( KASE.NE.0 ) THEN
148:          IF( UPPER ) THEN
149: *
150: *           Multiply by inv(U').
151: *
152:             CALL CLATBS( 'Upper', 'Conjugate transpose', 'Non-unit',
153:      $                   NORMIN, N, KD, AB, LDAB, WORK, SCALEL, RWORK,
154:      $                   INFO )
155:             NORMIN = 'Y'
156: *
157: *           Multiply by inv(U).
158: *
159:             CALL CLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
160:      $                   KD, AB, LDAB, WORK, SCALEU, RWORK, INFO )
161:          ELSE
162: *
163: *           Multiply by inv(L).
164: *
165:             CALL CLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
166:      $                   KD, AB, LDAB, WORK, SCALEL, RWORK, INFO )
167:             NORMIN = 'Y'
168: *
169: *           Multiply by inv(L').
170: *
171:             CALL CLATBS( 'Lower', 'Conjugate transpose', 'Non-unit',
172:      $                   NORMIN, N, KD, AB, LDAB, WORK, SCALEU, RWORK,
173:      $                   INFO )
174:          END IF
175: *
176: *        Multiply by 1/SCALE if doing so will not cause overflow.
177: *
178:          SCALE = SCALEL*SCALEU
179:          IF( SCALE.NE.ONE ) THEN
180:             IX = ICAMAX( N, WORK, 1 )
181:             IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
182:      $         GO TO 20
183:             CALL CSRSCL( N, SCALE, WORK, 1 )
184:          END IF
185:          GO TO 10
186:       END IF
187: *
188: *     Compute the estimate of the reciprocal condition number.
189: *
190:       IF( AINVNM.NE.ZERO )
191:      $   RCOND = ( ONE / AINVNM ) / ANORM
192: *
193:    20 CONTINUE
194: *
195:       RETURN
196: *
197: *     End of CPBCON
198: *
199:       END
200: