1 SUBROUTINE SGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, 2 $ DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, 3 $ WORK, IWORK, INFO ) 4* 5* -- LAPACK routine (version 3.0) -- 6* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 7* Courant Institute, Argonne National Lab, and Rice University 8* June 30, 1999 9* 10* .. Scalar Arguments .. 11 CHARACTER FACT, TRANS 12 INTEGER INFO, LDB, LDX, N, NRHS 13 REAL RCOND 14* .. 15* .. Array Arguments .. 16 INTEGER IPIV( * ), IWORK( * ) 17 REAL B( LDB, * ), BERR( * ), D( * ), DF( * ), 18 $ DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ), 19 $ FERR( * ), WORK( * ), X( LDX, * ) 20* .. 21* 22* Purpose 23* ======= 24* 25* SGTSVX uses the LU factorization to compute the solution to a real 26* system of linear equations A * X = B or A**T * X = B, 27* where A is a tridiagonal matrix of order N and X and B are N-by-NRHS 28* matrices. 29* 30* Error bounds on the solution and a condition estimate are also 31* provided. 32* 33* Description 34* =========== 35* 36* The following steps are performed: 37* 38* 1. If FACT = 'N', the LU decomposition is used to factor the matrix A 39* as A = L * U, where L is a product of permutation and unit lower 40* bidiagonal matrices and U is upper triangular with nonzeros in 41* only the main diagonal and first two superdiagonals. 42* 43* 2. If some U(i,i)=0, so that U is exactly singular, then the routine 44* returns with INFO = i. Otherwise, the factored form of A is used 45* to estimate the condition number of the matrix A. If the 46* reciprocal of the condition number is less than machine precision, 47* INFO = N+1 is returned as a warning, but the routine still goes on 48* to solve for X and compute error bounds as described below. 49* 50* 3. The system of equations is solved for X using the factored form 51* of A. 52* 53* 4. Iterative refinement is applied to improve the computed solution 54* matrix and calculate error bounds and backward error estimates 55* for it. 56* 57* Arguments 58* ========= 59* 60* FACT (input) CHARACTER*1 61* Specifies whether or not the factored form of A has been 62* supplied on entry. 63* = 'F': DLF, DF, DUF, DU2, and IPIV contain the factored 64* form of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV 65* will not be modified. 66* = 'N': The matrix will be copied to DLF, DF, and DUF 67* and factored. 68* 69* TRANS (input) CHARACTER*1 70* Specifies the form of the system of equations: 71* = 'N': A * X = B (No transpose) 72* = 'T': A**T * X = B (Transpose) 73* = 'C': A**H * X = B (Conjugate transpose = Transpose) 74* 75* N (input) INTEGER 76* The order of the matrix A. N >= 0. 77* 78* NRHS (input) INTEGER 79* The number of right hand sides, i.e., the number of columns 80* of the matrix B. NRHS >= 0. 81* 82* DL (input) REAL array, dimension (N-1) 83* The (n-1) subdiagonal elements of A. 84* 85* D (input) REAL array, dimension (N) 86* The n diagonal elements of A. 87* 88* DU (input) REAL array, dimension (N-1) 89* The (n-1) superdiagonal elements of A. 90* 91* DLF (input or output) REAL array, dimension (N-1) 92* If FACT = 'F', then DLF is an input argument and on entry 93* contains the (n-1) multipliers that define the matrix L from 94* the LU factorization of A as computed by SGTTRF. 95* 96* If FACT = 'N', then DLF is an output argument and on exit 97* contains the (n-1) multipliers that define the matrix L from 98* the LU factorization of A. 99* 100* DF (input or output) REAL array, dimension (N) 101* If FACT = 'F', then DF is an input argument and on entry 102* contains the n diagonal elements of the upper triangular 103* matrix U from the LU factorization of A. 104* 105* If FACT = 'N', then DF is an output argument and on exit 106* contains the n diagonal elements of the upper triangular 107* matrix U from the LU factorization of A. 108* 109* DUF (input or output) REAL array, dimension (N-1) 110* If FACT = 'F', then DUF is an input argument and on entry 111* contains the (n-1) elements of the first superdiagonal of U. 112* 113* If FACT = 'N', then DUF is an output argument and on exit 114* contains the (n-1) elements of the first superdiagonal of U. 115* 116* DU2 (input or output) REAL array, dimension (N-2) 117* If FACT = 'F', then DU2 is an input argument and on entry 118* contains the (n-2) elements of the second superdiagonal of 119* U. 120* 121* If FACT = 'N', then DU2 is an output argument and on exit 122* contains the (n-2) elements of the second superdiagonal of 123* U. 124* 125* IPIV (input or output) INTEGER array, dimension (N) 126* If FACT = 'F', then IPIV is an input argument and on entry 127* contains the pivot indices from the LU factorization of A as 128* computed by SGTTRF. 129* 130* If FACT = 'N', then IPIV is an output argument and on exit 131* contains the pivot indices from the LU factorization of A; 132* row i of the matrix was interchanged with row IPIV(i). 133* IPIV(i) will always be either i or i+1; IPIV(i) = i indicates 134* a row interchange was not required. 135* 136* B (input) REAL array, dimension (LDB,NRHS) 137* The N-by-NRHS right hand side matrix B. 138* 139* LDB (input) INTEGER 140* The leading dimension of the array B. LDB >= max(1,N). 141* 142* X (output) REAL array, dimension (LDX,NRHS) 143* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. 144* 145* LDX (input) INTEGER 146* The leading dimension of the array X. LDX >= max(1,N). 147* 148* RCOND (output) REAL 149* The estimate of the reciprocal condition number of the matrix 150* A. If RCOND is less than the machine precision (in 151* particular, if RCOND = 0), the matrix is singular to working 152* precision. This condition is indicated by a return code of 153* INFO > 0. 154* 155* FERR (output) REAL array, dimension (NRHS) 156* The estimated forward error bound for each solution vector 157* X(j) (the j-th column of the solution matrix X). 158* If XTRUE is the true solution corresponding to X(j), FERR(j) 159* is an estimated upper bound for the magnitude of the largest 160* element in (X(j) - XTRUE) divided by the magnitude of the 161* largest element in X(j). The estimate is as reliable as 162* the estimate for RCOND, and is almost always a slight 163* overestimate of the true error. 164* 165* BERR (output) REAL array, dimension (NRHS) 166* The componentwise relative backward error of each solution 167* vector X(j) (i.e., the smallest relative change in 168* any element of A or B that makes X(j) an exact solution). 169* 170* WORK (workspace) REAL array, dimension (3*N) 171* 172* IWORK (workspace) INTEGER array, dimension (N) 173* 174* INFO (output) INTEGER 175* = 0: successful exit 176* < 0: if INFO = -i, the i-th argument had an illegal value 177* > 0: if INFO = i, and i is 178* <= N: U(i,i) is exactly zero. The factorization 179* has not been completed unless i = N, but the 180* factor U is exactly singular, so the solution 181* and error bounds could not be computed. 182* RCOND = 0 is returned. 183* = N+1: U is nonsingular, but RCOND is less than machine 184* precision, meaning that the matrix is singular 185* to working precision. Nevertheless, the 186* solution and error bounds are computed because 187* there are a number of situations where the 188* computed solution can be more accurate than the 189* value of RCOND would suggest. 190* 191* ===================================================================== 192* 193* .. Parameters .. 194 REAL ZERO 195 PARAMETER ( ZERO = 0.0E+0 ) 196* .. 197* .. Local Scalars .. 198 LOGICAL NOFACT, NOTRAN 199 CHARACTER NORM 200 REAL ANORM 201* .. 202* .. External Functions .. 203 LOGICAL LSAME 204 REAL SLAMCH, SLANGT 205 EXTERNAL LSAME, SLAMCH, SLANGT 206* .. 207* .. External Subroutines .. 208 EXTERNAL SCOPY, SGTCON, SGTRFS, SGTTRF, SGTTRS, SLACPY, 209 $ XERBLA 210* .. 211* .. Intrinsic Functions .. 212 INTRINSIC MAX 213* .. 214* .. Executable Statements .. 215* 216 INFO = 0 217 NOFACT = LSAME( FACT, 'N' ) 218 NOTRAN = LSAME( TRANS, 'N' ) 219 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 220 INFO = -1 221 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 222 $ LSAME( TRANS, 'C' ) ) THEN 223 INFO = -2 224 ELSE IF( N.LT.0 ) THEN 225 INFO = -3 226 ELSE IF( NRHS.LT.0 ) THEN 227 INFO = -4 228 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 229 INFO = -14 230 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 231 INFO = -16 232 END IF 233 IF( INFO.NE.0 ) THEN 234 CALL XERBLA( 'SGTSVX', -INFO ) 235 RETURN 236 END IF 237* 238 IF( NOFACT ) THEN 239* 240* Compute the LU factorization of A. 241* 242 CALL SCOPY( N, D, 1, DF, 1 ) 243 IF( N.GT.1 ) THEN 244 CALL SCOPY( N-1, DL, 1, DLF, 1 ) 245 CALL SCOPY( N-1, DU, 1, DUF, 1 ) 246 END IF 247 CALL SGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO ) 248* 249* Return if INFO is non-zero. 250* 251 IF( INFO.NE.0 ) THEN 252 IF( INFO.GT.0 ) 253 $ RCOND = ZERO 254 RETURN 255 END IF 256 END IF 257* 258* Compute the norm of the matrix A. 259* 260 IF( NOTRAN ) THEN 261 NORM = '1' 262 ELSE 263 NORM = 'I' 264 END IF 265 ANORM = SLANGT( NORM, N, DL, D, DU ) 266* 267* Compute the reciprocal of the condition number of A. 268* 269 CALL SGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK, 270 $ IWORK, INFO ) 271* 272* Set INFO = N+1 if the matrix is singular to working precision. 273* 274 IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 275 $ INFO = N + 1 276* 277* Compute the solution vectors X. 278* 279 CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 280 CALL SGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX, 281 $ INFO ) 282* 283* Use iterative refinement to improve the computed solutions and 284* compute error bounds and backward error estimates for them. 285* 286 CALL SGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV, 287 $ B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO ) 288* 289 RETURN 290* 291* End of SGTSVX 292* 293 END 294