1: *> \brief \b DLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAGTF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagtf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagtf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagtf.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, N
25: * DOUBLE PRECISION LAMBDA, TOL
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IN( * )
29: * DOUBLE PRECISION A( * ), B( * ), C( * ), D( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
39: *> tridiagonal matrix and lambda is a scalar, as
40: *>
41: *> T - lambda*I = PLU,
42: *>
43: *> where P is a permutation matrix, L is a unit lower tridiagonal matrix
44: *> with at most one non-zero sub-diagonal elements per column and U is
45: *> an upper triangular matrix with at most two non-zero super-diagonal
46: *> elements per column.
47: *>
48: *> The factorization is obtained by Gaussian elimination with partial
49: *> pivoting and implicit row scaling.
50: *>
51: *> The parameter LAMBDA is included in the routine so that DLAGTF may
52: *> be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
53: *> inverse iteration.
54: *> \endverbatim
55: *
56: * Arguments:
57: * ==========
58: *
59: *> \param[in] N
60: *> \verbatim
61: *> N is INTEGER
62: *> The order of the matrix T.
63: *> \endverbatim
64: *>
65: *> \param[in,out] A
66: *> \verbatim
67: *> A is DOUBLE PRECISION array, dimension (N)
68: *> On entry, A must contain the diagonal elements of T.
69: *>
70: *> On exit, A is overwritten by the n diagonal elements of the
71: *> upper triangular matrix U of the factorization of T.
72: *> \endverbatim
73: *>
74: *> \param[in] LAMBDA
75: *> \verbatim
76: *> LAMBDA is DOUBLE PRECISION
77: *> On entry, the scalar lambda.
78: *> \endverbatim
79: *>
80: *> \param[in,out] B
81: *> \verbatim
82: *> B is DOUBLE PRECISION array, dimension (N-1)
83: *> On entry, B must contain the (n-1) super-diagonal elements of
84: *> T.
85: *>
86: *> On exit, B is overwritten by the (n-1) super-diagonal
87: *> elements of the matrix U of the factorization of T.
88: *> \endverbatim
89: *>
90: *> \param[in,out] C
91: *> \verbatim
92: *> C is DOUBLE PRECISION array, dimension (N-1)
93: *> On entry, C must contain the (n-1) sub-diagonal elements of
94: *> T.
95: *>
96: *> On exit, C is overwritten by the (n-1) sub-diagonal elements
97: *> of the matrix L of the factorization of T.
98: *> \endverbatim
99: *>
100: *> \param[in] TOL
101: *> \verbatim
102: *> TOL is DOUBLE PRECISION
103: *> On entry, a relative tolerance used to indicate whether or
104: *> not the matrix (T - lambda*I) is nearly singular. TOL should
105: *> normally be chose as approximately the largest relative error
106: *> in the elements of T. For example, if the elements of T are
107: *> correct to about 4 significant figures, then TOL should be
108: *> set to about 5*10**(-4). If TOL is supplied as less than eps,
109: *> where eps is the relative machine precision, then the value
110: *> eps is used in place of TOL.
111: *> \endverbatim
112: *>
113: *> \param[out] D
114: *> \verbatim
115: *> D is DOUBLE PRECISION array, dimension (N-2)
116: *> On exit, D is overwritten by the (n-2) second super-diagonal
117: *> elements of the matrix U of the factorization of T.
118: *> \endverbatim
119: *>
120: *> \param[out] IN
121: *> \verbatim
122: *> IN is INTEGER array, dimension (N)
123: *> On exit, IN contains details of the permutation matrix P. If
124: *> an interchange occurred at the kth step of the elimination,
125: *> then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
126: *> returns the smallest positive integer j such that
127: *>
128: *> abs( u(j,j) ) <= norm( (T - lambda*I)(j) )*TOL,
129: *>
130: *> where norm( A(j) ) denotes the sum of the absolute values of
131: *> the jth row of the matrix A. If no such j exists then IN(n)
132: *> is returned as zero. If IN(n) is returned as positive, then a
133: *> diagonal element of U is small, indicating that
134: *> (T - lambda*I) is singular or nearly singular,
135: *> \endverbatim
136: *>
137: *> \param[out] INFO
138: *> \verbatim
139: *> INFO is INTEGER
140: *> = 0: successful exit
141: *> < 0: if INFO = -k, the kth argument had an illegal value
142: *> \endverbatim
143: *
144: * Authors:
145: * ========
146: *
147: *> \author Univ. of Tennessee
148: *> \author Univ. of California Berkeley
149: *> \author Univ. of Colorado Denver
150: *> \author NAG Ltd.
151: *
152: *> \ingroup auxOTHERcomputational
153: *
154: * =====================================================================
155: SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
156: *
157: * -- LAPACK computational routine --
158: * -- LAPACK is a software package provided by Univ. of Tennessee, --
159: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160: *
161: * .. Scalar Arguments ..
162: INTEGER INFO, N
163: DOUBLE PRECISION LAMBDA, TOL
164: * ..
165: * .. Array Arguments ..
166: INTEGER IN( * )
167: DOUBLE PRECISION A( * ), B( * ), C( * ), D( * )
168: * ..
169: *
170: * =====================================================================
171: *
172: * .. Parameters ..
173: DOUBLE PRECISION ZERO
174: PARAMETER ( ZERO = 0.0D+0 )
175: * ..
176: * .. Local Scalars ..
177: INTEGER K
178: DOUBLE PRECISION EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
179: * ..
180: * .. Intrinsic Functions ..
181: INTRINSIC ABS, MAX
182: * ..
183: * .. External Functions ..
184: DOUBLE PRECISION DLAMCH
185: EXTERNAL DLAMCH
186: * ..
187: * .. External Subroutines ..
188: EXTERNAL XERBLA
189: * ..
190: * .. Executable Statements ..
191: *
192: INFO = 0
193: IF( N.LT.0 ) THEN
194: INFO = -1
195: CALL XERBLA( 'DLAGTF', -INFO )
196: RETURN
197: END IF
198: *
199: IF( N.EQ.0 )
200: $ RETURN
201: *
202: A( 1 ) = A( 1 ) - LAMBDA
203: IN( N ) = 0
204: IF( N.EQ.1 ) THEN
205: IF( A( 1 ).EQ.ZERO )
206: $ IN( 1 ) = 1
207: RETURN
208: END IF
209: *
210: EPS = DLAMCH( 'Epsilon' )
211: *
212: TL = MAX( TOL, EPS )
213: SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) )
214: DO 10 K = 1, N - 1
215: A( K+1 ) = A( K+1 ) - LAMBDA
216: SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) )
217: IF( K.LT.( N-1 ) )
218: $ SCALE2 = SCALE2 + ABS( B( K+1 ) )
219: IF( A( K ).EQ.ZERO ) THEN
220: PIV1 = ZERO
221: ELSE
222: PIV1 = ABS( A( K ) ) / SCALE1
223: END IF
224: IF( C( K ).EQ.ZERO ) THEN
225: IN( K ) = 0
226: PIV2 = ZERO
227: SCALE1 = SCALE2
228: IF( K.LT.( N-1 ) )
229: $ D( K ) = ZERO
230: ELSE
231: PIV2 = ABS( C( K ) ) / SCALE2
232: IF( PIV2.LE.PIV1 ) THEN
233: IN( K ) = 0
234: SCALE1 = SCALE2
235: C( K ) = C( K ) / A( K )
236: A( K+1 ) = A( K+1 ) - C( K )*B( K )
237: IF( K.LT.( N-1 ) )
238: $ D( K ) = ZERO
239: ELSE
240: IN( K ) = 1
241: MULT = A( K ) / C( K )
242: A( K ) = C( K )
243: TEMP = A( K+1 )
244: A( K+1 ) = B( K ) - MULT*TEMP
245: IF( K.LT.( N-1 ) ) THEN
246: D( K ) = B( K+1 )
247: B( K+1 ) = -MULT*D( K )
248: END IF
249: B( K ) = TEMP
250: C( K ) = MULT
251: END IF
252: END IF
253: IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) )
254: $ IN( N ) = K
255: 10 CONTINUE
256: IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) )
257: $ IN( N ) = N
258: *
259: RETURN
260: *
261: * End of DLAGTF
262: *
263: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlagtf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack