1*> \brief <b> ZGTSVX computes the solution to system of linear equations A * X = B for GT matrices </b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZGTSVX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgtsvx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgtsvx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgtsvx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, 22* DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, 23* WORK, RWORK, INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER FACT, TRANS 27* INTEGER INFO, LDB, LDX, N, NRHS 28* DOUBLE PRECISION RCOND 29* .. 30* .. Array Arguments .. 31* INTEGER IPIV( * ) 32* DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ) 33* COMPLEX*16 B( LDB, * ), D( * ), DF( * ), DL( * ), 34* $ DLF( * ), DU( * ), DU2( * ), DUF( * ), 35* $ WORK( * ), X( LDX, * ) 36* .. 37* 38* 39*> \par Purpose: 40* ============= 41*> 42*> \verbatim 43*> 44*> ZGTSVX uses the LU factorization to compute the solution to a complex 45*> system of linear equations A * X = B, A**T * X = B, or A**H * X = B, 46*> where A is a tridiagonal matrix of order N and X and B are N-by-NRHS 47*> matrices. 48*> 49*> Error bounds on the solution and a condition estimate are also 50*> provided. 51*> \endverbatim 52* 53*> \par Description: 54* ================= 55*> 56*> \verbatim 57*> 58*> The following steps are performed: 59*> 60*> 1. If FACT = 'N', the LU decomposition is used to factor the matrix A 61*> as A = L * U, where L is a product of permutation and unit lower 62*> bidiagonal matrices and U is upper triangular with nonzeros in 63*> only the main diagonal and first two superdiagonals. 64*> 65*> 2. If some U(i,i)=0, so that U is exactly singular, then the routine 66*> returns with INFO = i. Otherwise, the factored form of A is used 67*> to estimate the condition number of the matrix A. If the 68*> reciprocal of the condition number is less than machine precision, 69*> INFO = N+1 is returned as a warning, but the routine still goes on 70*> to solve for X and compute error bounds as described below. 71*> 72*> 3. The system of equations is solved for X using the factored form 73*> of A. 74*> 75*> 4. Iterative refinement is applied to improve the computed solution 76*> matrix and calculate error bounds and backward error estimates 77*> for it. 78*> \endverbatim 79* 80* Arguments: 81* ========== 82* 83*> \param[in] FACT 84*> \verbatim 85*> FACT is CHARACTER*1 86*> Specifies whether or not the factored form of A has been 87*> supplied on entry. 88*> = 'F': DLF, DF, DUF, DU2, and IPIV contain the factored form 89*> of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV will not 90*> be modified. 91*> = 'N': The matrix will be copied to DLF, DF, and DUF 92*> and factored. 93*> \endverbatim 94*> 95*> \param[in] TRANS 96*> \verbatim 97*> TRANS is CHARACTER*1 98*> Specifies the form of the system of equations: 99*> = 'N': A * X = B (No transpose) 100*> = 'T': A**T * X = B (Transpose) 101*> = 'C': A**H * X = B (Conjugate transpose) 102*> \endverbatim 103*> 104*> \param[in] N 105*> \verbatim 106*> N is INTEGER 107*> The order of the matrix A. N >= 0. 108*> \endverbatim 109*> 110*> \param[in] NRHS 111*> \verbatim 112*> NRHS is INTEGER 113*> The number of right hand sides, i.e., the number of columns 114*> of the matrix B. NRHS >= 0. 115*> \endverbatim 116*> 117*> \param[in] DL 118*> \verbatim 119*> DL is COMPLEX*16 array, dimension (N-1) 120*> The (n-1) subdiagonal elements of A. 121*> \endverbatim 122*> 123*> \param[in] D 124*> \verbatim 125*> D is COMPLEX*16 array, dimension (N) 126*> The n diagonal elements of A. 127*> \endverbatim 128*> 129*> \param[in] DU 130*> \verbatim 131*> DU is COMPLEX*16 array, dimension (N-1) 132*> The (n-1) superdiagonal elements of A. 133*> \endverbatim 134*> 135*> \param[in,out] DLF 136*> \verbatim 137*> DLF is COMPLEX*16 array, dimension (N-1) 138*> If FACT = 'F', then DLF is an input argument and on entry 139*> contains the (n-1) multipliers that define the matrix L from 140*> the LU factorization of A as computed by ZGTTRF. 141*> 142*> If FACT = 'N', then DLF is an output argument and on exit 143*> contains the (n-1) multipliers that define the matrix L from 144*> the LU factorization of A. 145*> \endverbatim 146*> 147*> \param[in,out] DF 148*> \verbatim 149*> DF is COMPLEX*16 array, dimension (N) 150*> If FACT = 'F', then DF is an input argument and on entry 151*> contains the n diagonal elements of the upper triangular 152*> matrix U from the LU factorization of A. 153*> 154*> If FACT = 'N', then DF is an output argument and on exit 155*> contains the n diagonal elements of the upper triangular 156*> matrix U from the LU factorization of A. 157*> \endverbatim 158*> 159*> \param[in,out] DUF 160*> \verbatim 161*> DUF is COMPLEX*16 array, dimension (N-1) 162*> If FACT = 'F', then DUF is an input argument and on entry 163*> contains the (n-1) elements of the first superdiagonal of U. 164*> 165*> If FACT = 'N', then DUF is an output argument and on exit 166*> contains the (n-1) elements of the first superdiagonal of U. 167*> \endverbatim 168*> 169*> \param[in,out] DU2 170*> \verbatim 171*> DU2 is COMPLEX*16 array, dimension (N-2) 172*> If FACT = 'F', then DU2 is an input argument and on entry 173*> contains the (n-2) elements of the second superdiagonal of 174*> U. 175*> 176*> If FACT = 'N', then DU2 is an output argument and on exit 177*> contains the (n-2) elements of the second superdiagonal of 178*> U. 179*> \endverbatim 180*> 181*> \param[in,out] IPIV 182*> \verbatim 183*> IPIV is INTEGER array, dimension (N) 184*> If FACT = 'F', then IPIV is an input argument and on entry 185*> contains the pivot indices from the LU factorization of A as 186*> computed by ZGTTRF. 187*> 188*> If FACT = 'N', then IPIV is an output argument and on exit 189*> contains the pivot indices from the LU factorization of A; 190*> row i of the matrix was interchanged with row IPIV(i). 191*> IPIV(i) will always be either i or i+1; IPIV(i) = i indicates 192*> a row interchange was not required. 193*> \endverbatim 194*> 195*> \param[in] B 196*> \verbatim 197*> B is COMPLEX*16 array, dimension (LDB,NRHS) 198*> The N-by-NRHS right hand side matrix B. 199*> \endverbatim 200*> 201*> \param[in] LDB 202*> \verbatim 203*> LDB is INTEGER 204*> The leading dimension of the array B. LDB >= max(1,N). 205*> \endverbatim 206*> 207*> \param[out] X 208*> \verbatim 209*> X is COMPLEX*16 array, dimension (LDX,NRHS) 210*> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. 211*> \endverbatim 212*> 213*> \param[in] LDX 214*> \verbatim 215*> LDX is INTEGER 216*> The leading dimension of the array X. LDX >= max(1,N). 217*> \endverbatim 218*> 219*> \param[out] RCOND 220*> \verbatim 221*> RCOND is DOUBLE PRECISION 222*> The estimate of the reciprocal condition number of the matrix 223*> A. If RCOND is less than the machine precision (in 224*> particular, if RCOND = 0), the matrix is singular to working 225*> precision. This condition is indicated by a return code of 226*> INFO > 0. 227*> \endverbatim 228*> 229*> \param[out] FERR 230*> \verbatim 231*> FERR is DOUBLE PRECISION array, dimension (NRHS) 232*> The estimated forward error bound for each solution vector 233*> X(j) (the j-th column of the solution matrix X). 234*> If XTRUE is the true solution corresponding to X(j), FERR(j) 235*> is an estimated upper bound for the magnitude of the largest 236*> element in (X(j) - XTRUE) divided by the magnitude of the 237*> largest element in X(j). The estimate is as reliable as 238*> the estimate for RCOND, and is almost always a slight 239*> overestimate of the true error. 240*> \endverbatim 241*> 242*> \param[out] BERR 243*> \verbatim 244*> BERR is DOUBLE PRECISION array, dimension (NRHS) 245*> The componentwise relative backward error of each solution 246*> vector X(j) (i.e., the smallest relative change in 247*> any element of A or B that makes X(j) an exact solution). 248*> \endverbatim 249*> 250*> \param[out] WORK 251*> \verbatim 252*> WORK is COMPLEX*16 array, dimension (2*N) 253*> \endverbatim 254*> 255*> \param[out] RWORK 256*> \verbatim 257*> RWORK is DOUBLE PRECISION array, dimension (N) 258*> \endverbatim 259*> 260*> \param[out] INFO 261*> \verbatim 262*> INFO is INTEGER 263*> = 0: successful exit 264*> < 0: if INFO = -i, the i-th argument had an illegal value 265*> > 0: if INFO = i, and i is 266*> <= N: U(i,i) is exactly zero. The factorization 267*> has not been completed unless i = N, but the 268*> factor U is exactly singular, so the solution 269*> and error bounds could not be computed. 270*> RCOND = 0 is returned. 271*> = N+1: U is nonsingular, but RCOND is less than machine 272*> precision, meaning that the matrix is singular 273*> to working precision. Nevertheless, the 274*> solution and error bounds are computed because 275*> there are a number of situations where the 276*> computed solution can be more accurate than the 277*> value of RCOND would suggest. 278*> \endverbatim 279* 280* Authors: 281* ======== 282* 283*> \author Univ. of Tennessee 284*> \author Univ. of California Berkeley 285*> \author Univ. of Colorado Denver 286*> \author NAG Ltd. 287* 288*> \ingroup complex16GTsolve 289* 290* ===================================================================== 291 SUBROUTINE ZGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, 292 $ DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, 293 $ WORK, RWORK, INFO ) 294* 295* -- LAPACK driver routine -- 296* -- LAPACK is a software package provided by Univ. of Tennessee, -- 297* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 298* 299* .. Scalar Arguments .. 300 CHARACTER FACT, TRANS 301 INTEGER INFO, LDB, LDX, N, NRHS 302 DOUBLE PRECISION RCOND 303* .. 304* .. Array Arguments .. 305 INTEGER IPIV( * ) 306 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ) 307 COMPLEX*16 B( LDB, * ), D( * ), DF( * ), DL( * ), 308 $ DLF( * ), DU( * ), DU2( * ), DUF( * ), 309 $ WORK( * ), X( LDX, * ) 310* .. 311* 312* ===================================================================== 313* 314* .. Parameters .. 315 DOUBLE PRECISION ZERO 316 PARAMETER ( ZERO = 0.0D+0 ) 317* .. 318* .. Local Scalars .. 319 LOGICAL NOFACT, NOTRAN 320 CHARACTER NORM 321 DOUBLE PRECISION ANORM 322* .. 323* .. External Functions .. 324 LOGICAL LSAME 325 DOUBLE PRECISION DLAMCH, ZLANGT 326 EXTERNAL LSAME, DLAMCH, ZLANGT 327* .. 328* .. External Subroutines .. 329 EXTERNAL XERBLA, ZCOPY, ZGTCON, ZGTRFS, ZGTTRF, ZGTTRS, 330 $ ZLACPY 331* .. 332* .. Intrinsic Functions .. 333 INTRINSIC MAX 334* .. 335* .. Executable Statements .. 336* 337 INFO = 0 338 NOFACT = LSAME( FACT, 'N' ) 339 NOTRAN = LSAME( TRANS, 'N' ) 340 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 341 INFO = -1 342 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 343 $ LSAME( TRANS, 'C' ) ) THEN 344 INFO = -2 345 ELSE IF( N.LT.0 ) THEN 346 INFO = -3 347 ELSE IF( NRHS.LT.0 ) THEN 348 INFO = -4 349 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 350 INFO = -14 351 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 352 INFO = -16 353 END IF 354 IF( INFO.NE.0 ) THEN 355 CALL XERBLA( 'ZGTSVX', -INFO ) 356 RETURN 357 END IF 358* 359 IF( NOFACT ) THEN 360* 361* Compute the LU factorization of A. 362* 363 CALL ZCOPY( N, D, 1, DF, 1 ) 364 IF( N.GT.1 ) THEN 365 CALL ZCOPY( N-1, DL, 1, DLF, 1 ) 366 CALL ZCOPY( N-1, DU, 1, DUF, 1 ) 367 END IF 368 CALL ZGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO ) 369* 370* Return if INFO is non-zero. 371* 372 IF( INFO.GT.0 )THEN 373 RCOND = ZERO 374 RETURN 375 END IF 376 END IF 377* 378* Compute the norm of the matrix A. 379* 380 IF( NOTRAN ) THEN 381 NORM = '1' 382 ELSE 383 NORM = 'I' 384 END IF 385 ANORM = ZLANGT( NORM, N, DL, D, DU ) 386* 387* Compute the reciprocal of the condition number of A. 388* 389 CALL ZGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK, 390 $ INFO ) 391* 392* Compute the solution vectors X. 393* 394 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 395 CALL ZGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX, 396 $ INFO ) 397* 398* Use iterative refinement to improve the computed solutions and 399* compute error bounds and backward error estimates for them. 400* 401 CALL ZGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV, 402 $ B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO ) 403* 404* Set INFO = N+1 if the matrix is singular to working precision. 405* 406 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 407 $ INFO = N + 1 408* 409 RETURN 410* 411* End of ZGTSVX 412* 413 END 414