LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
cptsvx.f
Go to the documentation of this file.
1 *> \brief <b> CPTSVX computes the solution to system of linear equations A * X = B for PT 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 CPTSVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cptsvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cptsvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cptsvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
22 * RCOND, FERR, BERR, WORK, RWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER FACT
26 * INTEGER INFO, LDB, LDX, N, NRHS
27 * REAL RCOND
28 * ..
29 * .. Array Arguments ..
30 * REAL BERR( * ), D( * ), DF( * ), FERR( * ),
31 * $ RWORK( * )
32 * COMPLEX B( LDB, * ), E( * ), EF( * ), WORK( * ),
33 * $ X( LDX, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> CPTSVX uses the factorization A = L*D*L**H to compute the solution
43 *> to a complex system of linear equations A*X = B, where A is an
44 *> N-by-N Hermitian positive definite tridiagonal matrix and X and B
45 *> are N-by-NRHS matrices.
46 *>
47 *> Error bounds on the solution and a condition estimate are also
48 *> provided.
49 *> \endverbatim
50 *
51 *> \par Description:
52 * =================
53 *>
54 *> \verbatim
55 *>
56 *> The following steps are performed:
57 *>
58 *> 1. If FACT = 'N', the matrix A is factored as A = L*D*L**H, where L
59 *> is a unit lower bidiagonal matrix and D is diagonal. The
60 *> factorization can also be regarded as having the form
61 *> A = U**H*D*U.
62 *>
63 *> 2. If the leading i-by-i principal minor is not positive definite,
64 *> then the routine returns with INFO = i. Otherwise, the factored
65 *> form of A is used to estimate the condition number of the matrix
66 *> A. If the reciprocal of the condition number is less than machine
67 *> precision, INFO = N+1 is returned as a warning, but the routine
68 *> still goes on to solve for X and compute error bounds as
69 *> 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 the matrix
86 *> A is supplied on entry.
87 *> = 'F': On entry, DF and EF contain the factored form of A.
88 *> D, E, DF, and EF will not be modified.
89 *> = 'N': The matrix A will be copied to DF and EF and
90 *> factored.
91 *> \endverbatim
92 *>
93 *> \param[in] N
94 *> \verbatim
95 *> N is INTEGER
96 *> The order of the matrix A. N >= 0.
97 *> \endverbatim
98 *>
99 *> \param[in] NRHS
100 *> \verbatim
101 *> NRHS is INTEGER
102 *> The number of right hand sides, i.e., the number of columns
103 *> of the matrices B and X. NRHS >= 0.
104 *> \endverbatim
105 *>
106 *> \param[in] D
107 *> \verbatim
108 *> D is REAL array, dimension (N)
109 *> The n diagonal elements of the tridiagonal matrix A.
110 *> \endverbatim
111 *>
112 *> \param[in] E
113 *> \verbatim
114 *> E is COMPLEX array, dimension (N-1)
115 *> The (n-1) subdiagonal elements of the tridiagonal matrix A.
116 *> \endverbatim
117 *>
118 *> \param[in,out] DF
119 *> \verbatim
120 *> DF is REAL array, dimension (N)
121 *> If FACT = 'F', then DF is an input argument and on entry
122 *> contains the n diagonal elements of the diagonal matrix D
123 *> from the L*D*L**H factorization of A.
124 *> If FACT = 'N', then DF is an output argument and on exit
125 *> contains the n diagonal elements of the diagonal matrix D
126 *> from the L*D*L**H factorization of A.
127 *> \endverbatim
128 *>
129 *> \param[in,out] EF
130 *> \verbatim
131 *> EF is COMPLEX array, dimension (N-1)
132 *> If FACT = 'F', then EF is an input argument and on entry
133 *> contains the (n-1) subdiagonal elements of the unit
134 *> bidiagonal factor L from the L*D*L**H factorization of A.
135 *> If FACT = 'N', then EF is an output argument and on exit
136 *> contains the (n-1) subdiagonal elements of the unit
137 *> bidiagonal factor L from the L*D*L**H factorization of A.
138 *> \endverbatim
139 *>
140 *> \param[in] B
141 *> \verbatim
142 *> B is COMPLEX array, dimension (LDB,NRHS)
143 *> The N-by-NRHS right hand side matrix B.
144 *> \endverbatim
145 *>
146 *> \param[in] LDB
147 *> \verbatim
148 *> LDB is INTEGER
149 *> The leading dimension of the array B. LDB >= max(1,N).
150 *> \endverbatim
151 *>
152 *> \param[out] X
153 *> \verbatim
154 *> X is COMPLEX array, dimension (LDX,NRHS)
155 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
156 *> \endverbatim
157 *>
158 *> \param[in] LDX
159 *> \verbatim
160 *> LDX is INTEGER
161 *> The leading dimension of the array X. LDX >= max(1,N).
162 *> \endverbatim
163 *>
164 *> \param[out] RCOND
165 *> \verbatim
166 *> RCOND is REAL
167 *> The reciprocal condition number of the matrix A. If RCOND
168 *> is less than the machine precision (in particular, if
169 *> RCOND = 0), the matrix is singular to working precision.
170 *> This condition is indicated by a return code of INFO > 0.
171 *> \endverbatim
172 *>
173 *> \param[out] FERR
174 *> \verbatim
175 *> FERR is REAL array, dimension (NRHS)
176 *> The forward error bound for each solution vector
177 *> X(j) (the j-th column of the solution matrix X).
178 *> If XTRUE is the true solution corresponding to X(j), FERR(j)
179 *> is an estimated upper bound for the magnitude of the largest
180 *> element in (X(j) - XTRUE) divided by the magnitude of the
181 *> largest element in X(j).
182 *> \endverbatim
183 *>
184 *> \param[out] BERR
185 *> \verbatim
186 *> BERR is REAL array, dimension (NRHS)
187 *> The componentwise relative backward error of each solution
188 *> vector X(j) (i.e., the smallest relative change in any
189 *> element of A or B that makes X(j) an exact solution).
190 *> \endverbatim
191 *>
192 *> \param[out] WORK
193 *> \verbatim
194 *> WORK is COMPLEX array, dimension (N)
195 *> \endverbatim
196 *>
197 *> \param[out] RWORK
198 *> \verbatim
199 *> RWORK is REAL array, dimension (N)
200 *> \endverbatim
201 *>
202 *> \param[out] INFO
203 *> \verbatim
204 *> INFO is INTEGER
205 *> = 0: successful exit
206 *> < 0: if INFO = -i, the i-th argument had an illegal value
207 *> > 0: if INFO = i, and i is
208 *> <= N: the leading minor of order i of A is
209 *> not positive definite, so the factorization
210 *> could not be completed, and the solution has not
211 *> been computed. RCOND = 0 is returned.
212 *> = N+1: U is nonsingular, but RCOND is less than machine
213 *> precision, meaning that the matrix is singular
214 *> to working precision. Nevertheless, the
215 *> solution and error bounds are computed because
216 *> there are a number of situations where the
217 *> computed solution can be more accurate than the
218 *> value of RCOND would suggest.
219 *> \endverbatim
220 *
221 * Authors:
222 * ========
223 *
224 *> \author Univ. of Tennessee
225 *> \author Univ. of California Berkeley
226 *> \author Univ. of Colorado Denver
227 *> \author NAG Ltd.
228 *
229 *> \date September 2012
230 *
231 *> \ingroup complexPTsolve
232 *
233 * =====================================================================
234  SUBROUTINE cptsvx( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
235  $ rcond, ferr, berr, work, rwork, info )
236 *
237 * -- LAPACK driver routine (version 3.4.2) --
238 * -- LAPACK is a software package provided by Univ. of Tennessee, --
239 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
240 * September 2012
241 *
242 * .. Scalar Arguments ..
243  CHARACTER fact
244  INTEGER info, ldb, ldx, n, nrhs
245  REAL rcond
246 * ..
247 * .. Array Arguments ..
248  REAL berr( * ), d( * ), df( * ), ferr( * ),
249  $ rwork( * )
250  COMPLEX b( ldb, * ), e( * ), ef( * ), work( * ),
251  $ x( ldx, * )
252 * ..
253 *
254 * =====================================================================
255 *
256 * .. Parameters ..
257  REAL zero
258  parameter( zero = 0.0e+0 )
259 * ..
260 * .. Local Scalars ..
261  LOGICAL nofact
262  REAL anorm
263 * ..
264 * .. External Functions ..
265  LOGICAL lsame
266  REAL clanht, slamch
267  EXTERNAL lsame, clanht, slamch
268 * ..
269 * .. External Subroutines ..
270  EXTERNAL ccopy, clacpy, cptcon, cptrfs, cpttrf, cpttrs,
271  $ scopy, xerbla
272 * ..
273 * .. Intrinsic Functions ..
274  INTRINSIC max
275 * ..
276 * .. Executable Statements ..
277 *
278 * Test the input parameters.
279 *
280  info = 0
281  nofact = lsame( fact, 'N' )
282  IF( .NOT.nofact .AND. .NOT.lsame( fact, 'F' ) ) THEN
283  info = -1
284  ELSE IF( n.LT.0 ) THEN
285  info = -2
286  ELSE IF( nrhs.LT.0 ) THEN
287  info = -3
288  ELSE IF( ldb.LT.max( 1, n ) ) THEN
289  info = -9
290  ELSE IF( ldx.LT.max( 1, n ) ) THEN
291  info = -11
292  END IF
293  IF( info.NE.0 ) THEN
294  CALL xerbla( 'CPTSVX', -info )
295  RETURN
296  END IF
297 *
298  IF( nofact ) THEN
299 *
300 * Compute the L*D*L**H (or U**H*D*U) factorization of A.
301 *
302  CALL scopy( n, d, 1, df, 1 )
303  IF( n.GT.1 )
304  $ CALL ccopy( n-1, e, 1, ef, 1 )
305  CALL cpttrf( n, df, ef, info )
306 *
307 * Return if INFO is non-zero.
308 *
309  IF( info.GT.0 )THEN
310  rcond = zero
311  RETURN
312  END IF
313  END IF
314 *
315 * Compute the norm of the matrix A.
316 *
317  anorm = clanht( '1', n, d, e )
318 *
319 * Compute the reciprocal of the condition number of A.
320 *
321  CALL cptcon( n, df, ef, anorm, rcond, rwork, info )
322 *
323 * Compute the solution vectors X.
324 *
325  CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
326  CALL cpttrs( 'Lower', n, nrhs, df, ef, x, ldx, info )
327 *
328 * Use iterative refinement to improve the computed solutions and
329 * compute error bounds and backward error estimates for them.
330 *
331  CALL cptrfs( 'Lower', n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr,
332  $ berr, work, rwork, info )
333 *
334 * Set INFO = N+1 if the matrix is singular to working precision.
335 *
336  IF( rcond.LT.slamch( 'Epsilon' ) )
337  $ info = n + 1
338 *
339  RETURN
340 *
341 * End of CPTSVX
342 *
343  END