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: * =====================================================================
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: *
296: * -- LAPACK computational routine (version 3.4.0) --
297: * -- LAPACK is a software package provided by Univ. of Tennessee, --
298: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
299: * November 2011
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>