001:       SUBROUTINE DLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU,
002:      $                                NRHS, AB, LDAB, AFB, LDAFB, IPIV,
003:      $                                COLEQU, C, B, LDB, Y, LDY,
004:      $                                BERR_OUT, N_NORMS, ERRS_N, ERRS_C,
005:      $                                RES, AYB, DY, Y_TAIL, RCOND,
006:      $                                ITHRESH, RTHRESH, DZ_UB,
007:      $                                IGNORE_CWISE, INFO )
008: *
009: *     -- LAPACK routine (version 3.2)                                 --
010: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
011: *     -- Jason Riedy of Univ. of California Berkeley.                 --
012: *     -- November 2008                                                --
013: *
014: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
015: *     -- Univ. of California Berkeley and NAG Ltd.                    --
016: *
017:       IMPLICIT NONE
018: *     ..
019: *     .. Scalar Arguments ..
020:       INTEGER            INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS,
021:      $                   PREC_TYPE, TRANS_TYPE, N_NORMS, ITHRESH
022:       LOGICAL            COLEQU, IGNORE_CWISE
023:       DOUBLE PRECISION   RTHRESH, DZ_UB
024: *     ..
025: *     .. Array Arguments ..
026:       INTEGER            IPIV( * )
027:       DOUBLE PRECISION   AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
028:      $                   Y( LDY, * ), RES(*), DY(*), Y_TAIL(*)
029:       DOUBLE PRECISION   C( * ), AYB(*), RCOND, BERR_OUT(*),
030:      $                   ERRS_N( NRHS, * ), ERRS_C( NRHS, * )
031: *     ..
032: *     .. Local Scalars ..
033:       CHARACTER          TRANS
034:       INTEGER            CNT, I, J, M, X_STATE, Z_STATE, Y_PREC_STATE
035:       DOUBLE PRECISION   YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
036:      $                   DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX,
037:      $                   DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z,
038:      $                   EPS, HUGEVAL, INCR_THRESH
039:       LOGICAL            INCR_PREC
040: *     ..
041: *     .. Parameters ..
042:       INTEGER            UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
043:      $                   NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
044:      $                   EXTRA_Y
045:       PARAMETER          ( UNSTABLE_STATE = 0, WORKING_STATE = 1,
046:      $                   CONV_STATE = 2, NOPROG_STATE = 3 )
047:       PARAMETER          ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1,
048:      $                   EXTRA_Y = 2 )
049:       INTEGER            FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
050:       INTEGER            RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
051:       INTEGER            CMP_ERR_I, PIV_GROWTH_I
052:       PARAMETER          ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
053:      $                   BERR_I = 3 )
054:       PARAMETER          ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
055:       PARAMETER          ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
056:      $                   PIV_GROWTH_I = 9 )
057:       INTEGER            LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
058:      $                   LA_LINRX_CWISE_I
059:       PARAMETER          ( LA_LINRX_ITREF_I = 1,
060:      $                   LA_LINRX_ITHRESH_I = 2 )
061:       PARAMETER          ( LA_LINRX_CWISE_I = 3 )
062:       INTEGER            LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
063:      $                   LA_LINRX_RCOND_I
064:       PARAMETER          ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
065:       PARAMETER          ( LA_LINRX_RCOND_I = 3 )
066:       INTEGER            LA_LINRX_MAX_N_ERRS
067:       PARAMETER          ( LA_LINRX_MAX_N_ERRS = 3 )
068: *     ..
069: *     .. External Subroutines ..
070:       EXTERNAL           DAXPY, DCOPY, DGBTRS, DGBMV, BLAS_DGBMV_X,
071:      $                   BLAS_DGBMV2_X, DLA_GBAMV, DLA_WWADDW, DLAMCH,
072:      $                   CHLA_TRANSTYPE, DLA_LIN_BERR
073:       DOUBLE PRECISION   DLAMCH
074:       CHARACTER          CHLA_TRANSTYPE
075: *     ..
076: *     .. Intrinsic Functions ..
077:       INTRINSIC          ABS, MAX, MIN
078: *     ..
079: *     .. Executable Statements ..
080: *
081:       IF (INFO.NE.0) RETURN
082:       TRANS = CHLA_TRANSTYPE(TRANS_TYPE)
083:       EPS = DLAMCH( 'Epsilon' )
084:       HUGEVAL = DLAMCH( 'Overflow' )
085: *     Force HUGEVAL to Inf
086:       HUGEVAL = HUGEVAL * HUGEVAL
087: *     Using HUGEVAL may lead to spurious underflows.
088:       INCR_THRESH = DBLE( N ) * EPS
089:       M = KL+KU+1
090: 
091:       DO J = 1, NRHS
092:          Y_PREC_STATE = EXTRA_RESIDUAL
093:          IF ( Y_PREC_STATE .EQ. EXTRA_Y ) THEN
094:             DO I = 1, N
095:                Y_TAIL( I ) = 0.0D+0
096:             END DO
097:          END IF
098: 
099:          DXRAT = 0.0D+0
100:          DXRATMAX = 0.0D+0
101:          DZRAT = 0.0D+0
102:          DZRATMAX = 0.0D+0
103:          FINAL_DX_X = HUGEVAL
104:          FINAL_DZ_Z = HUGEVAL
105:          PREVNORMDX = HUGEVAL
106:          PREV_DZ_Z = HUGEVAL
107:          DZ_Z = HUGEVAL
108:          DX_X = HUGEVAL
109: 
110:          X_STATE = WORKING_STATE
111:          Z_STATE = UNSTABLE_STATE
112:          INCR_PREC = .FALSE.
113: 
114:          DO CNT = 1, ITHRESH
115: *
116: *        Compute residual RES = B_s - op(A_s) * Y,
117: *            op(A) = A, A**T, or A**H depending on TRANS (and type).
118: *
119:             CALL DCOPY( N, B( 1, J ), 1, RES, 1 )
120:             IF ( Y_PREC_STATE .EQ. BASE_RESIDUAL ) THEN
121:                CALL DGBMV( TRANS, M, N, KL, KU, -1.0D+0, AB, LDAB,
122:      $              Y( 1, J ), 1, 1.0D+0, RES, 1 )
123:             ELSE IF ( Y_PREC_STATE .EQ. EXTRA_RESIDUAL ) THEN
124:                CALL BLAS_DGBMV_X( TRANS_TYPE, N, N, KL, KU,
125:      $              -1.0D+0, AB, LDAB, Y( 1, J ), 1, 1.0D+0, RES, 1,
126:      $              PREC_TYPE )
127:             ELSE
128:                CALL BLAS_DGBMV2_X( TRANS_TYPE, N, N, KL, KU, -1.0D+0,
129:      $              AB, LDAB, Y( 1, J ), Y_TAIL, 1, 1.0D+0, RES, 1,
130:      $              PREC_TYPE )
131:             END IF
132: 
133: !        XXX: RES is no longer needed.
134:             CALL DCOPY( N, RES, 1, DY, 1 )
135:             CALL DGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV, DY, N,
136:      $           INFO )
137: *
138: *         Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
139: *
140:             NORMX = 0.0D+0
141:             NORMY = 0.0D+0
142:             NORMDX = 0.0D+0
143:             DZ_Z = 0.0D+0
144:             YMIN = HUGEVAL
145: 
146:             DO I = 1, N
147:                YK = ABS( Y( I, J ) )
148:                DYK = ABS( DY( I ) )
149: 
150:                IF ( YK .NE. 0.0D+0 ) THEN
151:                   DZ_Z = MAX( DZ_Z, DYK / YK )
152:                ELSE IF ( DYK .NE. 0.0D+0 ) THEN
153:                   DZ_Z = HUGEVAL
154:                END IF
155: 
156:                YMIN = MIN( YMIN, YK )
157: 
158:                NORMY = MAX( NORMY, YK )
159: 
160:                IF ( COLEQU ) THEN
161:                   NORMX = MAX( NORMX, YK * C( I ) )
162:                   NORMDX = MAX( NORMDX, DYK * C( I ) )
163:                ELSE
164:                   NORMX = NORMY
165:                   NORMDX = MAX( NORMDX, DYK )
166:                END IF
167:             END DO
168: 
169:             IF ( NORMX .NE. 0.0D+0 ) THEN
170:                DX_X = NORMDX / NORMX
171:             ELSE IF ( NORMDX .EQ. 0.0D+0 ) THEN
172:                DX_X = 0.0D+0
173:             ELSE
174:                DX_X = HUGEVAL
175:             END IF
176: 
177:             DXRAT = NORMDX / PREVNORMDX
178:             DZRAT = DZ_Z / PREV_DZ_Z
179: *
180: *         Check termination criteria.
181: *
182:             IF ( .NOT.IGNORE_CWISE
183:      $           .AND. YMIN*RCOND .LT. INCR_THRESH*NORMY
184:      $           .AND. Y_PREC_STATE .LT. EXTRA_Y )
185:      $           INCR_PREC = .TRUE.
186: 
187:             IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH )
188:      $           X_STATE = WORKING_STATE
189:             IF ( X_STATE .EQ. WORKING_STATE ) THEN
190:                IF ( DX_X .LE. EPS ) THEN
191:                   X_STATE = CONV_STATE
192:                ELSE IF ( DXRAT .GT. RTHRESH ) THEN
193:                   IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
194:                      INCR_PREC = .TRUE.
195:                   ELSE
196:                      X_STATE = NOPROG_STATE
197:                   END IF
198:                ELSE
199:                   IF ( DXRAT .GT. DXRATMAX ) DXRATMAX = DXRAT
200:                END IF
201:                IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X
202:             END IF
203: 
204:             IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB )
205:      $           Z_STATE = WORKING_STATE
206:             IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH )
207:      $           Z_STATE = WORKING_STATE
208:             IF ( Z_STATE .EQ. WORKING_STATE ) THEN
209:                IF ( DZ_Z .LE. EPS ) THEN
210:                   Z_STATE = CONV_STATE
211:                ELSE IF ( DZ_Z .GT. DZ_UB ) THEN
212:                   Z_STATE = UNSTABLE_STATE
213:                   DZRATMAX = 0.0D+0
214:                   FINAL_DZ_Z = HUGEVAL
215:                ELSE IF ( DZRAT .GT. RTHRESH ) THEN
216:                   IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
217:                      INCR_PREC = .TRUE.
218:                   ELSE
219:                      Z_STATE = NOPROG_STATE
220:                   END IF
221:                ELSE
222:                   IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT
223:                END IF
224:                IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
225:             END IF
226: *
227: *           Exit if both normwise and componentwise stopped working,
228: *           but if componentwise is unstable, let it go at least two
229: *           iterations.
230: *
231:             IF ( X_STATE.NE.WORKING_STATE ) THEN
232:                IF ( IGNORE_CWISE ) GOTO 666
233:                IF ( Z_STATE.EQ.NOPROG_STATE .OR. Z_STATE.EQ.CONV_STATE )
234:      $              GOTO 666
235:                IF ( Z_STATE.EQ.UNSTABLE_STATE .AND. CNT.GT.1 ) GOTO 666
236:             END IF
237: 
238:             IF ( INCR_PREC ) THEN
239:                INCR_PREC = .FALSE.
240:                Y_PREC_STATE = Y_PREC_STATE + 1
241:                DO I = 1, N
242:                   Y_TAIL( I ) = 0.0D+0
243:                END DO
244:             END IF
245: 
246:             PREVNORMDX = NORMDX
247:             PREV_DZ_Z = DZ_Z
248: *
249: *           Update soluton.
250: *
251:             IF (Y_PREC_STATE .LT. EXTRA_Y) THEN
252:                CALL DAXPY( N, 1.0D+0, DY, 1, Y(1,J), 1 )
253:             ELSE
254:                CALL DLA_WWADDW( N, Y(1,J), Y_TAIL, DY )
255:             END IF
256: 
257:          END DO
258: *        Target of "IF (Z_STOP .AND. X_STOP)".  Sun's f77 won't EXIT.
259:  666     CONTINUE
260: *
261: *     Set final_* when cnt hits ithresh.
262: *
263:          IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X
264:          IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
265: *
266: *     Compute error bounds.
267: *
268:          IF ( N_NORMS .GE. 1 ) THEN
269:             ERRS_N( J, LA_LINRX_ERR_I ) = FINAL_DX_X / (1 - DXRATMAX)
270:          END IF
271:          IF (N_NORMS .GE. 2) THEN
272:             ERRS_C( J, LA_LINRX_ERR_I ) = FINAL_DZ_Z / (1 - DZRATMAX)
273:          END IF
274: *
275: *     Compute componentwise relative backward error from formula
276: *         max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
277: *     where abs(Z) is the componentwise absolute value of the matrix
278: *     or vector Z.
279: *
280: *        Compute residual RES = B_s - op(A_s) * Y,
281: *            op(A) = A, A**T, or A**H depending on TRANS (and type).
282: *
283:          CALL DCOPY( N, B( 1, J ), 1, RES, 1 )
284:          CALL DGBMV(TRANS, N, N, KL, KU, -1.0D+0, AB, LDAB, Y(1,J),
285:      $        1, 1.0D+0, RES, 1 )
286: 
287:          DO I = 1, N
288:             AYB( I ) = ABS( B( I, J ) )
289:          END DO
290: *
291: *     Compute abs(op(A_s))*abs(Y) + abs(B_s).
292: *
293:         CALL DLA_GBAMV( TRANS_TYPE, N, N, KL, KU, 1.0D+0,
294:      $        AB, LDAB, Y(1, J), 1, 1.0D+0, AYB, 1 )
295: 
296:          CALL DLA_LIN_BERR( N, N, 1, RES, AYB, BERR_OUT( J ) )
297: *
298: *     End of loop for each RHS
299: *
300:       END DO
301: *
302:       RETURN
303:       END
304: