LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zptsvx.f
Go to the documentation of this file.
1 *> \brief <b> ZPTSVX 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 ZPTSVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zptsvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zptsvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zptsvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZPTSVX( 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 * DOUBLE PRECISION RCOND
28 * ..
29 * .. Array Arguments ..
30 * DOUBLE PRECISION BERR( * ), D( * ), DF( * ), FERR( * ),
31 * $ RWORK( * )
32 * COMPLEX*16 B( LDB, * ), E( * ), EF( * ), WORK( * ),
33 * $ X( LDX, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> ZPTSVX 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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION 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*16 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*16 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*16 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 array, dimension (N)
195 *> \endverbatim
196 *>
197 *> \param[out] RWORK
198 *> \verbatim
199 *> RWORK is DOUBLE PRECISION 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 *> \ingroup complex16PTsolve
230 *
231 * =====================================================================
232  SUBROUTINE zptsvx( 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  DOUBLE PRECISION RCOND
243 * ..
244 * .. Array Arguments ..
245  DOUBLE PRECISION BERR( * ), D( * ), DF( * ), FERR( * ),
246  $ rwork( * )
247  COMPLEX*16 B( LDB, * ), E( * ), EF( * ), WORK( * ),
248  $ x( ldx, * )
249 * ..
250 *
251 * =====================================================================
252 *
253 * .. Parameters ..
254  DOUBLE PRECISION ZERO
255  parameter( zero = 0.0d+0 )
256 * ..
257 * .. Local Scalars ..
258  LOGICAL NOFACT
259  DOUBLE PRECISION ANORM
260 * ..
261 * .. External Functions ..
262  LOGICAL LSAME
263  DOUBLE PRECISION DLAMCH, ZLANHT
264  EXTERNAL lsame, dlamch, zlanht
265 * ..
266 * .. External Subroutines ..
267  EXTERNAL dcopy, xerbla, zcopy, zlacpy, zptcon, zptrfs,
268  $ zpttrf, zpttrs
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( 'ZPTSVX', -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 dcopy( n, d, 1, df, 1 )
300  IF( n.GT.1 )
301  $ CALL zcopy( n-1, e, 1, ef, 1 )
302  CALL zpttrf( 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 = zlanht( '1', n, d, e )
315 *
316 * Compute the reciprocal of the condition number of A.
317 *
318  CALL zptcon( n, df, ef, anorm, rcond, rwork, info )
319 *
320 * Compute the solution vectors X.
321 *
322  CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
323  CALL zpttrs( '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 zptrfs( '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.dlamch( 'Epsilon' ) )
334  $ info = n + 1
335 *
336  RETURN
337 *
338 * End of ZPTSVX
339 *
340  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zpttrs(UPLO, N, NRHS, D, E, B, LDB, INFO)
ZPTTRS
Definition: zpttrs.f:121
subroutine zptrfs(UPLO, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZPTRFS
Definition: zptrfs.f:183
subroutine zpttrf(N, D, E, INFO)
ZPTTRF
Definition: zpttrf.f:92
subroutine zptcon(N, D, E, ANORM, RCOND, RWORK, INFO)
ZPTCON
Definition: zptcon.f:119
subroutine zptsvx(FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO)
ZPTSVX computes the solution to system of linear equations A * X = B for PT matrices
Definition: zptsvx.f:234
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82