LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
cgtsvx.f
Go to the documentation of this file.
1 *> \brief <b> CGTSVX 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 CGTSVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgtsvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgtsvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgtsvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGTSVX( 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 * REAL RCOND
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IPIV( * )
32 * REAL BERR( * ), FERR( * ), RWORK( * )
33 * COMPLEX 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 *> CGTSVX 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 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 array, dimension (N)
126 *> The n diagonal elements of A.
127 *> \endverbatim
128 *>
129 *> \param[in] DU
130 *> \verbatim
131 *> DU is COMPLEX 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 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 CGTTRF.
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 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 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 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 CGTTRF.
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 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 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 REAL
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 REAL 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 REAL 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 array, dimension (2*N)
253 *> \endverbatim
254 *>
255 *> \param[out] RWORK
256 *> \verbatim
257 *> RWORK is REAL 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 complexGTsolve
291 *
292 * =====================================================================
293  SUBROUTINE cgtsvx( 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  REAL rcond
306 * ..
307 * .. Array Arguments ..
308  INTEGER ipiv( * )
309  REAL berr( * ), ferr( * ), rwork( * )
310  COMPLEX b( ldb, * ), d( * ), df( * ), dl( * ),
311  $ dlf( * ), du( * ), du2( * ), duf( * ),
312  $ work( * ), x( ldx, * )
313 * ..
314 *
315 * =====================================================================
316 *
317 * .. Parameters ..
318  REAL zero
319  parameter( zero = 0.0e+0 )
320 * ..
321 * .. Local Scalars ..
322  LOGICAL nofact, notran
323  CHARACTER norm
324  REAL anorm
325 * ..
326 * .. External Functions ..
327  LOGICAL lsame
328  REAL clangt, slamch
329  EXTERNAL lsame, clangt, slamch
330 * ..
331 * .. External Subroutines ..
332  EXTERNAL ccopy, cgtcon, cgtrfs, cgttrf, cgttrs, clacpy,
333  $ xerbla
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( 'CGTSVX', -info )
359  RETURN
360  END IF
361 *
362  IF( nofact ) THEN
363 *
364 * Compute the LU factorization of A.
365 *
366  CALL ccopy( n, d, 1, df, 1 )
367  IF( n.GT.1 ) THEN
368  CALL ccopy( n-1, dl, 1, dlf, 1 )
369  CALL ccopy( n-1, du, 1, duf, 1 )
370  END IF
371  CALL cgttrf( 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 = clangt( norm, n, dl, d, du )
389 *
390 * Compute the reciprocal of the condition number of A.
391 *
392  CALL cgtcon( norm, n, dlf, df, duf, du2, ipiv, anorm, rcond, work,
393  $ info )
394 *
395 * Compute the solution vectors X.
396 *
397  CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
398  CALL cgttrs( 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 cgtrfs( 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.slamch( 'Epsilon' ) )
410  $ info = n + 1
411 *
412  RETURN
413 *
414 * End of CGTSVX
415 *
416  END