Annotation of rpl/lapack/lapack/zgtsvx.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b ZGTSVX
! 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 or output) 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 or output) 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 or output) 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 or output) 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 or output) 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: *> \date November 2011
! 289: *
! 290: *> \ingroup complex16OTHERcomputational
! 291: *
! 292: * =====================================================================
1.1 bertrand 293: SUBROUTINE ZGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
294: $ DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
295: $ WORK, RWORK, INFO )
296: *
1.8 ! bertrand 297: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 298: * -- LAPACK is a software package provided by Univ. of Tennessee, --
299: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 300: * November 2011
1.1 bertrand 301: *
302: * .. Scalar Arguments ..
303: CHARACTER FACT, TRANS
304: INTEGER INFO, LDB, LDX, N, NRHS
305: DOUBLE PRECISION RCOND
306: * ..
307: * .. Array Arguments ..
308: INTEGER IPIV( * )
309: DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
310: COMPLEX*16 B( LDB, * ), D( * ), DF( * ), DL( * ),
311: $ DLF( * ), DU( * ), DU2( * ), DUF( * ),
312: $ WORK( * ), X( LDX, * )
313: * ..
314: *
315: * =====================================================================
316: *
317: * .. Parameters ..
318: DOUBLE PRECISION ZERO
319: PARAMETER ( ZERO = 0.0D+0 )
320: * ..
321: * .. Local Scalars ..
322: LOGICAL NOFACT, NOTRAN
323: CHARACTER NORM
324: DOUBLE PRECISION ANORM
325: * ..
326: * .. External Functions ..
327: LOGICAL LSAME
328: DOUBLE PRECISION DLAMCH, ZLANGT
329: EXTERNAL LSAME, DLAMCH, ZLANGT
330: * ..
331: * .. External Subroutines ..
332: EXTERNAL XERBLA, ZCOPY, ZGTCON, ZGTRFS, ZGTTRF, ZGTTRS,
333: $ ZLACPY
334: * ..
335: * .. Intrinsic Functions ..
336: INTRINSIC MAX
337: * ..
338: * .. Executable Statements ..
339: *
340: INFO = 0
341: NOFACT = LSAME( FACT, 'N' )
342: NOTRAN = LSAME( TRANS, 'N' )
343: IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
344: INFO = -1
345: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
346: $ LSAME( TRANS, 'C' ) ) THEN
347: INFO = -2
348: ELSE IF( N.LT.0 ) THEN
349: INFO = -3
350: ELSE IF( NRHS.LT.0 ) THEN
351: INFO = -4
352: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
353: INFO = -14
354: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
355: INFO = -16
356: END IF
357: IF( INFO.NE.0 ) THEN
358: CALL XERBLA( 'ZGTSVX', -INFO )
359: RETURN
360: END IF
361: *
362: IF( NOFACT ) THEN
363: *
364: * Compute the LU factorization of A.
365: *
366: CALL ZCOPY( N, D, 1, DF, 1 )
367: IF( N.GT.1 ) THEN
368: CALL ZCOPY( N-1, DL, 1, DLF, 1 )
369: CALL ZCOPY( N-1, DU, 1, DUF, 1 )
370: END IF
371: CALL ZGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO )
372: *
373: * Return if INFO is non-zero.
374: *
375: IF( INFO.GT.0 )THEN
376: RCOND = ZERO
377: RETURN
378: END IF
379: END IF
380: *
381: * Compute the norm of the matrix A.
382: *
383: IF( NOTRAN ) THEN
384: NORM = '1'
385: ELSE
386: NORM = 'I'
387: END IF
388: ANORM = ZLANGT( NORM, N, DL, D, DU )
389: *
390: * Compute the reciprocal of the condition number of A.
391: *
392: CALL ZGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK,
393: $ INFO )
394: *
395: * Compute the solution vectors X.
396: *
397: CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
398: CALL ZGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX,
399: $ INFO )
400: *
401: * Use iterative refinement to improve the computed solutions and
402: * compute error bounds and backward error estimates for them.
403: *
404: CALL ZGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV,
405: $ B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
406: *
407: * Set INFO = N+1 if the matrix is singular to working precision.
408: *
409: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
410: $ INFO = N + 1
411: *
412: RETURN
413: *
414: * End of ZGTSVX
415: *
416: END
CVSweb interface <joel.bertrand@systella.fr>