LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 principal minor of order i is not positive,
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 principal minor of order i of A
209*> is not positive, so the factorization could not
210*> be completed, and the solution has not been
211*> 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*> \ingroup ptsvx
230*
231* =====================================================================
232 SUBROUTINE cptsvx( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
233 $ RCOND, FERR, BERR, WORK, RWORK, INFO )
234*
235* -- LAPACK driver routine --
236* -- LAPACK is a software package provided by Univ. of Tennessee, --
237* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
238*
239* .. Scalar Arguments ..
240 CHARACTER FACT
241 INTEGER INFO, LDB, LDX, N, NRHS
242 REAL RCOND
243* ..
244* .. Array Arguments ..
245 REAL BERR( * ), D( * ), DF( * ), FERR( * ),
246 $ rwork( * )
247 COMPLEX B( LDB, * ), E( * ), EF( * ), WORK( * ),
248 $ x( ldx, * )
249* ..
250*
251* =====================================================================
252*
253* .. Parameters ..
254 REAL ZERO
255 parameter( zero = 0.0e+0 )
256* ..
257* .. Local Scalars ..
258 LOGICAL NOFACT
259 REAL ANORM
260* ..
261* .. External Functions ..
262 LOGICAL LSAME
263 REAL CLANHT, SLAMCH
264 EXTERNAL lsame, clanht, slamch
265* ..
266* .. External Subroutines ..
267 EXTERNAL ccopy, clacpy, cptcon, cptrfs, cpttrf, cpttrs,
268 $ scopy, xerbla
269* ..
270* .. Intrinsic Functions ..
271 INTRINSIC max
272* ..
273* .. Executable Statements ..
274*
275* Test the input parameters.
276*
277 info = 0
278 nofact = lsame( fact, 'N' )
279 IF( .NOT.nofact .AND. .NOT.lsame( fact, 'F' ) ) THEN
280 info = -1
281 ELSE IF( n.LT.0 ) THEN
282 info = -2
283 ELSE IF( nrhs.LT.0 ) THEN
284 info = -3
285 ELSE IF( ldb.LT.max( 1, n ) ) THEN
286 info = -9
287 ELSE IF( ldx.LT.max( 1, n ) ) THEN
288 info = -11
289 END IF
290 IF( info.NE.0 ) THEN
291 CALL xerbla( 'CPTSVX', -info )
292 RETURN
293 END IF
294*
295 IF( nofact ) THEN
296*
297* Compute the L*D*L**H (or U**H*D*U) factorization of A.
298*
299 CALL scopy( n, d, 1, df, 1 )
300 IF( n.GT.1 )
301 $ CALL ccopy( n-1, e, 1, ef, 1 )
302 CALL cpttrf( n, df, ef, info )
303*
304* Return if INFO is non-zero.
305*
306 IF( info.GT.0 )THEN
307 rcond = zero
308 RETURN
309 END IF
310 END IF
311*
312* Compute the norm of the matrix A.
313*
314 anorm = clanht( '1', n, d, e )
315*
316* Compute the reciprocal of the condition number of A.
317*
318 CALL cptcon( n, df, ef, anorm, rcond, rwork, info )
319*
320* Compute the solution vectors X.
321*
322 CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
323 CALL cpttrs( 'Lower', n, nrhs, df, ef, x, ldx, info )
324*
325* Use iterative refinement to improve the computed solutions and
326* compute error bounds and backward error estimates for them.
327*
328 CALL cptrfs( 'Lower', n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr,
329 $ berr, work, rwork, info )
330*
331* Set INFO = N+1 if the matrix is singular to working precision.
332*
333 IF( rcond.LT.slamch( 'Epsilon' ) )
334 $ info = n + 1
335*
336 RETURN
337*
338* End of CPTSVX
339*
340 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine cptcon(n, d, e, anorm, rcond, rwork, info)
CPTCON
Definition cptcon.f:119
subroutine cptrfs(uplo, n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr, berr, work, rwork, info)
CPTRFS
Definition cptrfs.f:183
subroutine cptsvx(fact, n, nrhs, d, e, df, ef, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info)
CPTSVX computes the solution to system of linear equations A * X = B for PT matrices
Definition cptsvx.f:234
subroutine cpttrf(n, d, e, info)
CPTTRF
Definition cpttrf.f:92
subroutine cpttrs(uplo, n, nrhs, d, e, b, ldb, info)
CPTTRS
Definition cpttrs.f:121