1: *> \brief \b DLAGTF
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) ).le. 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: *> .lt. 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: *> \date November 2011
153: *
154: *> \ingroup auxOTHERcomputational
155: *
156: * =====================================================================
157: SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
158: *
159: * -- LAPACK computational routine (version 3.4.0) --
160: * -- LAPACK is a software package provided by Univ. of Tennessee, --
161: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162: * November 2011
163: *
164: * .. Scalar Arguments ..
165: INTEGER INFO, N
166: DOUBLE PRECISION LAMBDA, TOL
167: * ..
168: * .. Array Arguments ..
169: INTEGER IN( * )
170: DOUBLE PRECISION A( * ), B( * ), C( * ), D( * )
171: * ..
172: *
173: * =====================================================================
174: *
175: * .. Parameters ..
176: DOUBLE PRECISION ZERO
177: PARAMETER ( ZERO = 0.0D+0 )
178: * ..
179: * .. Local Scalars ..
180: INTEGER K
181: DOUBLE PRECISION EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
182: * ..
183: * .. Intrinsic Functions ..
184: INTRINSIC ABS, MAX
185: * ..
186: * .. External Functions ..
187: DOUBLE PRECISION DLAMCH
188: EXTERNAL DLAMCH
189: * ..
190: * .. External Subroutines ..
191: EXTERNAL XERBLA
192: * ..
193: * .. Executable Statements ..
194: *
195: INFO = 0
196: IF( N.LT.0 ) THEN
197: INFO = -1
198: CALL XERBLA( 'DLAGTF', -INFO )
199: RETURN
200: END IF
201: *
202: IF( N.EQ.0 )
203: $ RETURN
204: *
205: A( 1 ) = A( 1 ) - LAMBDA
206: IN( N ) = 0
207: IF( N.EQ.1 ) THEN
208: IF( A( 1 ).EQ.ZERO )
209: $ IN( 1 ) = 1
210: RETURN
211: END IF
212: *
213: EPS = DLAMCH( 'Epsilon' )
214: *
215: TL = MAX( TOL, EPS )
216: SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) )
217: DO 10 K = 1, N - 1
218: A( K+1 ) = A( K+1 ) - LAMBDA
219: SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) )
220: IF( K.LT.( N-1 ) )
221: $ SCALE2 = SCALE2 + ABS( B( K+1 ) )
222: IF( A( K ).EQ.ZERO ) THEN
223: PIV1 = ZERO
224: ELSE
225: PIV1 = ABS( A( K ) ) / SCALE1
226: END IF
227: IF( C( K ).EQ.ZERO ) THEN
228: IN( K ) = 0
229: PIV2 = ZERO
230: SCALE1 = SCALE2
231: IF( K.LT.( N-1 ) )
232: $ D( K ) = ZERO
233: ELSE
234: PIV2 = ABS( C( K ) ) / SCALE2
235: IF( PIV2.LE.PIV1 ) THEN
236: IN( K ) = 0
237: SCALE1 = SCALE2
238: C( K ) = C( K ) / A( K )
239: A( K+1 ) = A( K+1 ) - C( K )*B( K )
240: IF( K.LT.( N-1 ) )
241: $ D( K ) = ZERO
242: ELSE
243: IN( K ) = 1
244: MULT = A( K ) / C( K )
245: A( K ) = C( K )
246: TEMP = A( K+1 )
247: A( K+1 ) = B( K ) - MULT*TEMP
248: IF( K.LT.( N-1 ) ) THEN
249: D( K ) = B( K+1 )
250: B( K+1 ) = -MULT*D( K )
251: END IF
252: B( K ) = TEMP
253: C( K ) = MULT
254: END IF
255: END IF
256: IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) )
257: $ IN( N ) = K
258: 10 CONTINUE
259: IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) )
260: $ IN( N ) = N
261: *
262: RETURN
263: *
264: * End of DLAGTF
265: *
266: END
CVSweb interface <joel.bertrand@systella.fr>