1: *> \brief \b DLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAGTS + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagts.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagts.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagts.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, JOB, N
25: * DOUBLE PRECISION TOL
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IN( * )
29: * DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLAGTS may be used to solve one of the systems of equations
39: *>
40: *> (T - lambda*I)*x = y or (T - lambda*I)**T*x = y,
41: *>
42: *> where T is an n by n tridiagonal matrix, for x, following the
43: *> factorization of (T - lambda*I) as
44: *>
45: *> (T - lambda*I) = P*L*U ,
46: *>
47: *> by routine DLAGTF. The choice of equation to be solved is
48: *> controlled by the argument JOB, and in each case there is an option
49: *> to perturb zero or very small diagonal elements of U, this option
50: *> being intended for use in applications such as inverse iteration.
51: *> \endverbatim
52: *
53: * Arguments:
54: * ==========
55: *
56: *> \param[in] JOB
57: *> \verbatim
58: *> JOB is INTEGER
59: *> Specifies the job to be performed by DLAGTS as follows:
60: *> = 1: The equations (T - lambda*I)x = y are to be solved,
61: *> but diagonal elements of U are not to be perturbed.
62: *> = -1: The equations (T - lambda*I)x = y are to be solved
63: *> and, if overflow would otherwise occur, the diagonal
64: *> elements of U are to be perturbed. See argument TOL
65: *> below.
66: *> = 2: The equations (T - lambda*I)**Tx = y are to be solved,
67: *> but diagonal elements of U are not to be perturbed.
68: *> = -2: The equations (T - lambda*I)**Tx = y are to be solved
69: *> and, if overflow would otherwise occur, the diagonal
70: *> elements of U are to be perturbed. See argument TOL
71: *> below.
72: *> \endverbatim
73: *>
74: *> \param[in] N
75: *> \verbatim
76: *> N is INTEGER
77: *> The order of the matrix T.
78: *> \endverbatim
79: *>
80: *> \param[in] A
81: *> \verbatim
82: *> A is DOUBLE PRECISION array, dimension (N)
83: *> On entry, A must contain the diagonal elements of U as
84: *> returned from DLAGTF.
85: *> \endverbatim
86: *>
87: *> \param[in] B
88: *> \verbatim
89: *> B is DOUBLE PRECISION array, dimension (N-1)
90: *> On entry, B must contain the first super-diagonal elements of
91: *> U as returned from DLAGTF.
92: *> \endverbatim
93: *>
94: *> \param[in] C
95: *> \verbatim
96: *> C is DOUBLE PRECISION array, dimension (N-1)
97: *> On entry, C must contain the sub-diagonal elements of L as
98: *> returned from DLAGTF.
99: *> \endverbatim
100: *>
101: *> \param[in] D
102: *> \verbatim
103: *> D is DOUBLE PRECISION array, dimension (N-2)
104: *> On entry, D must contain the second super-diagonal elements
105: *> of U as returned from DLAGTF.
106: *> \endverbatim
107: *>
108: *> \param[in] IN
109: *> \verbatim
110: *> IN is INTEGER array, dimension (N)
111: *> On entry, IN must contain details of the matrix P as returned
112: *> from DLAGTF.
113: *> \endverbatim
114: *>
115: *> \param[in,out] Y
116: *> \verbatim
117: *> Y is DOUBLE PRECISION array, dimension (N)
118: *> On entry, the right hand side vector y.
119: *> On exit, Y is overwritten by the solution vector x.
120: *> \endverbatim
121: *>
122: *> \param[in,out] TOL
123: *> \verbatim
124: *> TOL is DOUBLE PRECISION
125: *> On entry, with JOB .lt. 0, TOL should be the minimum
126: *> perturbation to be made to very small diagonal elements of U.
127: *> TOL should normally be chosen as about eps*norm(U), where eps
128: *> is the relative machine precision, but if TOL is supplied as
129: *> non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
130: *> If JOB .gt. 0 then TOL is not referenced.
131: *>
132: *> On exit, TOL is changed as described above, only if TOL is
133: *> non-positive on entry. Otherwise TOL is unchanged.
134: *> \endverbatim
135: *>
136: *> \param[out] INFO
137: *> \verbatim
138: *> INFO is INTEGER
139: *> = 0 : successful exit
140: *> .lt. 0: if INFO = -i, the i-th argument had an illegal value
141: *> .gt. 0: overflow would occur when computing the INFO(th)
142: *> element of the solution vector x. This can only occur
143: *> when JOB is supplied as positive and either means
144: *> that a diagonal element of U is very small, or that
145: *> the elements of the right-hand side vector y are very
146: *> large.
147: *> \endverbatim
148: *
149: * Authors:
150: * ========
151: *
152: *> \author Univ. of Tennessee
153: *> \author Univ. of California Berkeley
154: *> \author Univ. of Colorado Denver
155: *> \author NAG Ltd.
156: *
157: *> \date September 2012
158: *
159: *> \ingroup auxOTHERauxiliary
160: *
161: * =====================================================================
162: SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
163: *
164: * -- LAPACK auxiliary routine (version 3.4.2) --
165: * -- LAPACK is a software package provided by Univ. of Tennessee, --
166: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167: * September 2012
168: *
169: * .. Scalar Arguments ..
170: INTEGER INFO, JOB, N
171: DOUBLE PRECISION TOL
172: * ..
173: * .. Array Arguments ..
174: INTEGER IN( * )
175: DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * )
176: * ..
177: *
178: * =====================================================================
179: *
180: * .. Parameters ..
181: DOUBLE PRECISION ONE, ZERO
182: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
183: * ..
184: * .. Local Scalars ..
185: INTEGER K
186: DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
187: * ..
188: * .. Intrinsic Functions ..
189: INTRINSIC ABS, MAX, SIGN
190: * ..
191: * .. External Functions ..
192: DOUBLE PRECISION DLAMCH
193: EXTERNAL DLAMCH
194: * ..
195: * .. External Subroutines ..
196: EXTERNAL XERBLA
197: * ..
198: * .. Executable Statements ..
199: *
200: INFO = 0
201: IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
202: INFO = -1
203: ELSE IF( N.LT.0 ) THEN
204: INFO = -2
205: END IF
206: IF( INFO.NE.0 ) THEN
207: CALL XERBLA( 'DLAGTS', -INFO )
208: RETURN
209: END IF
210: *
211: IF( N.EQ.0 )
212: $ RETURN
213: *
214: EPS = DLAMCH( 'Epsilon' )
215: SFMIN = DLAMCH( 'Safe minimum' )
216: BIGNUM = ONE / SFMIN
217: *
218: IF( JOB.LT.0 ) THEN
219: IF( TOL.LE.ZERO ) THEN
220: TOL = ABS( A( 1 ) )
221: IF( N.GT.1 )
222: $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
223: DO 10 K = 3, N
224: TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
225: $ ABS( D( K-2 ) ) )
226: 10 CONTINUE
227: TOL = TOL*EPS
228: IF( TOL.EQ.ZERO )
229: $ TOL = EPS
230: END IF
231: END IF
232: *
233: IF( ABS( JOB ).EQ.1 ) THEN
234: DO 20 K = 2, N
235: IF( IN( K-1 ).EQ.0 ) THEN
236: Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
237: ELSE
238: TEMP = Y( K-1 )
239: Y( K-1 ) = Y( K )
240: Y( K ) = TEMP - C( K-1 )*Y( K )
241: END IF
242: 20 CONTINUE
243: IF( JOB.EQ.1 ) THEN
244: DO 30 K = N, 1, -1
245: IF( K.LE.N-2 ) THEN
246: TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
247: ELSE IF( K.EQ.N-1 ) THEN
248: TEMP = Y( K ) - B( K )*Y( K+1 )
249: ELSE
250: TEMP = Y( K )
251: END IF
252: AK = A( K )
253: ABSAK = ABS( AK )
254: IF( ABSAK.LT.ONE ) THEN
255: IF( ABSAK.LT.SFMIN ) THEN
256: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
257: $ THEN
258: INFO = K
259: RETURN
260: ELSE
261: TEMP = TEMP*BIGNUM
262: AK = AK*BIGNUM
263: END IF
264: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
265: INFO = K
266: RETURN
267: END IF
268: END IF
269: Y( K ) = TEMP / AK
270: 30 CONTINUE
271: ELSE
272: DO 50 K = N, 1, -1
273: IF( K.LE.N-2 ) THEN
274: TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
275: ELSE IF( K.EQ.N-1 ) THEN
276: TEMP = Y( K ) - B( K )*Y( K+1 )
277: ELSE
278: TEMP = Y( K )
279: END IF
280: AK = A( K )
281: PERT = SIGN( TOL, AK )
282: 40 CONTINUE
283: ABSAK = ABS( AK )
284: IF( ABSAK.LT.ONE ) THEN
285: IF( ABSAK.LT.SFMIN ) THEN
286: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
287: $ THEN
288: AK = AK + PERT
289: PERT = 2*PERT
290: GO TO 40
291: ELSE
292: TEMP = TEMP*BIGNUM
293: AK = AK*BIGNUM
294: END IF
295: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
296: AK = AK + PERT
297: PERT = 2*PERT
298: GO TO 40
299: END IF
300: END IF
301: Y( K ) = TEMP / AK
302: 50 CONTINUE
303: END IF
304: ELSE
305: *
306: * Come to here if JOB = 2 or -2
307: *
308: IF( JOB.EQ.2 ) THEN
309: DO 60 K = 1, N
310: IF( K.GE.3 ) THEN
311: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
312: ELSE IF( K.EQ.2 ) THEN
313: TEMP = Y( K ) - B( K-1 )*Y( K-1 )
314: ELSE
315: TEMP = Y( K )
316: END IF
317: AK = A( K )
318: ABSAK = ABS( AK )
319: IF( ABSAK.LT.ONE ) THEN
320: IF( ABSAK.LT.SFMIN ) THEN
321: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
322: $ THEN
323: INFO = K
324: RETURN
325: ELSE
326: TEMP = TEMP*BIGNUM
327: AK = AK*BIGNUM
328: END IF
329: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
330: INFO = K
331: RETURN
332: END IF
333: END IF
334: Y( K ) = TEMP / AK
335: 60 CONTINUE
336: ELSE
337: DO 80 K = 1, N
338: IF( K.GE.3 ) THEN
339: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
340: ELSE IF( K.EQ.2 ) THEN
341: TEMP = Y( K ) - B( K-1 )*Y( K-1 )
342: ELSE
343: TEMP = Y( K )
344: END IF
345: AK = A( K )
346: PERT = SIGN( TOL, AK )
347: 70 CONTINUE
348: ABSAK = ABS( AK )
349: IF( ABSAK.LT.ONE ) THEN
350: IF( ABSAK.LT.SFMIN ) THEN
351: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
352: $ THEN
353: AK = AK + PERT
354: PERT = 2*PERT
355: GO TO 70
356: ELSE
357: TEMP = TEMP*BIGNUM
358: AK = AK*BIGNUM
359: END IF
360: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
361: AK = AK + PERT
362: PERT = 2*PERT
363: GO TO 70
364: END IF
365: END IF
366: Y( K ) = TEMP / AK
367: 80 CONTINUE
368: END IF
369: *
370: DO 90 K = N, 2, -1
371: IF( IN( K-1 ).EQ.0 ) THEN
372: Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
373: ELSE
374: TEMP = Y( K-1 )
375: Y( K-1 ) = Y( K )
376: Y( K ) = TEMP - C( K-1 )*Y( K )
377: END IF
378: 90 CONTINUE
379: END IF
380: *
381: * End of DLAGTS
382: *
383: END
CVSweb interface <joel.bertrand@systella.fr>
![[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