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