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