001:SUBROUTINEZGBSVX( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, 002: $ LDAFB, IPIV, EQUED, R, C, B, LDB, X, LDX, 003: $ RCOND, FERR, BERR, WORK, RWORK, INFO ) 004:*005:* -- LAPACK driver routine (version 3.2) --006:* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..007:* November 2006008:*009:* .. Scalar Arguments ..010: CHARACTER EQUED, FACT, TRANS 011: INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS 012: DOUBLE PRECISION RCOND 013:* ..014:* .. Array Arguments ..015: INTEGERIPIV( * ) 016: DOUBLE PRECISIONBERR( * ),C( * ),FERR( * ),R( * ), 017: $RWORK( * ) 018: COMPLEX*16AB( LDAB, * ),AFB( LDAFB, * ),B( LDB, * ), 019: $WORK( * ),X( LDX, * ) 020:* ..021:*022:* Purpose023:* =======024:*025:* ZGBSVX uses the LU factorization to compute the solution to a complex026:* system of linear equations A * X = B, A**T * X = B, or A**H * X = B,027:* where A is a band matrix of order N with KL subdiagonals and KU028:* superdiagonals, and X and B are N-by-NRHS matrices.029:*030:* Error bounds on the solution and a condition estimate are also031:* provided.032:*033:* Description034:* ===========035:*036:* The following steps are performed by this subroutine:037:*038:* 1. If FACT = 'E', real scaling factors are computed to equilibrate039:* the system:040:* TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B041:* TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B042:* TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B043:* Whether or not the system will be equilibrated depends on the044:* scaling of the matrix A, but if equilibration is used, A is045:* overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')046:* or diag(C)*B (if TRANS = 'T' or 'C').047:*048:* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the049:* matrix A (after equilibration if FACT = 'E') as050:* A = L * U,051:* where L is a product of permutation and unit lower triangular052:* matrices with KL subdiagonals, and U is upper triangular with053:* KL+KU superdiagonals.054:*055:* 3. If some U(i,i)=0, so that U is exactly singular, then the routine056:* returns with INFO = i. Otherwise, the factored form of A is used057:* to estimate the condition number of the matrix A. If the058:* reciprocal of the condition number is less than machine precision,059:* INFO = N+1 is returned as a warning, but the routine still goes on060:* to solve for X and compute error bounds as described below.061:*062:* 4. The system of equations is solved for X using the factored form063:* of A.064:*065:* 5. Iterative refinement is applied to improve the computed solution066:* matrix and calculate error bounds and backward error estimates067:* for it.068:*069:* 6. If equilibration was used, the matrix X is premultiplied by070:* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so071:* that it solves the original system before equilibration.072:*073:* Arguments074:* =========075:*076:* FACT (input) CHARACTER*1077:* Specifies whether or not the factored form of the matrix A is078:* supplied on entry, and if not, whether the matrix A should be079:* equilibrated before it is factored.080:* = 'F': On entry, AFB and IPIV contain the factored form of081:* A. If EQUED is not 'N', the matrix A has been082:* equilibrated with scaling factors given by R and C.083:* AB, AFB, and IPIV are not modified.084:* = 'N': The matrix A will be copied to AFB and factored.085:* = 'E': The matrix A will be equilibrated if necessary, then086:* copied to AFB and factored.087:*088:* TRANS (input) CHARACTER*1089:* Specifies the form of the system of equations.090:* = 'N': A * X = B (No transpose)091:* = 'T': A**T * X = B (Transpose)092:* = 'C': A**H * X = B (Conjugate transpose)093:*094:* N (input) INTEGER095:* The number of linear equations, i.e., the order of the096:* matrix A. N >= 0.097:*098:* KL (input) INTEGER099:* The number of subdiagonals within the band of A. KL >= 0.100:*101:* KU (input) INTEGER102:* The number of superdiagonals within the band of A. KU >= 0.103:*104:* NRHS (input) INTEGER105:* The number of right hand sides, i.e., the number of columns106:* of the matrices B and X. NRHS >= 0.107:*108:* AB (input/output) COMPLEX*16 array, dimension (LDAB,N)109:* On entry, the matrix A in band storage, in rows 1 to KL+KU+1.110:* The j-th column of A is stored in the j-th column of the111:* array AB as follows:112:* AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)113:*114:* If FACT = 'F' and EQUED is not 'N', then A must have been115:* equilibrated by the scaling factors in R and/or C. AB is not116:* modified if FACT = 'F' or 'N', or if FACT = 'E' and117:* EQUED = 'N' on exit.118:*119:* On exit, if EQUED .ne. 'N', A is scaled as follows:120:* EQUED = 'R': A := diag(R) * A121:* EQUED = 'C': A := A * diag(C)122:* EQUED = 'B': A := diag(R) * A * diag(C).123:*124:* LDAB (input) INTEGER125:* The leading dimension of the array AB. LDAB >= KL+KU+1.126:*127:* AFB (input or output) COMPLEX*16 array, dimension (LDAFB,N)128:* If FACT = 'F', then AFB is an input argument and on entry129:* contains details of the LU factorization of the band matrix130:* A, as computed by ZGBTRF. U is stored as an upper triangular131:* band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,132:* and the multipliers used during the factorization are stored133:* in rows KL+KU+2 to 2*KL+KU+1. If EQUED .ne. 'N', then AFB is134:* the factored form of the equilibrated matrix A.135:*136:* If FACT = 'N', then AFB is an output argument and on exit137:* returns details of the LU factorization of A.138:*139:* If FACT = 'E', then AFB is an output argument and on exit140:* returns details of the LU factorization of the equilibrated141:* matrix A (see the description of AB for the form of the142:* equilibrated matrix).143:*144:* LDAFB (input) INTEGER145:* The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.146:*147:* IPIV (input or output) INTEGER array, dimension (N)148:* If FACT = 'F', then IPIV is an input argument and on entry149:* contains the pivot indices from the factorization A = L*U150:* as computed by ZGBTRF; row i of the matrix was interchanged151:* with row IPIV(i).152:*153:* If FACT = 'N', then IPIV is an output argument and on exit154:* contains the pivot indices from the factorization A = L*U155:* of the original matrix A.156:*157:* If FACT = 'E', then IPIV is an output argument and on exit158:* contains the pivot indices from the factorization A = L*U159:* of the equilibrated matrix A.160:*161:* EQUED (input or output) CHARACTER*1162:* Specifies the form of equilibration that was done.163:* = 'N': No equilibration (always true if FACT = 'N').164:* = 'R': Row equilibration, i.e., A has been premultiplied by165:* diag(R).166:* = 'C': Column equilibration, i.e., A has been postmultiplied167:* by diag(C).168:* = 'B': Both row and column equilibration, i.e., A has been169:* replaced by diag(R) * A * diag(C).170:* EQUED is an input argument if FACT = 'F'; otherwise, it is an171:* output argument.172:*173:* R (input or output) DOUBLE PRECISION array, dimension (N)174:* The row scale factors for A. If EQUED = 'R' or 'B', A is175:* multiplied on the left by diag(R); if EQUED = 'N' or 'C', R176:* is not accessed. R is an input argument if FACT = 'F';177:* otherwise, R is an output argument. If FACT = 'F' and178:* EQUED = 'R' or 'B', each element of R must be positive.179:*180:* C (input or output) DOUBLE PRECISION array, dimension (N)181:* The column scale factors for A. If EQUED = 'C' or 'B', A is182:* multiplied on the right by diag(C); if EQUED = 'N' or 'R', C183:* is not accessed. C is an input argument if FACT = 'F';184:* otherwise, C is an output argument. If FACT = 'F' and185:* EQUED = 'C' or 'B', each element of C must be positive.186:*187:* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)188:* On entry, the right hand side matrix B.189:* On exit,190:* if EQUED = 'N', B is not modified;191:* if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by192:* diag(R)*B;193:* if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is194:* overwritten by diag(C)*B.195:*196:* LDB (input) INTEGER197:* The leading dimension of the array B. LDB >= max(1,N).198:*199:* X (output) COMPLEX*16 array, dimension (LDX,NRHS)200:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X201:* to the original system of equations. Note that A and B are202:* modified on exit if EQUED .ne. 'N', and the solution to the203:* equilibrated system is inv(diag(C))*X if TRANS = 'N' and204:* EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'205:* and EQUED = 'R' or 'B'.206:*207:* LDX (input) INTEGER208:* The leading dimension of the array X. LDX >= max(1,N).209:*210:* RCOND (output) DOUBLE PRECISION211:* The estimate of the reciprocal condition number of the matrix212:* A after equilibration (if done). If RCOND is less than the213:* machine precision (in particular, if RCOND = 0), the matrix214:* is singular to working precision. This condition is215:* indicated by a return code of INFO > 0.216:*217:* FERR (output) DOUBLE PRECISION array, dimension (NRHS)218:* The estimated forward error bound for each solution vector219:* X(j) (the j-th column of the solution matrix X).220:* If XTRUE is the true solution corresponding to X(j), FERR(j)221:* is an estimated upper bound for the magnitude of the largest222:* element in (X(j) - XTRUE) divided by the magnitude of the223:* largest element in X(j). The estimate is as reliable as224:* the estimate for RCOND, and is almost always a slight225:* overestimate of the true error.226:*227:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)228:* The componentwise relative backward error of each solution229:* vector X(j) (i.e., the smallest relative change in230:* any element of A or B that makes X(j) an exact solution).231:*232:* WORK (workspace) COMPLEX*16 array, dimension (2*N)233:*234:* RWORK (workspace/output) DOUBLE PRECISION array, dimension (N)235:* On exit, RWORK(1) contains the reciprocal pivot growth236:* factor norm(A)/norm(U). The "max absolute element" norm is237:* used. If RWORK(1) is much less than 1, then the stability238:* of the LU factorization of the (equilibrated) matrix A239:* could be poor. This also means that the solution X, condition240:* estimator RCOND, and forward error bound FERR could be241:* unreliable. If factorization fails with 0<INFO<=N, then242:* RWORK(1) contains the reciprocal pivot growth factor for the243:* leading INFO columns of A.244:*245:* INFO (output) INTEGER246:* = 0: successful exit247:* < 0: if INFO = -i, the i-th argument had an illegal value248:* > 0: if INFO = i, and i is249:* <= N: U(i,i) is exactly zero. The factorization250:* has been completed, but the factor U is exactly251:* singular, so the solution and error bounds252:* could not be computed. RCOND = 0 is returned.253:* = N+1: U is nonsingular, but RCOND is less than machine254:* precision, meaning that the matrix is singular255:* to working precision. Nevertheless, the256:* solution and error bounds are computed because257:* there are a number of situations where the258:* computed solution can be more accurate than the259:* value of RCOND would suggest.260:*261:* =====================================================================262:* Moved setting of INFO = N+1 so INFO does not subsequently get263:* overwritten. Sven, 17 Mar 05.264:* =====================================================================265:*266:* .. Parameters ..267: DOUBLE PRECISION ZERO, ONE 268:PARAMETER( ZERO = 0.0D+0, ONE = 1.0D+0 ) 269:* ..270:* .. Local Scalars ..271:LOGICALCOLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU 272: CHARACTER NORM 273: INTEGER I, INFEQU, J, J1, J2 274: DOUBLE PRECISION AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN, 275: $ ROWCND, RPVGRW, SMLNUM 276:* ..277:* .. External Functions ..278:LOGICALLSAME 279: DOUBLE PRECISION DLAMCH, ZLANGB, ZLANTB 280:EXTERNALLSAME, DLAMCH, ZLANGB, ZLANTB 281:* ..282:* .. External Subroutines ..283:EXTERNALXERBLA, ZCOPY, ZGBCON, ZGBEQU, ZGBRFS, ZGBTRF, 284: $ ZGBTRS, ZLACPY, ZLAQGB 285:* ..286:* .. Intrinsic Functions ..287:INTRINSICABS, MAX, MIN 288:* ..289:* .. Executable Statements ..290:*291: INFO = 0 292: NOFACT =LSAME( FACT, 'N' ) 293: EQUIL =LSAME( FACT, 'E' ) 294: NOTRAN =LSAME( TRANS, 'N' ) 295:IF( NOFACT .OR. EQUIL )THEN296: EQUED = 'N' 297: ROWEQU = .FALSE. 298: COLEQU = .FALSE. 299:ELSE300: ROWEQU =LSAME( EQUED, 'R' ) .OR.LSAME( EQUED, 'B' ) 301: COLEQU =LSAME( EQUED, 'C' ) .OR.LSAME( EQUED, 'B' ) 302: SMLNUM =DLAMCH( 'Safe minimum' ) 303: BIGNUM = ONE / SMLNUM 304:ENDIF305:*306:* Test the input parameters.307:*308:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 309: $THEN310: INFO = -1 311:ELSEIF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 312: $LSAME( TRANS, 'C' ) )THEN313: INFO = -2 314:ELSEIF( N.LT.0 )THEN315: INFO = -3 316:ELSEIF( KL.LT.0 )THEN317: INFO = -4 318:ELSEIF( KU.LT.0 )THEN319: INFO = -5 320:ELSEIF( NRHS.LT.0 )THEN321: INFO = -6 322:ELSEIF( LDAB.LT.KL+KU+1 )THEN323: INFO = -8 324:ELSEIF( LDAFB.LT.2*KL+KU+1 )THEN325: INFO = -10 326:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 327: $ ( ROWEQU .OR. COLEQU .OR.LSAME( EQUED, 'N' ) ) )THEN328: INFO = -12 329:ELSE330:IF( ROWEQU )THEN331: RCMIN = BIGNUM 332: RCMAX = ZERO 333:DO10 J = 1, N 334: RCMIN =MIN( RCMIN,R( J ) ) 335: RCMAX =MAX( RCMAX,R( J ) ) 336: 10CONTINUE337:IF( RCMIN.LE.ZERO )THEN338: INFO = -13 339:ELSEIF( N.GT.0 )THEN340: ROWCND =MAX( RCMIN, SMLNUM ) /MIN( RCMAX, BIGNUM ) 341:ELSE342: ROWCND = ONE 343:ENDIF344:ENDIF345:IF( COLEQU .AND. INFO.EQ.0 )THEN346: RCMIN = BIGNUM 347: RCMAX = ZERO 348:DO20 J = 1, N 349: RCMIN =MIN( RCMIN,C( J ) ) 350: RCMAX =MAX( RCMAX,C( J ) ) 351: 20CONTINUE352:IF( RCMIN.LE.ZERO )THEN353: INFO = -14 354:ELSEIF( N.GT.0 )THEN355: COLCND =MAX( RCMIN, SMLNUM ) /MIN( RCMAX, BIGNUM ) 356:ELSE357: COLCND = ONE 358:ENDIF359:ENDIF360:IF( INFO.EQ.0 )THEN361:IF( LDB.LT.MAX( 1, N ) )THEN362: INFO = -16 363:ELSEIF( LDX.LT.MAX( 1, N ) )THEN364: INFO = -18 365:ENDIF366:ENDIF367:ENDIF368:*369:IF( INFO.NE.0 )THEN370:CALLXERBLA( 'ZGBSVX', -INFO ) 371:RETURN372:ENDIF373:*374:IF( EQUIL )THEN375:*376:* Compute row and column scalings to equilibrate the matrix A.377:*378:CALLZGBEQU( N, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, 379: $ AMAX, INFEQU ) 380:IF( INFEQU.EQ.0 )THEN381:*382:* Equilibrate the matrix.383:*384:CALLZLAQGB( N, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, 385: $ AMAX, EQUED ) 386: ROWEQU =LSAME( EQUED, 'R' ) .OR.LSAME( EQUED, 'B' ) 387: COLEQU =LSAME( EQUED, 'C' ) .OR.LSAME( EQUED, 'B' ) 388:ENDIF389:ENDIF390:*391:* Scale the right hand side.392:*393:IF( NOTRAN )THEN394:IF( ROWEQU )THEN395:DO40 J = 1, NRHS 396:DO30 I = 1, N 397:B( I, J ) =R( I )*B( I, J ) 398: 30CONTINUE399: 40CONTINUE400:ENDIF401:ELSEIF( COLEQU )THEN402:DO60 J = 1, NRHS 403:DO50 I = 1, N 404:B( I, J ) =C( I )*B( I, J ) 405: 50CONTINUE406: 60CONTINUE407:ENDIF408:*409:IF( NOFACT .OR. EQUIL )THEN410:*411:* Compute the LU factorization of the band matrix A.412:*413:DO70 J = 1, N 414: J1 =MAX( J-KU, 1 ) 415: J2 =MIN( J+KL, N ) 416:CALLZCOPY( J2-J1+1,AB( KU+1-J+J1, J ), 1, 417: $AFB( KL+KU+1-J+J1, J ), 1 ) 418: 70CONTINUE419:*420:CALLZGBTRF( N, N, KL, KU, AFB, LDAFB, IPIV, INFO ) 421:*422:* Return if INFO is non-zero.423:*424:IF( INFO.GT.0 )THEN425:*426:* Compute the reciprocal pivot growth factor of the427:* leading rank-deficient INFO columns of A.428:*429: ANORM = ZERO 430:DO90 J = 1, INFO 431:DO80 I =MAX( KU+2-J, 1 ),MIN( N+KU+1-J, KL+KU+1 ) 432: ANORM =MAX( ANORM,ABS(AB( I, J ) ) ) 433: 80CONTINUE434: 90CONTINUE435: RPVGRW =ZLANTB( 'M', 'U', 'N', INFO,MIN( INFO-1, KL+KU ), 436: $AFB(MAX( 1, KL+KU+2-INFO ), 1 ), LDAFB, 437: $ RWORK ) 438:IF( RPVGRW.EQ.ZERO )THEN439: RPVGRW = ONE 440:ELSE441: RPVGRW = ANORM / RPVGRW 442:ENDIF443:RWORK( 1 ) = RPVGRW 444: RCOND = ZERO 445:RETURN446:ENDIF447:ENDIF448:*449:* Compute the norm of the matrix A and the450:* reciprocal pivot growth factor RPVGRW.451:*452:IF( NOTRAN )THEN453: NORM = '1' 454:ELSE455: NORM = 'I' 456:ENDIF457: ANORM =ZLANGB( NORM, N, KL, KU, AB, LDAB, RWORK ) 458: RPVGRW =ZLANTB( 'M', 'U', 'N', N, KL+KU, AFB, LDAFB, RWORK ) 459:IF( RPVGRW.EQ.ZERO )THEN460: RPVGRW = ONE 461:ELSE462: RPVGRW =ZLANGB( 'M', N, KL, KU, AB, LDAB, RWORK ) / RPVGRW 463:ENDIF464:*465:* Compute the reciprocal of the condition number of A.466:*467:CALLZGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND, 468: $ WORK, RWORK, INFO ) 469:*470:* Compute the solution matrix X.471:*472:CALLZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 473:CALLZGBTRS( TRANS, N, KL, KU, NRHS, AFB, LDAFB, IPIV, X, LDX, 474: $ INFO ) 475:*476:* Use iterative refinement to improve the computed solution and477:* compute error bounds and backward error estimates for it.478:*479:CALLZGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, 480: $ B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO ) 481:*482:* Transform the solution matrix X to a solution of the original483:* system.484:*485:IF( NOTRAN )THEN486:IF( COLEQU )THEN487:DO110 J = 1, NRHS 488:DO100 I = 1, N 489:X( I, J ) =C( I )*X( I, J ) 490: 100CONTINUE491: 110CONTINUE492:DO120 J = 1, NRHS 493:FERR( J ) =FERR( J ) / COLCND 494: 120CONTINUE495:ENDIF496:ELSEIF( ROWEQU )THEN497:DO140 J = 1, NRHS 498:DO130 I = 1, N 499:X( I, J ) =R( I )*X( I, J ) 500: 130CONTINUE501: 140CONTINUE502:DO150 J = 1, NRHS 503:FERR( J ) =FERR( J ) / ROWCND 504: 150CONTINUE505:ENDIF506:*507:* Set INFO = N+1 if the matrix is singular to working precision.508:*509:IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 510: $ INFO = N + 1 511:*512:RWORK( 1 ) = RPVGRW 513:RETURN514:*515:* End of ZGBSVX516:*517:END518: