001:SUBROUTINESGBSVX( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, 002: $ LDAFB, IPIV, EQUED, R, C, B, LDB, X, LDX, 003: $ RCOND, FERR, BERR, WORK, IWORK, 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: REAL RCOND 013:* ..014:* .. Array Arguments ..015: INTEGERIPIV( * ),IWORK( * ) 016: REALAB( LDAB, * ),AFB( LDAFB, * ),B( LDB, * ), 017: $BERR( * ),C( * ),FERR( * ),R( * ), 018: $WORK( * ),X( LDX, * ) 019:* ..020:*021:* Purpose022:* =======023:*024:* SGBSVX uses the LU factorization to compute the solution to a real025:* system of linear equations A * X = B, A**T * X = B, or A**H * X = B,026:* where A is a band matrix of order N with KL subdiagonals and KU027:* superdiagonals, and X and B are N-by-NRHS matrices.028:*029:* Error bounds on the solution and a condition estimate are also030:* provided.031:*032:* Description033:* ===========034:*035:* The following steps are performed by this subroutine:036:*037:* 1. If FACT = 'E', real scaling factors are computed to equilibrate038:* the system:039:* TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B040:* TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B041:* TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B042:* Whether or not the system will be equilibrated depends on the043:* scaling of the matrix A, but if equilibration is used, A is044:* overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')045:* or diag(C)*B (if TRANS = 'T' or 'C').046:*047:* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the048:* matrix A (after equilibration if FACT = 'E') as049:* A = L * U,050:* where L is a product of permutation and unit lower triangular051:* matrices with KL subdiagonals, and U is upper triangular with052:* KL+KU superdiagonals.053:*054:* 3. If some U(i,i)=0, so that U is exactly singular, then the routine055:* returns with INFO = i. Otherwise, the factored form of A is used056:* to estimate the condition number of the matrix A. If the057:* reciprocal of the condition number is less than machine precision,058:* INFO = N+1 is returned as a warning, but the routine still goes on059:* to solve for X and compute error bounds as described below.060:*061:* 4. The system of equations is solved for X using the factored form062:* of A.063:*064:* 5. Iterative refinement is applied to improve the computed solution065:* matrix and calculate error bounds and backward error estimates066:* for it.067:*068:* 6. If equilibration was used, the matrix X is premultiplied by069:* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so070:* that it solves the original system before equilibration.071:*072:* Arguments073:* =========074:*075:* FACT (input) CHARACTER*1076:* Specifies whether or not the factored form of the matrix A is077:* supplied on entry, and if not, whether the matrix A should be078:* equilibrated before it is factored.079:* = 'F': On entry, AFB and IPIV contain the factored form of080:* A. If EQUED is not 'N', the matrix A has been081:* equilibrated with scaling factors given by R and C.082:* AB, AFB, and IPIV are not modified.083:* = 'N': The matrix A will be copied to AFB and factored.084:* = 'E': The matrix A will be equilibrated if necessary, then085:* copied to AFB and factored.086:*087:* TRANS (input) CHARACTER*1088:* Specifies the form of the system of equations.089:* = 'N': A * X = B (No transpose)090:* = 'T': A**T * X = B (Transpose)091:* = 'C': A**H * X = B (Transpose)092:*093:* N (input) INTEGER094:* The number of linear equations, i.e., the order of the095:* matrix A. N >= 0.096:*097:* KL (input) INTEGER098:* The number of subdiagonals within the band of A. KL >= 0.099:*100:* KU (input) INTEGER101:* The number of superdiagonals within the band of A. KU >= 0.102:*103:* NRHS (input) INTEGER104:* The number of right hand sides, i.e., the number of columns105:* of the matrices B and X. NRHS >= 0.106:*107:* AB (input/output) REAL array, dimension (LDAB,N)108:* On entry, the matrix A in band storage, in rows 1 to KL+KU+1.109:* The j-th column of A is stored in the j-th column of the110:* array AB as follows:111:* AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)112:*113:* If FACT = 'F' and EQUED is not 'N', then A must have been114:* equilibrated by the scaling factors in R and/or C. AB is not115:* modified if FACT = 'F' or 'N', or if FACT = 'E' and116:* EQUED = 'N' on exit.117:*118:* On exit, if EQUED .ne. 'N', A is scaled as follows:119:* EQUED = 'R': A := diag(R) * A120:* EQUED = 'C': A := A * diag(C)121:* EQUED = 'B': A := diag(R) * A * diag(C).122:*123:* LDAB (input) INTEGER124:* The leading dimension of the array AB. LDAB >= KL+KU+1.125:*126:* AFB (input or output) REAL array, dimension (LDAFB,N)127:* If FACT = 'F', then AFB is an input argument and on entry128:* contains details of the LU factorization of the band matrix129:* A, as computed by SGBTRF. U is stored as an upper triangular130:* band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,131:* and the multipliers used during the factorization are stored132:* in rows KL+KU+2 to 2*KL+KU+1. If EQUED .ne. 'N', then AFB is133:* the factored form of the equilibrated matrix A.134:*135:* If FACT = 'N', then AFB is an output argument and on exit136:* returns details of the LU factorization of A.137:*138:* If FACT = 'E', then AFB is an output argument and on exit139:* returns details of the LU factorization of the equilibrated140:* matrix A (see the description of AB for the form of the141:* equilibrated matrix).142:*143:* LDAFB (input) INTEGER144:* The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.145:*146:* IPIV (input or output) INTEGER array, dimension (N)147:* If FACT = 'F', then IPIV is an input argument and on entry148:* contains the pivot indices from the factorization A = L*U149:* as computed by SGBTRF; row i of the matrix was interchanged150:* with row IPIV(i).151:*152:* If FACT = 'N', then IPIV is an output argument and on exit153:* contains the pivot indices from the factorization A = L*U154:* of the original matrix A.155:*156:* If FACT = 'E', then IPIV is an output argument and on exit157:* contains the pivot indices from the factorization A = L*U158:* of the equilibrated matrix A.159:*160:* EQUED (input or output) CHARACTER*1161:* Specifies the form of equilibration that was done.162:* = 'N': No equilibration (always true if FACT = 'N').163:* = 'R': Row equilibration, i.e., A has been premultiplied by164:* diag(R).165:* = 'C': Column equilibration, i.e., A has been postmultiplied166:* by diag(C).167:* = 'B': Both row and column equilibration, i.e., A has been168:* replaced by diag(R) * A * diag(C).169:* EQUED is an input argument if FACT = 'F'; otherwise, it is an170:* output argument.171:*172:* R (input or output) REAL array, dimension (N)173:* The row scale factors for A. If EQUED = 'R' or 'B', A is174:* multiplied on the left by diag(R); if EQUED = 'N' or 'C', R175:* is not accessed. R is an input argument if FACT = 'F';176:* otherwise, R is an output argument. If FACT = 'F' and177:* EQUED = 'R' or 'B', each element of R must be positive.178:*179:* C (input or output) REAL array, dimension (N)180:* The column scale factors for A. If EQUED = 'C' or 'B', A is181:* multiplied on the right by diag(C); if EQUED = 'N' or 'R', C182:* is not accessed. C is an input argument if FACT = 'F';183:* otherwise, C is an output argument. If FACT = 'F' and184:* EQUED = 'C' or 'B', each element of C must be positive.185:*186:* B (input/output) REAL array, dimension (LDB,NRHS)187:* On entry, the right hand side matrix B.188:* On exit,189:* if EQUED = 'N', B is not modified;190:* if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by191:* diag(R)*B;192:* if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is193:* overwritten by diag(C)*B.194:*195:* LDB (input) INTEGER196:* The leading dimension of the array B. LDB >= max(1,N).197:*198:* X (output) REAL array, dimension (LDX,NRHS)199:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X200:* to the original system of equations. Note that A and B are201:* modified on exit if EQUED .ne. 'N', and the solution to the202:* equilibrated system is inv(diag(C))*X if TRANS = 'N' and203:* EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'204:* and EQUED = 'R' or 'B'.205:*206:* LDX (input) INTEGER207:* The leading dimension of the array X. LDX >= max(1,N).208:*209:* RCOND (output) REAL210:* The estimate of the reciprocal condition number of the matrix211:* A after equilibration (if done). If RCOND is less than the212:* machine precision (in particular, if RCOND = 0), the matrix213:* is singular to working precision. This condition is214:* indicated by a return code of INFO > 0.215:*216:* FERR (output) REAL array, dimension (NRHS)217:* The estimated forward error bound for each solution vector218:* X(j) (the j-th column of the solution matrix X).219:* If XTRUE is the true solution corresponding to X(j), FERR(j)220:* is an estimated upper bound for the magnitude of the largest221:* element in (X(j) - XTRUE) divided by the magnitude of the222:* largest element in X(j). The estimate is as reliable as223:* the estimate for RCOND, and is almost always a slight224:* overestimate of the true error.225:*226:* BERR (output) REAL array, dimension (NRHS)227:* The componentwise relative backward error of each solution228:* vector X(j) (i.e., the smallest relative change in229:* any element of A or B that makes X(j) an exact solution).230:*231:* WORK (workspace/output) REAL array, dimension (3*N)232:* On exit, WORK(1) contains the reciprocal pivot growth233:* factor norm(A)/norm(U). The "max absolute element" norm is234:* used. If WORK(1) is much less than 1, then the stability235:* of the LU factorization of the (equilibrated) matrix A236:* could be poor. This also means that the solution X, condition237:* estimator RCOND, and forward error bound FERR could be238:* unreliable. If factorization fails with 0<INFO<=N, then239:* WORK(1) contains the reciprocal pivot growth factor for the240:* leading INFO columns of A.241:*242:* IWORK (workspace) INTEGER array, dimension (N)243:*244:* INFO (output) INTEGER245:* = 0: successful exit246:* < 0: if INFO = -i, the i-th argument had an illegal value247:* > 0: if INFO = i, and i is248:* <= N: U(i,i) is exactly zero. The factorization249:* has been completed, but the factor U is exactly250:* singular, so the solution and error bounds251:* could not be computed. RCOND = 0 is returned.252:* = N+1: U is nonsingular, but RCOND is less than machine253:* precision, meaning that the matrix is singular254:* to working precision. Nevertheless, the255:* solution and error bounds are computed because256:* there are a number of situations where the257:* computed solution can be more accurate than the258:*259:* value of RCOND would suggest.260:* =====================================================================261:* Moved setting of INFO = N+1 so INFO does not subsequently get262:* overwritten. Sven, 17 Mar 05.263:* =====================================================================264:*265:* .. Parameters ..266: REAL ZERO, ONE 267:PARAMETER( ZERO = 0.0E+0, ONE = 1.0E+0 ) 268:* ..269:* .. Local Scalars ..270:LOGICALCOLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU 271: CHARACTER NORM 272: INTEGER I, INFEQU, J, J1, J2 273: REAL AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN, 274: $ ROWCND, RPVGRW, SMLNUM 275:* ..276:* .. External Functions ..277:LOGICALLSAME 278: REAL SLAMCH, SLANGB, SLANTB 279:EXTERNALLSAME, SLAMCH, SLANGB, SLANTB 280:* ..281:* .. External Subroutines ..282:EXTERNALSCOPY, SGBCON, SGBEQU, SGBRFS, SGBTRF, SGBTRS, 283: $ SLACPY, SLAQGB, XERBLA 284:* ..285:* .. Intrinsic Functions ..286:INTRINSICABS, MAX, MIN 287:* ..288:* .. Executable Statements ..289:*290: INFO = 0 291: NOFACT =LSAME( FACT, 'N' ) 292: EQUIL =LSAME( FACT, 'E' ) 293: NOTRAN =LSAME( TRANS, 'N' ) 294:IF( NOFACT .OR. EQUIL )THEN295: EQUED = 'N' 296: ROWEQU = .FALSE. 297: COLEQU = .FALSE. 298:ELSE299: ROWEQU =LSAME( EQUED, 'R' ) .OR.LSAME( EQUED, 'B' ) 300: COLEQU =LSAME( EQUED, 'C' ) .OR.LSAME( EQUED, 'B' ) 301: SMLNUM =SLAMCH( 'Safe minimum' ) 302: BIGNUM = ONE / SMLNUM 303:ENDIF304:*305:* Test the input parameters.306:*307:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 308: $THEN309: INFO = -1 310:ELSEIF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 311: $LSAME( TRANS, 'C' ) )THEN312: INFO = -2 313:ELSEIF( N.LT.0 )THEN314: INFO = -3 315:ELSEIF( KL.LT.0 )THEN316: INFO = -4 317:ELSEIF( KU.LT.0 )THEN318: INFO = -5 319:ELSEIF( NRHS.LT.0 )THEN320: INFO = -6 321:ELSEIF( LDAB.LT.KL+KU+1 )THEN322: INFO = -8 323:ELSEIF( LDAFB.LT.2*KL+KU+1 )THEN324: INFO = -10 325:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 326: $ ( ROWEQU .OR. COLEQU .OR.LSAME( EQUED, 'N' ) ) )THEN327: INFO = -12 328:ELSE329:IF( ROWEQU )THEN330: RCMIN = BIGNUM 331: RCMAX = ZERO 332:DO10 J = 1, N 333: RCMIN =MIN( RCMIN,R( J ) ) 334: RCMAX =MAX( RCMAX,R( J ) ) 335: 10CONTINUE336:IF( RCMIN.LE.ZERO )THEN337: INFO = -13 338:ELSEIF( N.GT.0 )THEN339: ROWCND =MAX( RCMIN, SMLNUM ) /MIN( RCMAX, BIGNUM ) 340:ELSE341: ROWCND = ONE 342:ENDIF343:ENDIF344:IF( COLEQU .AND. INFO.EQ.0 )THEN345: RCMIN = BIGNUM 346: RCMAX = ZERO 347:DO20 J = 1, N 348: RCMIN =MIN( RCMIN,C( J ) ) 349: RCMAX =MAX( RCMAX,C( J ) ) 350: 20CONTINUE351:IF( RCMIN.LE.ZERO )THEN352: INFO = -14 353:ELSEIF( N.GT.0 )THEN354: COLCND =MAX( RCMIN, SMLNUM ) /MIN( RCMAX, BIGNUM ) 355:ELSE356: COLCND = ONE 357:ENDIF358:ENDIF359:IF( INFO.EQ.0 )THEN360:IF( LDB.LT.MAX( 1, N ) )THEN361: INFO = -16 362:ELSEIF( LDX.LT.MAX( 1, N ) )THEN363: INFO = -18 364:ENDIF365:ENDIF366:ENDIF367:*368:IF( INFO.NE.0 )THEN369:CALLXERBLA( 'SGBSVX', -INFO ) 370:RETURN371:ENDIF372:*373:IF( EQUIL )THEN374:*375:* Compute row and column scalings to equilibrate the matrix A.376:*377:CALLSGBEQU( N, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, 378: $ AMAX, INFEQU ) 379:IF( INFEQU.EQ.0 )THEN380:*381:* Equilibrate the matrix.382:*383:CALLSLAQGB( N, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, 384: $ AMAX, EQUED ) 385: ROWEQU =LSAME( EQUED, 'R' ) .OR.LSAME( EQUED, 'B' ) 386: COLEQU =LSAME( EQUED, 'C' ) .OR.LSAME( EQUED, 'B' ) 387:ENDIF388:ENDIF389:*390:* Scale the right hand side.391:*392:IF( NOTRAN )THEN393:IF( ROWEQU )THEN394:DO40 J = 1, NRHS 395:DO30 I = 1, N 396:B( I, J ) =R( I )*B( I, J ) 397: 30CONTINUE398: 40CONTINUE399:ENDIF400:ELSEIF( COLEQU )THEN401:DO60 J = 1, NRHS 402:DO50 I = 1, N 403:B( I, J ) =C( I )*B( I, J ) 404: 50CONTINUE405: 60CONTINUE406:ENDIF407:*408:IF( NOFACT .OR. EQUIL )THEN409:*410:* Compute the LU factorization of the band matrix A.411:*412:DO70 J = 1, N 413: J1 =MAX( J-KU, 1 ) 414: J2 =MIN( J+KL, N ) 415:CALLSCOPY( J2-J1+1,AB( KU+1-J+J1, J ), 1, 416: $AFB( KL+KU+1-J+J1, J ), 1 ) 417: 70CONTINUE418:*419:CALLSGBTRF( N, N, KL, KU, AFB, LDAFB, IPIV, INFO ) 420:*421:* Return if INFO is non-zero.422:*423:IF( INFO.GT.0 )THEN424:*425:* Compute the reciprocal pivot growth factor of the426:* leading rank-deficient INFO columns of A.427:*428: ANORM = ZERO 429:DO90 J = 1, INFO 430:DO80 I =MAX( KU+2-J, 1 ),MIN( N+KU+1-J, KL+KU+1 ) 431: ANORM =MAX( ANORM,ABS(AB( I, J ) ) ) 432: 80CONTINUE433: 90CONTINUE434: RPVGRW =SLANTB( 'M', 'U', 'N', INFO,MIN( INFO-1, KL+KU ), 435: $AFB(MAX( 1, KL+KU+2-INFO ), 1 ), LDAFB, 436: $ WORK ) 437:IF( RPVGRW.EQ.ZERO )THEN438: RPVGRW = ONE 439:ELSE440: RPVGRW = ANORM / RPVGRW 441:ENDIF442:WORK( 1 ) = RPVGRW 443: RCOND = ZERO 444:RETURN445:ENDIF446:ENDIF447:*448:* Compute the norm of the matrix A and the449:* reciprocal pivot growth factor RPVGRW.450:*451:IF( NOTRAN )THEN452: NORM = '1' 453:ELSE454: NORM = 'I' 455:ENDIF456: ANORM =SLANGB( NORM, N, KL, KU, AB, LDAB, WORK ) 457: RPVGRW =SLANTB( 'M', 'U', 'N', N, KL+KU, AFB, LDAFB, WORK ) 458:IF( RPVGRW.EQ.ZERO )THEN459: RPVGRW = ONE 460:ELSE461: RPVGRW =SLANGB( 'M', N, KL, KU, AB, LDAB, WORK ) / RPVGRW 462:ENDIF463:*464:* Compute the reciprocal of the condition number of A.465:*466:CALLSGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND, 467: $ WORK, IWORK, INFO ) 468:*469:* Compute the solution matrix X.470:*471:CALLSLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 472:CALLSGBTRS( TRANS, N, KL, KU, NRHS, AFB, LDAFB, IPIV, X, LDX, 473: $ INFO ) 474:*475:* Use iterative refinement to improve the computed solution and476:* compute error bounds and backward error estimates for it.477:*478:CALLSGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, 479: $ B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO ) 480:*481:* Transform the solution matrix X to a solution of the original482:* system.483:*484:IF( NOTRAN )THEN485:IF( COLEQU )THEN486:DO110 J = 1, NRHS 487:DO100 I = 1, N 488:X( I, J ) =C( I )*X( I, J ) 489: 100CONTINUE490: 110CONTINUE491:DO120 J = 1, NRHS 492:FERR( J ) =FERR( J ) / COLCND 493: 120CONTINUE494:ENDIF495:ELSEIF( ROWEQU )THEN496:DO140 J = 1, NRHS 497:DO130 I = 1, N 498:X( I, J ) =R( I )*X( I, J ) 499: 130CONTINUE500: 140CONTINUE501:DO150 J = 1, NRHS 502:FERR( J ) =FERR( J ) / ROWCND 503: 150CONTINUE504:ENDIF505:*506:* Set INFO = N+1 if the matrix is singular to working precision.507:*508:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 509: $ INFO = N + 1 510:*511:WORK( 1 ) = RPVGRW 512:RETURN513:*514:* End of SGBSVX515:*516:END517: