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