001:       SUBROUTINE SGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
002:      $                   WORK, IWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     Modified to call SLACN2 in place of SLACON, 7 Feb 03, SJH.
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          NORM
012:       INTEGER            INFO, KL, KU, LDAB, N
013:       REAL               ANORM, RCOND
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            IPIV( * ), IWORK( * )
017:       REAL               AB( LDAB, * ), WORK( * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  SGBCON estimates the reciprocal of the condition number of a real
024: *  general band matrix A, in either the 1-norm or the infinity-norm,
025: *  using the LU factorization computed by SGBTRF.
026: *
027: *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
028: *  condition number is computed as
029: *     RCOND = 1 / ( norm(A) * norm(inv(A)) ).
030: *
031: *  Arguments
032: *  =========
033: *
034: *  NORM    (input) CHARACTER*1
035: *          Specifies whether the 1-norm condition number or the
036: *          infinity-norm condition number is required:
037: *          = '1' or 'O':  1-norm;
038: *          = 'I':         Infinity-norm.
039: *
040: *  N       (input) INTEGER
041: *          The order of the matrix A.  N >= 0.
042: *
043: *  KL      (input) INTEGER
044: *          The number of subdiagonals within the band of A.  KL >= 0.
045: *
046: *  KU      (input) INTEGER
047: *          The number of superdiagonals within the band of A.  KU >= 0.
048: *
049: *  AB      (input) REAL array, dimension (LDAB,N)
050: *          Details of the LU factorization of the band matrix A, as
051: *          computed by SGBTRF.  U is stored as an upper triangular band
052: *          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
053: *          the multipliers used during the factorization are stored in
054: *          rows KL+KU+2 to 2*KL+KU+1.
055: *
056: *  LDAB    (input) INTEGER
057: *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
058: *
059: *  IPIV    (input) INTEGER array, dimension (N)
060: *          The pivot indices; for 1 <= i <= N, row i of the matrix was
061: *          interchanged with row IPIV(i).
062: *
063: *  ANORM   (input) REAL
064: *          If NORM = '1' or 'O', the 1-norm of the original matrix A.
065: *          If NORM = 'I', the infinity-norm of the original matrix A.
066: *
067: *  RCOND   (output) REAL
068: *          The reciprocal of the condition number of the matrix A,
069: *          computed as RCOND = 1/(norm(A) * norm(inv(A))).
070: *
071: *  WORK    (workspace) REAL array, dimension (3*N)
072: *
073: *  IWORK   (workspace) INTEGER array, dimension (N)
074: *
075: *  INFO    (output) INTEGER
076: *          = 0:  successful exit
077: *          < 0: if INFO = -i, the i-th argument had an illegal value
078: *
079: *  =====================================================================
080: *
081: *     .. Parameters ..
082:       REAL               ONE, ZERO
083:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
084: *     ..
085: *     .. Local Scalars ..
086:       LOGICAL            LNOTI, ONENRM
087:       CHARACTER          NORMIN
088:       INTEGER            IX, J, JP, KASE, KASE1, KD, LM
089:       REAL               AINVNM, SCALE, SMLNUM, T
090: *     ..
091: *     .. Local Arrays ..
092:       INTEGER            ISAVE( 3 )
093: *     ..
094: *     .. External Functions ..
095:       LOGICAL            LSAME
096:       INTEGER            ISAMAX
097:       REAL               SDOT, SLAMCH
098:       EXTERNAL           LSAME, ISAMAX, SDOT, SLAMCH
099: *     ..
100: *     .. External Subroutines ..
101:       EXTERNAL           SAXPY, SLACN2, SLATBS, SRSCL, XERBLA
102: *     ..
103: *     .. Intrinsic Functions ..
104:       INTRINSIC          ABS, MIN
105: *     ..
106: *     .. Executable Statements ..
107: *
108: *     Test the input parameters.
109: *
110:       INFO = 0
111:       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
112:       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
113:          INFO = -1
114:       ELSE IF( N.LT.0 ) THEN
115:          INFO = -2
116:       ELSE IF( KL.LT.0 ) THEN
117:          INFO = -3
118:       ELSE IF( KU.LT.0 ) THEN
119:          INFO = -4
120:       ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
121:          INFO = -6
122:       ELSE IF( ANORM.LT.ZERO ) THEN
123:          INFO = -8
124:       END IF
125:       IF( INFO.NE.0 ) THEN
126:          CALL XERBLA( 'SGBCON', -INFO )
127:          RETURN
128:       END IF
129: *
130: *     Quick return if possible
131: *
132:       RCOND = ZERO
133:       IF( N.EQ.0 ) THEN
134:          RCOND = ONE
135:          RETURN
136:       ELSE IF( ANORM.EQ.ZERO ) THEN
137:          RETURN
138:       END IF
139: *
140:       SMLNUM = SLAMCH( 'Safe minimum' )
141: *
142: *     Estimate the norm of inv(A).
143: *
144:       AINVNM = ZERO
145:       NORMIN = 'N'
146:       IF( ONENRM ) THEN
147:          KASE1 = 1
148:       ELSE
149:          KASE1 = 2
150:       END IF
151:       KD = KL + KU + 1
152:       LNOTI = KL.GT.0
153:       KASE = 0
154:    10 CONTINUE
155:       CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
156:       IF( KASE.NE.0 ) THEN
157:          IF( KASE.EQ.KASE1 ) THEN
158: *
159: *           Multiply by inv(L).
160: *
161:             IF( LNOTI ) THEN
162:                DO 20 J = 1, N - 1
163:                   LM = MIN( KL, N-J )
164:                   JP = IPIV( J )
165:                   T = WORK( JP )
166:                   IF( JP.NE.J ) THEN
167:                      WORK( JP ) = WORK( J )
168:                      WORK( J ) = T
169:                   END IF
170:                   CALL SAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
171:    20          CONTINUE
172:             END IF
173: *
174: *           Multiply by inv(U).
175: *
176:             CALL SLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
177:      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
178:      $                   INFO )
179:          ELSE
180: *
181: *           Multiply by inv(U').
182: *
183:             CALL SLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
184:      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
185:      $                   INFO )
186: *
187: *           Multiply by inv(L').
188: *
189:             IF( LNOTI ) THEN
190:                DO 30 J = N - 1, 1, -1
191:                   LM = MIN( KL, N-J )
192:                   WORK( J ) = WORK( J ) - SDOT( LM, AB( KD+1, J ), 1,
193:      $                        WORK( J+1 ), 1 )
194:                   JP = IPIV( J )
195:                   IF( JP.NE.J ) THEN
196:                      T = WORK( JP )
197:                      WORK( JP ) = WORK( J )
198:                      WORK( J ) = T
199:                   END IF
200:    30          CONTINUE
201:             END IF
202:          END IF
203: *
204: *        Divide X by 1/SCALE if doing so will not cause overflow.
205: *
206:          NORMIN = 'Y'
207:          IF( SCALE.NE.ONE ) THEN
208:             IX = ISAMAX( N, WORK, 1 )
209:             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
210:      $         GO TO 40
211:             CALL SRSCL( N, SCALE, WORK, 1 )
212:          END IF
213:          GO TO 10
214:       END IF
215: *
216: *     Compute the estimate of the reciprocal condition number.
217: *
218:       IF( AINVNM.NE.ZERO )
219:      $   RCOND = ( ONE / AINVNM ) / ANORM
220: *
221:    40 CONTINUE
222:       RETURN
223: *
224: *     End of SGBCON
225: *
226:       END
227: