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