LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlagtf.f
Go to the documentation of this file.
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) ).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 September 2012
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.2) --
160 * -- LAPACK is a software package provided by Univ. of Tennessee, --
161 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162 * September 2012
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