LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zgtsvx.f
Go to the documentation of this file.
1 *> \brief <b> ZGTSVX 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 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 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 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 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 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 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 September 2012
289 *
290 *> \ingroup complex16GTsolve
291 *
292 * =====================================================================
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 *
297 * -- LAPACK driver routine (version 3.4.2) --
298 * -- LAPACK is a software package provided by Univ. of Tennessee, --
299 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
300 * September 2012
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
subroutine zgttrs(TRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB, INFO)
ZGTTRS
Definition: zgttrs.f:140
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zgtsvx(FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO)
ZGTSVX computes the solution to system of linear equations A * X = B for GT matrices ...
Definition: zgtsvx.f:296
subroutine zgtcon(NORM, N, DL, D, DU, DU2, IPIV, ANORM, RCOND, WORK, INFO)
ZGTCON
Definition: zgtcon.f:143
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgttrf(N, DL, D, DU, DU2, IPIV, INFO)
ZGTTRF
Definition: zgttrf.f:126
subroutine zgtrfs(TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZGTRFS
Definition: zgtrfs.f:212