001:SUBROUTINECPBSVX( FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, 002: $ EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR, 003: $ 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, UPLO 011: INTEGER INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS 012: REAL RCOND 013:* ..014:* .. Array Arguments ..015: REALBERR( * ),FERR( * ),RWORK( * ),S( * ) 016: COMPLEXAB( LDAB, * ),AFB( LDAFB, * ),B( LDB, * ), 017: $WORK( * ),X( LDX, * ) 018:* ..019:*020:* Purpose021:* =======022:*023:* CPBSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to024:* compute the solution to a complex system of linear equations025:* A * X = B,026:* where A is an N-by-N Hermitian positive definite band matrix and X027:* 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:036:*037:* 1. If FACT = 'E', real scaling factors are computed to equilibrate038:* the system:039:* diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B040:* Whether or not the system will be equilibrated depends on the041:* scaling of the matrix A, but if equilibration is used, A is042:* overwritten by diag(S)*A*diag(S) and B by diag(S)*B.043:*044:* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to045:* factor the matrix A (after equilibration if FACT = 'E') as046:* A = U**H * U, if UPLO = 'U', or047:* A = L * L**H, if UPLO = 'L',048:* where U is an upper triangular band matrix, and L is a lower049:* triangular band matrix.050:*051:* 3. If the leading i-by-i principal minor is not positive definite,052:* then the routine returns with INFO = i. Otherwise, the factored053:* form of A is used to estimate the condition number of the matrix054:* A. If the reciprocal of the condition number is less than machine055:* precision, INFO = N+1 is returned as a warning, but the routine056:* still goes on to solve for X and compute error bounds as057:* described below.058:*059:* 4. The system of equations is solved for X using the factored form060:* of A.061:*062:* 5. Iterative refinement is applied to improve the computed solution063:* matrix and calculate error bounds and backward error estimates064:* for it.065:*066:* 6. If equilibration was used, the matrix X is premultiplied by067:* diag(S) so that it solves the original system before068:* equilibration.069:*070:* Arguments071:* =========072:*073:* FACT (input) CHARACTER*1074:* Specifies whether or not the factored form of the matrix A is075:* supplied on entry, and if not, whether the matrix A should be076:* equilibrated before it is factored.077:* = 'F': On entry, AFB contains the factored form of A.078:* If EQUED = 'Y', the matrix A has been equilibrated079:* with scaling factors given by S. AB and AFB will not080:* be modified.081:* = 'N': The matrix A will be copied to AFB and factored.082:* = 'E': The matrix A will be equilibrated if necessary, then083:* copied to AFB and factored.084:*085:* UPLO (input) CHARACTER*1086:* = 'U': Upper triangle of A is stored;087:* = 'L': Lower triangle of A is stored.088:*089:* N (input) INTEGER090:* The number of linear equations, i.e., the order of the091:* matrix A. N >= 0.092:*093:* KD (input) INTEGER094:* The number of superdiagonals of the matrix A if UPLO = 'U',095:* or the number of subdiagonals if UPLO = 'L'. KD >= 0.096:*097:* NRHS (input) INTEGER098:* The number of right-hand sides, i.e., the number of columns099:* of the matrices B and X. NRHS >= 0.100:*101:* AB (input/output) COMPLEX array, dimension (LDAB,N)102:* On entry, the upper or lower triangle of the Hermitian band103:* matrix A, stored in the first KD+1 rows of the array, except104:* if FACT = 'F' and EQUED = 'Y', then A must contain the105:* equilibrated matrix diag(S)*A*diag(S). The j-th column of A106:* is stored in the j-th column of the array AB as follows:107:* if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;108:* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(N,j+KD).109:* See below for further details.110:*111:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by112:* diag(S)*A*diag(S).113:*114:* LDAB (input) INTEGER115:* The leading dimension of the array A. LDAB >= KD+1.116:*117:* AFB (input or output) COMPLEX array, dimension (LDAFB,N)118:* If FACT = 'F', then AFB is an input argument and on entry119:* contains the triangular factor U or L from the Cholesky120:* factorization A = U**H*U or A = L*L**H of the band matrix121:* A, in the same storage format as A (see AB). If EQUED = 'Y',122:* then AFB is the factored form of the equilibrated matrix A.123:*124:* If FACT = 'N', then AFB is an output argument and on exit125:* returns the triangular factor U or L from the Cholesky126:* factorization A = U**H*U or A = L*L**H.127:*128:* If FACT = 'E', then AFB is an output argument and on exit129:* returns the triangular factor U or L from the Cholesky130:* factorization A = U**H*U or A = L*L**H of the equilibrated131:* matrix A (see the description of A for the form of the132:* equilibrated matrix).133:*134:* LDAFB (input) INTEGER135:* The leading dimension of the array AFB. LDAFB >= KD+1.136:*137:* EQUED (input or output) CHARACTER*1138:* Specifies the form of equilibration that was done.139:* = 'N': No equilibration (always true if FACT = 'N').140:* = 'Y': Equilibration was done, i.e., A has been replaced by141:* diag(S) * A * diag(S).142:* EQUED is an input argument if FACT = 'F'; otherwise, it is an143:* output argument.144:*145:* S (input or output) REAL array, dimension (N)146:* The scale factors for A; not accessed if EQUED = 'N'. S is147:* an input argument if FACT = 'F'; otherwise, S is an output148:* argument. If FACT = 'F' and EQUED = 'Y', each element of S149:* must be positive.150:*151:* B (input/output) COMPLEX array, dimension (LDB,NRHS)152:* On entry, the N-by-NRHS right hand side matrix B.153:* On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',154:* B is overwritten by diag(S) * B.155:*156:* LDB (input) INTEGER157:* The leading dimension of the array B. LDB >= max(1,N).158:*159:* X (output) COMPLEX array, dimension (LDX,NRHS)160:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to161:* the original system of equations. Note that if EQUED = 'Y',162:* A and B are modified on exit, and the solution to the163:* equilibrated system is inv(diag(S))*X.164:*165:* LDX (input) INTEGER166:* The leading dimension of the array X. LDX >= max(1,N).167:*168:* RCOND (output) REAL169:* The estimate of the reciprocal condition number of the matrix170:* A after equilibration (if done). If RCOND is less than the171:* machine precision (in particular, if RCOND = 0), the matrix172:* is singular to working precision. This condition is173:* indicated by a return code of INFO > 0.174:*175:* FERR (output) REAL array, dimension (NRHS)176:* The estimated forward error bound for each solution vector177:* X(j) (the j-th column of the solution matrix X).178:* If XTRUE is the true solution corresponding to X(j), FERR(j)179:* is an estimated upper bound for the magnitude of the largest180:* element in (X(j) - XTRUE) divided by the magnitude of the181:* largest element in X(j). The estimate is as reliable as182:* the estimate for RCOND, and is almost always a slight183:* overestimate of the true error.184:*185:* BERR (output) REAL array, dimension (NRHS)186:* The componentwise relative backward error of each solution187:* vector X(j) (i.e., the smallest relative change in188:* any element of A or B that makes X(j) an exact solution).189:*190:* WORK (workspace) COMPLEX array, dimension (2*N)191:*192:* RWORK (workspace) REAL array, dimension (N)193:*194:* INFO (output) INTEGER195:* = 0: successful exit196:* < 0: if INFO = -i, the i-th argument had an illegal value197:* > 0: if INFO = i, and i is198:* <= N: the leading minor of order i of A is199:* not positive definite, so the factorization200:* could not be completed, and the solution has not201:* been computed. RCOND = 0 is returned.202:* = N+1: U is nonsingular, but RCOND is less than machine203:* precision, meaning that the matrix is singular204:* to working precision. Nevertheless, the205:* solution and error bounds are computed because206:* there are a number of situations where the207:* computed solution can be more accurate than the208:* value of RCOND would suggest.209:*210:* Further Details211:* ===============212:*213:* The band storage scheme is illustrated by the following example, when214:* N = 6, KD = 2, and UPLO = 'U':215:*216:* Two-dimensional storage of the Hermitian matrix A:217:*218:* a11 a12 a13219:* a22 a23 a24220:* a33 a34 a35221:* a44 a45 a46222:* a55 a56223:* (aij=conjg(aji)) a66224:*225:* Band storage of the upper triangle of A:226:*227:* * * a13 a24 a35 a46228:* * a12 a23 a34 a45 a56229:* a11 a22 a33 a44 a55 a66230:*231:* Similarly, if UPLO = 'L' the format of A is as follows:232:*233:* a11 a22 a33 a44 a55 a66234:* a21 a32 a43 a54 a65 *235:* a31 a42 a53 a64 * *236:*237:* Array elements marked * are not used by the routine.238:*239:* =====================================================================240:*241:* .. Parameters ..242: REAL ZERO, ONE 243:PARAMETER( ZERO = 0.0E+0, ONE = 1.0E+0 ) 244:* ..245:* .. Local Scalars ..246:LOGICALEQUIL, NOFACT, RCEQU, UPPER 247: INTEGER I, INFEQU, J, J1, J2 248: REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 249:* ..250:* .. External Functions ..251:LOGICALLSAME 252: REAL CLANHB, SLAMCH 253:EXTERNALLSAME, CLANHB, SLAMCH 254:* ..255:* .. External Subroutines ..256:EXTERNALCCOPY, CLACPY, CLAQHB, CPBCON, CPBEQU, CPBRFS, 257: $ CPBTRF, CPBTRS, XERBLA 258:* ..259:* .. Intrinsic Functions ..260:INTRINSICMAX, MIN 261:* ..262:* .. Executable Statements ..263:*264: INFO = 0 265: NOFACT =LSAME( FACT, 'N' ) 266: EQUIL =LSAME( FACT, 'E' ) 267: UPPER =LSAME( UPLO, 'U' ) 268:IF( NOFACT .OR. EQUIL )THEN269: EQUED = 'N' 270: RCEQU = .FALSE. 271:ELSE272: RCEQU =LSAME( EQUED, 'Y' ) 273: SMLNUM =SLAMCH( 'Safe minimum' ) 274: BIGNUM = ONE / SMLNUM 275:ENDIF276:*277:* Test the input parameters.278:*279:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 280: $THEN281: INFO = -1 282:ELSEIF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) )THEN283: INFO = -2 284:ELSEIF( N.LT.0 )THEN285: INFO = -3 286:ELSEIF( KD.LT.0 )THEN287: INFO = -4 288:ELSEIF( NRHS.LT.0 )THEN289: INFO = -5 290:ELSEIF( LDAB.LT.KD+1 )THEN291: INFO = -7 292:ELSEIF( LDAFB.LT.KD+1 )THEN293: INFO = -9 294:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 295: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN296: INFO = -10 297:ELSE298:IF( RCEQU )THEN299: SMIN = BIGNUM 300: SMAX = ZERO 301:DO10 J = 1, N 302: SMIN =MIN( SMIN,S( J ) ) 303: SMAX =MAX( SMAX,S( J ) ) 304: 10CONTINUE305:IF( SMIN.LE.ZERO )THEN306: INFO = -11 307:ELSEIF( N.GT.0 )THEN308: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 309:ELSE310: SCOND = ONE 311:ENDIF312:ENDIF313:IF( INFO.EQ.0 )THEN314:IF( LDB.LT.MAX( 1, N ) )THEN315: INFO = -13 316:ELSEIF( LDX.LT.MAX( 1, N ) )THEN317: INFO = -15 318:ENDIF319:ENDIF320:ENDIF321:*322:IF( INFO.NE.0 )THEN323:CALLXERBLA( 'CPBSVX', -INFO ) 324:RETURN325:ENDIF326:*327:IF( EQUIL )THEN328:*329:* Compute row and column scalings to equilibrate the matrix A.330:*331:CALLCPBEQU( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, INFEQU ) 332:IF( INFEQU.EQ.0 )THEN333:*334:* Equilibrate the matrix.335:*336:CALLCLAQHB( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED ) 337: RCEQU =LSAME( EQUED, 'Y' ) 338:ENDIF339:ENDIF340:*341:* Scale the right-hand side.342:*343:IF( RCEQU )THEN344:DO30 J = 1, NRHS 345:DO20 I = 1, N 346:B( I, J ) =S( I )*B( I, J ) 347: 20CONTINUE348: 30CONTINUE349:ENDIF350:*351:IF( NOFACT .OR. EQUIL )THEN352:*353:* Compute the Cholesky factorization A = U'*U or A = L*L'.354:*355:IF( UPPER )THEN356:DO40 J = 1, N 357: J1 =MAX( J-KD, 1 ) 358:CALLCCOPY( J-J1+1,AB( KD+1-J+J1, J ), 1, 359: $AFB( KD+1-J+J1, J ), 1 ) 360: 40CONTINUE361:ELSE362:DO50 J = 1, N 363: J2 =MIN( J+KD, N ) 364:CALLCCOPY( J2-J+1,AB( 1, J ), 1,AFB( 1, J ), 1 ) 365: 50CONTINUE366:ENDIF367:*368:CALLCPBTRF( UPLO, N, KD, AFB, LDAFB, INFO ) 369:*370:* Return if INFO is non-zero.371:*372:IF( INFO.GT.0 )THEN373: RCOND = ZERO 374:RETURN375:ENDIF376:ENDIF377:*378:* Compute the norm of the matrix A.379:*380: ANORM =CLANHB( '1', UPLO, N, KD, AB, LDAB, RWORK ) 381:*382:* Compute the reciprocal of the condition number of A.383:*384:CALLCPBCON( UPLO, N, KD, AFB, LDAFB, ANORM, RCOND, WORK, RWORK, 385: $ INFO ) 386:*387:* Compute the solution matrix X.388:*389:CALLCLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 390:CALLCPBTRS( UPLO, N, KD, NRHS, AFB, LDAFB, X, LDX, INFO ) 391:*392:* Use iterative refinement to improve the computed solution and393:* compute error bounds and backward error estimates for it.394:*395:CALLCPBRFS( UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, B, LDB, X, 396: $ LDX, FERR, BERR, WORK, RWORK, INFO ) 397:*398:* Transform the solution matrix X to a solution of the original399:* system.400:*401:IF( RCEQU )THEN402:DO70 J = 1, NRHS 403:DO60 I = 1, N 404:X( I, J ) =S( I )*X( I, J ) 405: 60CONTINUE406: 70CONTINUE407:DO80 J = 1, NRHS 408:FERR( J ) =FERR( J ) / SCOND 409: 80CONTINUE410:ENDIF411:*412:* Set INFO = N+1 if the matrix is singular to working precision.413:*414:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 415: $ INFO = N + 1 416:*417:RETURN418:*419:* End of CPBSVX420:*421:END422: