Annotation of rpl/lapack/lapack/zptsvx.f, revision 1.15
1.13 bertrand 1: *> \brief <b> ZPTSVX computes the solution to system of linear equations A * X = B for PT matrices</b>
1.9 bertrand 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
1.11 bertrand 120: *> DF is DOUBLE PRECISION array, dimension (N)
1.9 bertrand 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
1.11 bertrand 131: *> EF is COMPLEX*16 array, dimension (N-1)
1.9 bertrand 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: *
1.13 bertrand 229: *> \date September 2012
1.9 bertrand 230: *
1.13 bertrand 231: *> \ingroup complex16PTsolve
1.9 bertrand 232: *
233: * =====================================================================
1.1 bertrand 234: SUBROUTINE ZPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
235: $ RCOND, FERR, BERR, WORK, RWORK, INFO )
236: *
1.13 bertrand 237: * -- LAPACK driver routine (version 3.4.2) --
1.1 bertrand 238: * -- LAPACK is a software package provided by Univ. of Tennessee, --
239: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.13 bertrand 240: * September 2012
1.1 bertrand 241: *
242: * .. Scalar Arguments ..
243: CHARACTER FACT
244: INTEGER INFO, LDB, LDX, N, NRHS
245: DOUBLE PRECISION RCOND
246: * ..
247: * .. Array Arguments ..
248: DOUBLE PRECISION BERR( * ), D( * ), DF( * ), FERR( * ),
249: $ RWORK( * )
250: COMPLEX*16 B( LDB, * ), E( * ), EF( * ), WORK( * ),
251: $ X( LDX, * )
252: * ..
253: *
254: * =====================================================================
255: *
256: * .. Parameters ..
257: DOUBLE PRECISION ZERO
258: PARAMETER ( ZERO = 0.0D+0 )
259: * ..
260: * .. Local Scalars ..
261: LOGICAL NOFACT
262: DOUBLE PRECISION ANORM
263: * ..
264: * .. External Functions ..
265: LOGICAL LSAME
266: DOUBLE PRECISION DLAMCH, ZLANHT
267: EXTERNAL LSAME, DLAMCH, ZLANHT
268: * ..
269: * .. External Subroutines ..
270: EXTERNAL DCOPY, XERBLA, ZCOPY, ZLACPY, ZPTCON, ZPTRFS,
271: $ ZPTTRF, ZPTTRS
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( 'ZPTSVX', -INFO )
295: RETURN
296: END IF
297: *
298: IF( NOFACT ) THEN
299: *
1.8 bertrand 300: * Compute the L*D*L**H (or U**H*D*U) factorization of A.
1.1 bertrand 301: *
302: CALL DCOPY( N, D, 1, DF, 1 )
303: IF( N.GT.1 )
304: $ CALL ZCOPY( N-1, E, 1, EF, 1 )
305: CALL ZPTTRF( 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 = ZLANHT( '1', N, D, E )
318: *
319: * Compute the reciprocal of the condition number of A.
320: *
321: CALL ZPTCON( N, DF, EF, ANORM, RCOND, RWORK, INFO )
322: *
323: * Compute the solution vectors X.
324: *
325: CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
326: CALL ZPTTRS( '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 ZPTRFS( '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.DLAMCH( 'Epsilon' ) )
337: $ INFO = N + 1
338: *
339: RETURN
340: *
341: * End of ZPTSVX
342: *
343: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zptsvx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack