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:* -- LAPACK is a software package provided by Univ. of Tennessee, --007:* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--008:* November 2006009:*010:* .. Scalar Arguments ..011: CHARACTER EQUED, FACT, TRANS 012: INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS 013: DOUBLE PRECISION RCOND 014:* ..015:* .. Array Arguments ..016: INTEGERIPIV( * ) 017: DOUBLE PRECISIONBERR( * ),C( * ),FERR( * ),R( * ), 018: $RWORK( * ) 019: COMPLEX*16AB( LDAB, * ),AFB( LDAFB, * ),B( LDB, * ), 020: $WORK( * ),X( LDX, * ) 021:* ..022:*023:* Purpose024:* =======025:*026:* ZGBSVX uses the LU factorization to compute the solution to a complex027:* system of linear equations A * X = B, A**T * X = B, or A**H * X = B,028:* where A is a band matrix of order N with KL subdiagonals and KU029:* superdiagonals, and X and B are N-by-NRHS matrices.030:*031:* Error bounds on the solution and a condition estimate are also032:* provided.033:*034:* Description035:* ===========036:*037:* The following steps are performed by this subroutine:038:*039:* 1. If FACT = 'E', real scaling factors are computed to equilibrate040:* the system:041:* TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B042:* TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B043:* TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B044:* Whether or not the system will be equilibrated depends on the045:* scaling of the matrix A, but if equilibration is used, A is046:* overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')047:* or diag(C)*B (if TRANS = 'T' or 'C').048:*049:* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the050:* matrix A (after equilibration if FACT = 'E') as051:* A = L * U,052:* where L is a product of permutation and unit lower triangular053:* matrices with KL subdiagonals, and U is upper triangular with054:* KL+KU superdiagonals.055:*056:* 3. If some U(i,i)=0, so that U is exactly singular, then the routine057:* returns with INFO = i. Otherwise, the factored form of A is used058:* to estimate the condition number of the matrix A. If the059:* reciprocal of the condition number is less than machine precision,060:* INFO = N+1 is returned as a warning, but the routine still goes on061:* to solve for X and compute error bounds as described below.062:*063:* 4. The system of equations is solved for X using the factored form064:* of A.065:*066:* 5. Iterative refinement is applied to improve the computed solution067:* matrix and calculate error bounds and backward error estimates068:* for it.069:*070:* 6. If equilibration was used, the matrix X is premultiplied by071:* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so072:* that it solves the original system before equilibration.073:*074:* Arguments075:* =========076:*077:* FACT (input) CHARACTER*1078:* Specifies whether or not the factored form of the matrix A is079:* supplied on entry, and if not, whether the matrix A should be080:* equilibrated before it is factored.081:* = 'F': On entry, AFB and IPIV contain the factored form of082:* A. If EQUED is not 'N', the matrix A has been083:* equilibrated with scaling factors given by R and C.084:* AB, AFB, and IPIV are not modified.085:* = 'N': The matrix A will be copied to AFB and factored.086:* = 'E': The matrix A will be equilibrated if necessary, then087:* copied to AFB and factored.088:*089:* TRANS (input) CHARACTER*1090:* Specifies the form of the system of equations.091:* = 'N': A * X = B (No transpose)092:* = 'T': A**T * X = B (Transpose)093:* = 'C': A**H * X = B (Conjugate transpose)094:*095:* N (input) INTEGER096:* The number of linear equations, i.e., the order of the097:* matrix A. N >= 0.098:*099:* KL (input) INTEGER100:* The number of subdiagonals within the band of A. KL >= 0.101:*102:* KU (input) INTEGER103:* The number of superdiagonals within the band of A. KU >= 0.104:*105:* NRHS (input) INTEGER106:* The number of right hand sides, i.e., the number of columns107:* of the matrices B and X. NRHS >= 0.108:*109:* AB (input/output) COMPLEX*16 array, dimension (LDAB,N)110:* On entry, the matrix A in band storage, in rows 1 to KL+KU+1.111:* The j-th column of A is stored in the j-th column of the112:* array AB as follows:113:* AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)114:*115:* If FACT = 'F' and EQUED is not 'N', then A must have been116:* equilibrated by the scaling factors in R and/or C. AB is not117:* modified if FACT = 'F' or 'N', or if FACT = 'E' and118:* EQUED = 'N' on exit.119:*120:* On exit, if EQUED .ne. 'N', A is scaled as follows:121:* EQUED = 'R': A := diag(R) * A122:* EQUED = 'C': A := A * diag(C)123:* EQUED = 'B': A := diag(R) * A * diag(C).124:*125:* LDAB (input) INTEGER126:* The leading dimension of the array AB. LDAB >= KL+KU+1.127:*128:* AFB (input or output) COMPLEX*16 array, dimension (LDAFB,N)129:* If FACT = 'F', then AFB is an input argument and on entry130:* contains details of the LU factorization of the band matrix131:* A, as computed by ZGBTRF. U is stored as an upper triangular132:* band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,133:* and the multipliers used during the factorization are stored134:* in rows KL+KU+2 to 2*KL+KU+1. If EQUED .ne. 'N', then AFB is135:* the factored form of the equilibrated matrix A.136:*137:* If FACT = 'N', then AFB is an output argument and on exit138:* returns details of the LU factorization of A.139:*140:* If FACT = 'E', then AFB is an output argument and on exit141:* returns details of the LU factorization of the equilibrated142:* matrix A (see the description of AB for the form of the143:* equilibrated matrix).144:*145:* LDAFB (input) INTEGER146:* The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.147:*148:* IPIV (input or output) INTEGER array, dimension (N)149:* If FACT = 'F', then IPIV is an input argument and on entry150:* contains the pivot indices from the factorization A = L*U151:* as computed by ZGBTRF; row i of the matrix was interchanged152:* with row IPIV(i).153:*154:* If FACT = 'N', then IPIV is an output argument and on exit155:* contains the pivot indices from the factorization A = L*U156:* of the original matrix A.157:*158:* If FACT = 'E', then IPIV is an output argument and on exit159:* contains the pivot indices from the factorization A = L*U160:* of the equilibrated matrix A.161:*162:* EQUED (input or output) CHARACTER*1163:* Specifies the form of equilibration that was done.164:* = 'N': No equilibration (always true if FACT = 'N').165:* = 'R': Row equilibration, i.e., A has been premultiplied by166:* diag(R).167:* = 'C': Column equilibration, i.e., A has been postmultiplied168:* by diag(C).169:* = 'B': Both row and column equilibration, i.e., A has been170:* replaced by diag(R) * A * diag(C).171:* EQUED is an input argument if FACT = 'F'; otherwise, it is an172:* output argument.173:*174:* R (input or output) DOUBLE PRECISION array, dimension (N)175:* The row scale factors for A. If EQUED = 'R' or 'B', A is176:* multiplied on the left by diag(R); if EQUED = 'N' or 'C', R177:* is not accessed. R is an input argument if FACT = 'F';178:* otherwise, R is an output argument. If FACT = 'F' and179:* EQUED = 'R' or 'B', each element of R must be positive.180:*181:* C (input or output) DOUBLE PRECISION array, dimension (N)182:* The column scale factors for A. If EQUED = 'C' or 'B', A is183:* multiplied on the right by diag(C); if EQUED = 'N' or 'R', C184:* is not accessed. C is an input argument if FACT = 'F';185:* otherwise, C is an output argument. If FACT = 'F' and186:* EQUED = 'C' or 'B', each element of C must be positive.187:*188:* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)189:* On entry, the right hand side matrix B.190:* On exit,191:* if EQUED = 'N', B is not modified;192:* if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by193:* diag(R)*B;194:* if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is195:* overwritten by diag(C)*B.196:*197:* LDB (input) INTEGER198:* The leading dimension of the array B. LDB >= max(1,N).199:*200:* X (output) COMPLEX*16 array, dimension (LDX,NRHS)201:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X202:* to the original system of equations. Note that A and B are203:* modified on exit if EQUED .ne. 'N', and the solution to the204:* equilibrated system is inv(diag(C))*X if TRANS = 'N' and205:* EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'206:* and EQUED = 'R' or 'B'.207:*208:* LDX (input) INTEGER209:* The leading dimension of the array X. LDX >= max(1,N).210:*211:* RCOND (output) DOUBLE PRECISION212:* The estimate of the reciprocal condition number of the matrix213:* A after equilibration (if done). If RCOND is less than the214:* machine precision (in particular, if RCOND = 0), the matrix215:* is singular to working precision. This condition is216:* indicated by a return code of INFO > 0.217:*218:* FERR (output) DOUBLE PRECISION array, dimension (NRHS)219:* The estimated forward error bound for each solution vector220:* X(j) (the j-th column of the solution matrix X).221:* If XTRUE is the true solution corresponding to X(j), FERR(j)222:* is an estimated upper bound for the magnitude of the largest223:* element in (X(j) - XTRUE) divided by the magnitude of the224:* largest element in X(j). The estimate is as reliable as225:* the estimate for RCOND, and is almost always a slight226:* overestimate of the true error.227:*228:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)229:* The componentwise relative backward error of each solution230:* vector X(j) (i.e., the smallest relative change in231:* any element of A or B that makes X(j) an exact solution).232:*233:* WORK (workspace) COMPLEX*16 array, dimension (2*N)234:*235:* RWORK (workspace/output) DOUBLE PRECISION array, dimension (N)236:* On exit, RWORK(1) contains the reciprocal pivot growth237:* factor norm(A)/norm(U). The "max absolute element" norm is238:* used. If RWORK(1) is much less than 1, then the stability239:* of the LU factorization of the (equilibrated) matrix A240:* could be poor. This also means that the solution X, condition241:* estimator RCOND, and forward error bound FERR could be242:* unreliable. If factorization fails with 0<INFO<=N, then243:* RWORK(1) contains the reciprocal pivot growth factor for the244:* leading INFO columns of A.245:*246:* INFO (output) INTEGER247:* = 0: successful exit248:* < 0: if INFO = -i, the i-th argument had an illegal value249:* > 0: if INFO = i, and i is250:* <= N: U(i,i) is exactly zero. The factorization251:* has been completed, but the factor U is exactly252:* singular, so the solution and error bounds253:* could not be computed. RCOND = 0 is returned.254:* = N+1: U is nonsingular, but RCOND is less than machine255:* precision, meaning that the matrix is singular256:* to working precision. Nevertheless, the257:* solution and error bounds are computed because258:* there are a number of situations where the259:* computed solution can be more accurate than the260:* value of RCOND would suggest.261:*262:* =====================================================================263:* Moved setting of INFO = N+1 so INFO does not subsequently get264:* overwritten. Sven, 17 Mar 05.265:* =====================================================================266:*267:* .. Parameters ..268: DOUBLE PRECISION ZERO, ONE 269:PARAMETER( ZERO = 0.0D+0, ONE = 1.0D+0 ) 270:* ..271:* .. Local Scalars ..272:LOGICALCOLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU 273: CHARACTER NORM 274: INTEGER I, INFEQU, J, J1, J2 275: DOUBLE PRECISION AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN, 276: $ ROWCND, RPVGRW, SMLNUM 277:* ..278:* .. External Functions ..279:LOGICALLSAME 280: DOUBLE PRECISION DLAMCH, ZLANGB, ZLANTB 281:EXTERNALLSAME, DLAMCH, ZLANGB, ZLANTB 282:* ..283:* .. External Subroutines ..284:EXTERNALXERBLA, ZCOPY, ZGBCON, ZGBEQU, ZGBRFS, ZGBTRF, 285: $ ZGBTRS, ZLACPY, ZLAQGB 286:* ..287:* .. Intrinsic Functions ..288:INTRINSICABS, MAX, MIN 289:* ..290:* .. Executable Statements ..291:*292: INFO = 0 293: NOFACT =LSAME( FACT, 'N' ) 294: EQUIL =LSAME( FACT, 'E' ) 295: NOTRAN =LSAME( TRANS, 'N' ) 296:IF( NOFACT .OR. EQUIL )THEN297: EQUED = 'N' 298: ROWEQU = .FALSE. 299: COLEQU = .FALSE. 300:ELSE301: ROWEQU =LSAME( EQUED, 'R' ) .OR.LSAME( EQUED, 'B' ) 302: COLEQU =LSAME( EQUED, 'C' ) .OR.LSAME( EQUED, 'B' ) 303: SMLNUM =DLAMCH( 'Safe minimum' ) 304: BIGNUM = ONE / SMLNUM 305:ENDIF306:*307:* Test the input parameters.308:*309:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 310: $THEN311: INFO = -1 312:ELSEIF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 313: $LSAME( TRANS, 'C' ) )THEN314: INFO = -2 315:ELSEIF( N.LT.0 )THEN316: INFO = -3 317:ELSEIF( KL.LT.0 )THEN318: INFO = -4 319:ELSEIF( KU.LT.0 )THEN320: INFO = -5 321:ELSEIF( NRHS.LT.0 )THEN322: INFO = -6 323:ELSEIF( LDAB.LT.KL+KU+1 )THEN324: INFO = -8 325:ELSEIF( LDAFB.LT.2*KL+KU+1 )THEN326: INFO = -10 327:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 328: $ ( ROWEQU .OR. COLEQU .OR.LSAME( EQUED, 'N' ) ) )THEN329: INFO = -12 330:ELSE331:IF( ROWEQU )THEN332: RCMIN = BIGNUM 333: RCMAX = ZERO 334:DO10 J = 1, N 335: RCMIN =MIN( RCMIN,R( J ) ) 336: RCMAX =MAX( RCMAX,R( J ) ) 337: 10CONTINUE338:IF( RCMIN.LE.ZERO )THEN339: INFO = -13 340:ELSEIF( N.GT.0 )THEN341: ROWCND =MAX( RCMIN, SMLNUM ) /MIN( RCMAX, BIGNUM ) 342:ELSE343: ROWCND = ONE 344:ENDIF345:ENDIF346:IF( COLEQU .AND. INFO.EQ.0 )THEN347: RCMIN = BIGNUM 348: RCMAX = ZERO 349:DO20 J = 1, N 350: RCMIN =MIN( RCMIN,C( J ) ) 351: RCMAX =MAX( RCMAX,C( J ) ) 352: 20CONTINUE353:IF( RCMIN.LE.ZERO )THEN354: INFO = -14 355:ELSEIF( N.GT.0 )THEN356: COLCND =MAX( RCMIN, SMLNUM ) /MIN( RCMAX, BIGNUM ) 357:ELSE358: COLCND = ONE 359:ENDIF360:ENDIF361:IF( INFO.EQ.0 )THEN362:IF( LDB.LT.MAX( 1, N ) )THEN363: INFO = -16 364:ELSEIF( LDX.LT.MAX( 1, N ) )THEN365: INFO = -18 366:ENDIF367:ENDIF368:ENDIF369:*370:IF( INFO.NE.0 )THEN371:CALLXERBLA( 'ZGBSVX', -INFO ) 372:RETURN373:ENDIF374:*375:IF( EQUIL )THEN376:*377:* Compute row and column scalings to equilibrate the matrix A.378:*379:CALLZGBEQU( N, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, 380: $ AMAX, INFEQU ) 381:IF( INFEQU.EQ.0 )THEN382:*383:* Equilibrate the matrix.384:*385:CALLZLAQGB( N, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, 386: $ AMAX, EQUED ) 387: ROWEQU =LSAME( EQUED, 'R' ) .OR.LSAME( EQUED, 'B' ) 388: COLEQU =LSAME( EQUED, 'C' ) .OR.LSAME( EQUED, 'B' ) 389:ENDIF390:ENDIF391:*392:* Scale the right hand side.393:*394:IF( NOTRAN )THEN395:IF( ROWEQU )THEN396:DO40 J = 1, NRHS 397:DO30 I = 1, N 398:B( I, J ) =R( I )*B( I, J ) 399: 30CONTINUE400: 40CONTINUE401:ENDIF402:ELSEIF( COLEQU )THEN403:DO60 J = 1, NRHS 404:DO50 I = 1, N 405:B( I, J ) =C( I )*B( I, J ) 406: 50CONTINUE407: 60CONTINUE408:ENDIF409:*410:IF( NOFACT .OR. EQUIL )THEN411:*412:* Compute the LU factorization of the band matrix A.413:*414:DO70 J = 1, N 415: J1 =MAX( J-KU, 1 ) 416: J2 =MIN( J+KL, N ) 417:CALLZCOPY( J2-J1+1,AB( KU+1-J+J1, J ), 1, 418: $AFB( KL+KU+1-J+J1, J ), 1 ) 419: 70CONTINUE420:*421:CALLZGBTRF( N, N, KL, KU, AFB, LDAFB, IPIV, INFO ) 422:*423:* Return if INFO is non-zero.424:*425:IF( INFO.GT.0 )THEN426:*427:* Compute the reciprocal pivot growth factor of the428:* leading rank-deficient INFO columns of A.429:*430: ANORM = ZERO 431:DO90 J = 1, INFO 432:DO80 I =MAX( KU+2-J, 1 ),MIN( N+KU+1-J, KL+KU+1 ) 433: ANORM =MAX( ANORM,ABS(AB( I, J ) ) ) 434: 80CONTINUE435: 90CONTINUE436: RPVGRW =ZLANTB( 'M', 'U', 'N', INFO,MIN( INFO-1, KL+KU ), 437: $AFB(MAX( 1, KL+KU+2-INFO ), 1 ), LDAFB, 438: $ RWORK ) 439:IF( RPVGRW.EQ.ZERO )THEN440: RPVGRW = ONE 441:ELSE442: RPVGRW = ANORM / RPVGRW 443:ENDIF444:RWORK( 1 ) = RPVGRW 445: RCOND = ZERO 446:RETURN447:ENDIF448:ENDIF449:*450:* Compute the norm of the matrix A and the451:* reciprocal pivot growth factor RPVGRW.452:*453:IF( NOTRAN )THEN454: NORM = '1' 455:ELSE456: NORM = 'I' 457:ENDIF458: ANORM =ZLANGB( NORM, N, KL, KU, AB, LDAB, RWORK ) 459: RPVGRW =ZLANTB( 'M', 'U', 'N', N, KL+KU, AFB, LDAFB, RWORK ) 460:IF( RPVGRW.EQ.ZERO )THEN461: RPVGRW = ONE 462:ELSE463: RPVGRW =ZLANGB( 'M', N, KL, KU, AB, LDAB, RWORK ) / RPVGRW 464:ENDIF465:*466:* Compute the reciprocal of the condition number of A.467:*468:CALLZGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND, 469: $ WORK, RWORK, INFO ) 470:*471:* Compute the solution matrix X.472:*473:CALLZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 474:CALLZGBTRS( TRANS, N, KL, KU, NRHS, AFB, LDAFB, IPIV, X, LDX, 475: $ INFO ) 476:*477:* Use iterative refinement to improve the computed solution and478:* compute error bounds and backward error estimates for it.479:*480:CALLZGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, 481: $ B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO ) 482:*483:* Transform the solution matrix X to a solution of the original484:* system.485:*486:IF( NOTRAN )THEN487:IF( COLEQU )THEN488:DO110 J = 1, NRHS 489:DO100 I = 1, N 490:X( I, J ) =C( I )*X( I, J ) 491: 100CONTINUE492: 110CONTINUE493:DO120 J = 1, NRHS 494:FERR( J ) =FERR( J ) / COLCND 495: 120CONTINUE496:ENDIF497:ELSEIF( ROWEQU )THEN498:DO140 J = 1, NRHS 499:DO130 I = 1, N 500:X( I, J ) =R( I )*X( I, J ) 501: 130CONTINUE502: 140CONTINUE503:DO150 J = 1, NRHS 504:FERR( J ) =FERR( J ) / ROWCND 505: 150CONTINUE506:ENDIF507:*508:* Set INFO = N+1 if the matrix is singular to working precision.509:*510:IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 511: $ INFO = N + 1 512:*513:RWORK( 1 ) = RPVGRW 514:RETURN515:*516:* End of ZGBSVX517:*518:END519: