001:       SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, JOB, N
010:       DOUBLE PRECISION   TOL
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IN( * )
014:       DOUBLE PRECISION   A( * ), B( * ), C( * ), D( * ), Y( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DLAGTS may be used to solve one of the systems of equations
021: *
022: *     (T - lambda*I)*x = y   or   (T - lambda*I)'*x = y,
023: *
024: *  where T is an n by n tridiagonal matrix, for x, following the
025: *  factorization of (T - lambda*I) as
026: *
027: *     (T - lambda*I) = P*L*U ,
028: *
029: *  by routine DLAGTF. The choice of equation to be solved is
030: *  controlled by the argument JOB, and in each case there is an option
031: *  to perturb zero or very small diagonal elements of U, this option
032: *  being intended for use in applications such as inverse iteration.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  JOB     (input) INTEGER
038: *          Specifies the job to be performed by DLAGTS as follows:
039: *          =  1: The equations  (T - lambda*I)x = y  are to be solved,
040: *                but diagonal elements of U are not to be perturbed.
041: *          = -1: The equations  (T - lambda*I)x = y  are to be solved
042: *                and, if overflow would otherwise occur, the diagonal
043: *                elements of U are to be perturbed. See argument TOL
044: *                below.
045: *          =  2: The equations  (T - lambda*I)'x = y  are to be solved,
046: *                but diagonal elements of U are not to be perturbed.
047: *          = -2: The equations  (T - lambda*I)'x = y  are to be solved
048: *                and, if overflow would otherwise occur, the diagonal
049: *                elements of U are to be perturbed. See argument TOL
050: *                below.
051: *
052: *  N       (input) INTEGER
053: *          The order of the matrix T.
054: *
055: *  A       (input) DOUBLE PRECISION array, dimension (N)
056: *          On entry, A must contain the diagonal elements of U as
057: *          returned from DLAGTF.
058: *
059: *  B       (input) DOUBLE PRECISION array, dimension (N-1)
060: *          On entry, B must contain the first super-diagonal elements of
061: *          U as returned from DLAGTF.
062: *
063: *  C       (input) DOUBLE PRECISION array, dimension (N-1)
064: *          On entry, C must contain the sub-diagonal elements of L as
065: *          returned from DLAGTF.
066: *
067: *  D       (input) DOUBLE PRECISION array, dimension (N-2)
068: *          On entry, D must contain the second super-diagonal elements
069: *          of U as returned from DLAGTF.
070: *
071: *  IN      (input) INTEGER array, dimension (N)
072: *          On entry, IN must contain details of the matrix P as returned
073: *          from DLAGTF.
074: *
075: *  Y       (input/output) DOUBLE PRECISION array, dimension (N)
076: *          On entry, the right hand side vector y.
077: *          On exit, Y is overwritten by the solution vector x.
078: *
079: *  TOL     (input/output) DOUBLE PRECISION
080: *          On entry, with  JOB .lt. 0, TOL should be the minimum
081: *          perturbation to be made to very small diagonal elements of U.
082: *          TOL should normally be chosen as about eps*norm(U), where eps
083: *          is the relative machine precision, but if TOL is supplied as
084: *          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
085: *          If  JOB .gt. 0  then TOL is not referenced.
086: *
087: *          On exit, TOL is changed as described above, only if TOL is
088: *          non-positive on entry. Otherwise TOL is unchanged.
089: *
090: *  INFO    (output) INTEGER
091: *          = 0   : successful exit
092: *          .lt. 0: if INFO = -i, the i-th argument had an illegal value
093: *          .gt. 0: overflow would occur when computing the INFO(th)
094: *                  element of the solution vector x. This can only occur
095: *                  when JOB is supplied as positive and either means
096: *                  that a diagonal element of U is very small, or that
097: *                  the elements of the right-hand side vector y are very
098: *                  large.
099: *
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
306: