LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dgtsvx.f
Go to the documentation of this file.
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 *> \date September 2012
288 *
289 *> \ingroup doubleGTsolve
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 driver routine (version 3.4.2) --
297 * -- LAPACK is a software package provided by Univ. of Tennessee, --
298 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
299 * September 2012
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