![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlagts.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dlagts.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
2: *
3: * -- LAPACK auxiliary 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, JOB, N
10: DOUBLE PRECISION TOL
11: * ..
12: * .. Array Arguments ..
13: INTEGER IN( * )
14: DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DLAGTS may be used to solve one of the systems of equations
21: *
22: * (T - lambda*I)*x = y or (T - lambda*I)'*x = y,
23: *
24: * where T is an n by n tridiagonal matrix, for x, following the
25: * factorization of (T - lambda*I) as
26: *
27: * (T - lambda*I) = P*L*U ,
28: *
29: * by routine DLAGTF. The choice of equation to be solved is
30: * controlled by the argument JOB, and in each case there is an option
31: * to perturb zero or very small diagonal elements of U, this option
32: * being intended for use in applications such as inverse iteration.
33: *
34: * Arguments
35: * =========
36: *
37: * JOB (input) INTEGER
38: * Specifies the job to be performed by DLAGTS as follows:
39: * = 1: The equations (T - lambda*I)x = y are to be solved,
40: * but diagonal elements of U are not to be perturbed.
41: * = -1: The equations (T - lambda*I)x = y are to be solved
42: * and, if overflow would otherwise occur, the diagonal
43: * elements of U are to be perturbed. See argument TOL
44: * below.
45: * = 2: The equations (T - lambda*I)'x = y are to be solved,
46: * but diagonal elements of U are not to be perturbed.
47: * = -2: The equations (T - lambda*I)'x = y are to be solved
48: * and, if overflow would otherwise occur, the diagonal
49: * elements of U are to be perturbed. See argument TOL
50: * below.
51: *
52: * N (input) INTEGER
53: * The order of the matrix T.
54: *
55: * A (input) DOUBLE PRECISION array, dimension (N)
56: * On entry, A must contain the diagonal elements of U as
57: * returned from DLAGTF.
58: *
59: * B (input) DOUBLE PRECISION array, dimension (N-1)
60: * On entry, B must contain the first super-diagonal elements of
61: * U as returned from DLAGTF.
62: *
63: * C (input) DOUBLE PRECISION array, dimension (N-1)
64: * On entry, C must contain the sub-diagonal elements of L as
65: * returned from DLAGTF.
66: *
67: * D (input) DOUBLE PRECISION array, dimension (N-2)
68: * On entry, D must contain the second super-diagonal elements
69: * of U as returned from DLAGTF.
70: *
71: * IN (input) INTEGER array, dimension (N)
72: * On entry, IN must contain details of the matrix P as returned
73: * from DLAGTF.
74: *
75: * Y (input/output) DOUBLE PRECISION array, dimension (N)
76: * On entry, the right hand side vector y.
77: * On exit, Y is overwritten by the solution vector x.
78: *
79: * TOL (input/output) DOUBLE PRECISION
80: * On entry, with JOB .lt. 0, TOL should be the minimum
81: * perturbation to be made to very small diagonal elements of U.
82: * TOL should normally be chosen as about eps*norm(U), where eps
83: * is the relative machine precision, but if TOL is supplied as
84: * non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
85: * If JOB .gt. 0 then TOL is not referenced.
86: *
87: * On exit, TOL is changed as described above, only if TOL is
88: * non-positive on entry. Otherwise TOL is unchanged.
89: *
90: * INFO (output) INTEGER
91: * = 0 : successful exit
92: * .lt. 0: if INFO = -i, the i-th argument had an illegal value
93: * .gt. 0: overflow would occur when computing the INFO(th)
94: * element of the solution vector x. This can only occur
95: * when JOB is supplied as positive and either means
96: * that a diagonal element of U is very small, or that
97: * the elements of the right-hand side vector y are very
98: * large.
99: *
100: * =====================================================================
101: *
102: * .. Parameters ..
103: DOUBLE PRECISION ONE, ZERO
104: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
105: * ..
106: * .. Local Scalars ..
107: INTEGER K
108: DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
109: * ..
110: * .. Intrinsic Functions ..
111: INTRINSIC ABS, MAX, SIGN
112: * ..
113: * .. External Functions ..
114: DOUBLE PRECISION DLAMCH
115: EXTERNAL DLAMCH
116: * ..
117: * .. External Subroutines ..
118: EXTERNAL XERBLA
119: * ..
120: * .. Executable Statements ..
121: *
122: INFO = 0
123: IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
124: INFO = -1
125: ELSE IF( N.LT.0 ) THEN
126: INFO = -2
127: END IF
128: IF( INFO.NE.0 ) THEN
129: CALL XERBLA( 'DLAGTS', -INFO )
130: RETURN
131: END IF
132: *
133: IF( N.EQ.0 )
134: $ RETURN
135: *
136: EPS = DLAMCH( 'Epsilon' )
137: SFMIN = DLAMCH( 'Safe minimum' )
138: BIGNUM = ONE / SFMIN
139: *
140: IF( JOB.LT.0 ) THEN
141: IF( TOL.LE.ZERO ) THEN
142: TOL = ABS( A( 1 ) )
143: IF( N.GT.1 )
144: $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
145: DO 10 K = 3, N
146: TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
147: $ ABS( D( K-2 ) ) )
148: 10 CONTINUE
149: TOL = TOL*EPS
150: IF( TOL.EQ.ZERO )
151: $ TOL = EPS
152: END IF
153: END IF
154: *
155: IF( ABS( JOB ).EQ.1 ) THEN
156: DO 20 K = 2, N
157: IF( IN( K-1 ).EQ.0 ) THEN
158: Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
159: ELSE
160: TEMP = Y( K-1 )
161: Y( K-1 ) = Y( K )
162: Y( K ) = TEMP - C( K-1 )*Y( K )
163: END IF
164: 20 CONTINUE
165: IF( JOB.EQ.1 ) THEN
166: DO 30 K = N, 1, -1
167: IF( K.LE.N-2 ) THEN
168: TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
169: ELSE IF( K.EQ.N-1 ) THEN
170: TEMP = Y( K ) - B( K )*Y( K+1 )
171: ELSE
172: TEMP = Y( K )
173: END IF
174: AK = A( K )
175: ABSAK = ABS( AK )
176: IF( ABSAK.LT.ONE ) THEN
177: IF( ABSAK.LT.SFMIN ) THEN
178: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
179: $ THEN
180: INFO = K
181: RETURN
182: ELSE
183: TEMP = TEMP*BIGNUM
184: AK = AK*BIGNUM
185: END IF
186: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
187: INFO = K
188: RETURN
189: END IF
190: END IF
191: Y( K ) = TEMP / AK
192: 30 CONTINUE
193: ELSE
194: DO 50 K = N, 1, -1
195: IF( K.LE.N-2 ) THEN
196: TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
197: ELSE IF( K.EQ.N-1 ) THEN
198: TEMP = Y( K ) - B( K )*Y( K+1 )
199: ELSE
200: TEMP = Y( K )
201: END IF
202: AK = A( K )
203: PERT = SIGN( TOL, AK )
204: 40 CONTINUE
205: ABSAK = ABS( AK )
206: IF( ABSAK.LT.ONE ) THEN
207: IF( ABSAK.LT.SFMIN ) THEN
208: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
209: $ THEN
210: AK = AK + PERT
211: PERT = 2*PERT
212: GO TO 40
213: ELSE
214: TEMP = TEMP*BIGNUM
215: AK = AK*BIGNUM
216: END IF
217: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
218: AK = AK + PERT
219: PERT = 2*PERT
220: GO TO 40
221: END IF
222: END IF
223: Y( K ) = TEMP / AK
224: 50 CONTINUE
225: END IF
226: ELSE
227: *
228: * Come to here if JOB = 2 or -2
229: *
230: IF( JOB.EQ.2 ) THEN
231: DO 60 K = 1, N
232: IF( K.GE.3 ) THEN
233: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
234: ELSE IF( K.EQ.2 ) THEN
235: TEMP = Y( K ) - B( K-1 )*Y( K-1 )
236: ELSE
237: TEMP = Y( K )
238: END IF
239: AK = A( K )
240: ABSAK = ABS( AK )
241: IF( ABSAK.LT.ONE ) THEN
242: IF( ABSAK.LT.SFMIN ) THEN
243: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
244: $ THEN
245: INFO = K
246: RETURN
247: ELSE
248: TEMP = TEMP*BIGNUM
249: AK = AK*BIGNUM
250: END IF
251: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
252: INFO = K
253: RETURN
254: END IF
255: END IF
256: Y( K ) = TEMP / AK
257: 60 CONTINUE
258: ELSE
259: DO 80 K = 1, N
260: IF( K.GE.3 ) THEN
261: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
262: ELSE IF( K.EQ.2 ) THEN
263: TEMP = Y( K ) - B( K-1 )*Y( K-1 )
264: ELSE
265: TEMP = Y( K )
266: END IF
267: AK = A( K )
268: PERT = SIGN( TOL, AK )
269: 70 CONTINUE
270: ABSAK = ABS( AK )
271: IF( ABSAK.LT.ONE ) THEN
272: IF( ABSAK.LT.SFMIN ) THEN
273: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
274: $ THEN
275: AK = AK + PERT
276: PERT = 2*PERT
277: GO TO 70
278: ELSE
279: TEMP = TEMP*BIGNUM
280: AK = AK*BIGNUM
281: END IF
282: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
283: AK = AK + PERT
284: PERT = 2*PERT
285: GO TO 70
286: END IF
287: END IF
288: Y( K ) = TEMP / AK
289: 80 CONTINUE
290: END IF
291: *
292: DO 90 K = N, 2, -1
293: IF( IN( K-1 ).EQ.0 ) THEN
294: Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
295: ELSE
296: TEMP = Y( K-1 )
297: Y( K-1 ) = Y( K )
298: Y( K ) = TEMP - C( K-1 )*Y( K )
299: END IF
300: 90 CONTINUE
301: END IF
302: *
303: * End of DLAGTS
304: *
305: END