Annotation of rpl/lapack/lapack/dlagtf.f, revision 1.2
1.1 bertrand 1: SUBROUTINE DLAGTF( N, A, LAMBDA, B, C, TOL, D, IN, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER INFO, N
10: DOUBLE PRECISION LAMBDA, TOL
11: * ..
12: * .. Array Arguments ..
13: INTEGER IN( * )
14: DOUBLE PRECISION A( * ), B( * ), C( * ), D( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
21: * tridiagonal matrix and lambda is a scalar, as
22: *
23: * T - lambda*I = PLU,
24: *
25: * where P is a permutation matrix, L is a unit lower tridiagonal matrix
26: * with at most one non-zero sub-diagonal elements per column and U is
27: * an upper triangular matrix with at most two non-zero super-diagonal
28: * elements per column.
29: *
30: * The factorization is obtained by Gaussian elimination with partial
31: * pivoting and implicit row scaling.
32: *
33: * The parameter LAMBDA is included in the routine so that DLAGTF may
34: * be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
35: * inverse iteration.
36: *
37: * Arguments
38: * =========
39: *
40: * N (input) INTEGER
41: * The order of the matrix T.
42: *
43: * A (input/output) DOUBLE PRECISION array, dimension (N)
44: * On entry, A must contain the diagonal elements of T.
45: *
46: * On exit, A is overwritten by the n diagonal elements of the
47: * upper triangular matrix U of the factorization of T.
48: *
49: * LAMBDA (input) DOUBLE PRECISION
50: * On entry, the scalar lambda.
51: *
52: * B (input/output) DOUBLE PRECISION array, dimension (N-1)
53: * On entry, B must contain the (n-1) super-diagonal elements of
54: * T.
55: *
56: * On exit, B is overwritten by the (n-1) super-diagonal
57: * elements of the matrix U of the factorization of T.
58: *
59: * C (input/output) DOUBLE PRECISION array, dimension (N-1)
60: * On entry, C must contain the (n-1) sub-diagonal elements of
61: * T.
62: *
63: * On exit, C is overwritten by the (n-1) sub-diagonal elements
64: * of the matrix L of the factorization of T.
65: *
66: * TOL (input) DOUBLE PRECISION
67: * On entry, a relative tolerance used to indicate whether or
68: * not the matrix (T - lambda*I) is nearly singular. TOL should
69: * normally be chose as approximately the largest relative error
70: * in the elements of T. For example, if the elements of T are
71: * correct to about 4 significant figures, then TOL should be
72: * set to about 5*10**(-4). If TOL is supplied as less than eps,
73: * where eps is the relative machine precision, then the value
74: * eps is used in place of TOL.
75: *
76: * D (output) DOUBLE PRECISION array, dimension (N-2)
77: * On exit, D is overwritten by the (n-2) second super-diagonal
78: * elements of the matrix U of the factorization of T.
79: *
80: * IN (output) INTEGER array, dimension (N)
81: * On exit, IN contains details of the permutation matrix P. If
82: * an interchange occurred at the kth step of the elimination,
83: * then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
84: * returns the smallest positive integer j such that
85: *
86: * abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
87: *
88: * where norm( A(j) ) denotes the sum of the absolute values of
89: * the jth row of the matrix A. If no such j exists then IN(n)
90: * is returned as zero. If IN(n) is returned as positive, then a
91: * diagonal element of U is small, indicating that
92: * (T - lambda*I) is singular or nearly singular,
93: *
94: * INFO (output) INTEGER
95: * = 0 : successful exit
96: * .lt. 0: if INFO = -k, the kth argument had an illegal value
97: *
98: * =====================================================================
99: *
100: * .. Parameters ..
101: DOUBLE PRECISION ZERO
102: PARAMETER ( ZERO = 0.0D+0 )
103: * ..
104: * .. Local Scalars ..
105: INTEGER K
106: DOUBLE PRECISION EPS, MULT, PIV1, PIV2, SCALE1, SCALE2, TEMP, TL
107: * ..
108: * .. Intrinsic Functions ..
109: INTRINSIC ABS, MAX
110: * ..
111: * .. External Functions ..
112: DOUBLE PRECISION DLAMCH
113: EXTERNAL DLAMCH
114: * ..
115: * .. External Subroutines ..
116: EXTERNAL XERBLA
117: * ..
118: * .. Executable Statements ..
119: *
120: INFO = 0
121: IF( N.LT.0 ) THEN
122: INFO = -1
123: CALL XERBLA( 'DLAGTF', -INFO )
124: RETURN
125: END IF
126: *
127: IF( N.EQ.0 )
128: $ RETURN
129: *
130: A( 1 ) = A( 1 ) - LAMBDA
131: IN( N ) = 0
132: IF( N.EQ.1 ) THEN
133: IF( A( 1 ).EQ.ZERO )
134: $ IN( 1 ) = 1
135: RETURN
136: END IF
137: *
138: EPS = DLAMCH( 'Epsilon' )
139: *
140: TL = MAX( TOL, EPS )
141: SCALE1 = ABS( A( 1 ) ) + ABS( B( 1 ) )
142: DO 10 K = 1, N - 1
143: A( K+1 ) = A( K+1 ) - LAMBDA
144: SCALE2 = ABS( C( K ) ) + ABS( A( K+1 ) )
145: IF( K.LT.( N-1 ) )
146: $ SCALE2 = SCALE2 + ABS( B( K+1 ) )
147: IF( A( K ).EQ.ZERO ) THEN
148: PIV1 = ZERO
149: ELSE
150: PIV1 = ABS( A( K ) ) / SCALE1
151: END IF
152: IF( C( K ).EQ.ZERO ) THEN
153: IN( K ) = 0
154: PIV2 = ZERO
155: SCALE1 = SCALE2
156: IF( K.LT.( N-1 ) )
157: $ D( K ) = ZERO
158: ELSE
159: PIV2 = ABS( C( K ) ) / SCALE2
160: IF( PIV2.LE.PIV1 ) THEN
161: IN( K ) = 0
162: SCALE1 = SCALE2
163: C( K ) = C( K ) / A( K )
164: A( K+1 ) = A( K+1 ) - C( K )*B( K )
165: IF( K.LT.( N-1 ) )
166: $ D( K ) = ZERO
167: ELSE
168: IN( K ) = 1
169: MULT = A( K ) / C( K )
170: A( K ) = C( K )
171: TEMP = A( K+1 )
172: A( K+1 ) = B( K ) - MULT*TEMP
173: IF( K.LT.( N-1 ) ) THEN
174: D( K ) = B( K+1 )
175: B( K+1 ) = -MULT*D( K )
176: END IF
177: B( K ) = TEMP
178: C( K ) = MULT
179: END IF
180: END IF
181: IF( ( MAX( PIV1, PIV2 ).LE.TL ) .AND. ( IN( N ).EQ.0 ) )
182: $ IN( N ) = K
183: 10 CONTINUE
184: IF( ( ABS( A( N ) ).LE.SCALE1*TL ) .AND. ( IN( N ).EQ.0 ) )
185: $ IN( N ) = N
186: *
187: RETURN
188: *
189: * End of DLAGTF
190: *
191: 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