001:       SUBROUTINE SLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            INFO, JOB, N
009:       REAL               TOL
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IN( * )
013:       REAL               A( * ), B( * ), C( * ), D( * ), Y( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  SLAGTS may be used to solve one of the systems of equations
020: *
021: *     (T - lambda*I)*x = y   or   (T - lambda*I)'*x = y,
022: *
023: *  where T is an n by n tridiagonal matrix, for x, following the
024: *  factorization of (T - lambda*I) as
025: *
026: *     (T - lambda*I) = P*L*U ,
027: *
028: *  by routine SLAGTF. The choice of equation to be solved is
029: *  controlled by the argument JOB, and in each case there is an option
030: *  to perturb zero or very small diagonal elements of U, this option
031: *  being intended for use in applications such as inverse iteration.
032: *
033: *  Arguments
034: *  =========
035: *
036: *  JOB     (input) INTEGER
037: *          Specifies the job to be performed by SLAGTS as follows:
038: *          =  1: The equations  (T - lambda*I)x = y  are to be solved,
039: *                but diagonal elements of U are not to be perturbed.
040: *          = -1: The equations  (T - lambda*I)x = y  are to be solved
041: *                and, if overflow would otherwise occur, the diagonal
042: *                elements of U are to be perturbed. See argument TOL
043: *                below.
044: *          =  2: The equations  (T - lambda*I)'x = y  are to be solved,
045: *                but diagonal elements of U are not to be perturbed.
046: *          = -2: The equations  (T - lambda*I)'x = y  are to be solved
047: *                and, if overflow would otherwise occur, the diagonal
048: *                elements of U are to be perturbed. See argument TOL
049: *                below.
050: *
051: *  N       (input) INTEGER
052: *          The order of the matrix T.
053: *
054: *  A       (input) REAL array, dimension (N)
055: *          On entry, A must contain the diagonal elements of U as
056: *          returned from SLAGTF.
057: *
058: *  B       (input) REAL array, dimension (N-1)
059: *          On entry, B must contain the first super-diagonal elements of
060: *          U as returned from SLAGTF.
061: *
062: *  C       (input) REAL array, dimension (N-1)
063: *          On entry, C must contain the sub-diagonal elements of L as
064: *          returned from SLAGTF.
065: *
066: *  D       (input) REAL array, dimension (N-2)
067: *          On entry, D must contain the second super-diagonal elements
068: *          of U as returned from SLAGTF.
069: *
070: *  IN      (input) INTEGER array, dimension (N)
071: *          On entry, IN must contain details of the matrix P as returned
072: *          from SLAGTF.
073: *
074: *  Y       (input/output) REAL array, dimension (N)
075: *          On entry, the right hand side vector y.
076: *          On exit, Y is overwritten by the solution vector x.
077: *
078: *  TOL     (input/output) REAL
079: *          On entry, with  JOB .lt. 0, TOL should be the minimum
080: *          perturbation to be made to very small diagonal elements of U.
081: *          TOL should normally be chosen as about eps*norm(U), where eps
082: *          is the relative machine precision, but if TOL is supplied as
083: *          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
084: *          If  JOB .gt. 0  then TOL is not referenced.
085: *
086: *          On exit, TOL is changed as described above, only if TOL is
087: *          non-positive on entry. Otherwise TOL is unchanged.
088: *
089: *  INFO    (output) INTEGER
090: *          = 0   : successful exit
091: *          .lt. 0: if INFO = -i, the i-th argument had an illegal value
092: *          .gt. 0: overflow would occur when computing the INFO(th)
093: *                  element of the solution vector x. This can only occur
094: *                  when JOB is supplied as positive and either means
095: *                  that a diagonal element of U is very small, or that
096: *                  the elements of the right-hand side vector y are very
097: *                  large.
098: *
099: *  =====================================================================
100: *
101: *     .. Parameters ..
102:       REAL               ONE, ZERO
103:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
104: *     ..
105: *     .. Local Scalars ..
106:       INTEGER            K
107:       REAL               ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
108: *     ..
109: *     .. Intrinsic Functions ..
110:       INTRINSIC          ABS, MAX, SIGN
111: *     ..
112: *     .. External Functions ..
113:       REAL               SLAMCH
114:       EXTERNAL           SLAMCH
115: *     ..
116: *     .. External Subroutines ..
117:       EXTERNAL           XERBLA
118: *     ..
119: *     .. Executable Statements ..
120: *
121:       INFO = 0
122:       IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
123:          INFO = -1
124:       ELSE IF( N.LT.0 ) THEN
125:          INFO = -2
126:       END IF
127:       IF( INFO.NE.0 ) THEN
128:          CALL XERBLA( 'SLAGTS', -INFO )
129:          RETURN
130:       END IF
131: *
132:       IF( N.EQ.0 )
133:      $   RETURN
134: *
135:       EPS = SLAMCH( 'Epsilon' )
136:       SFMIN = SLAMCH( 'Safe minimum' )
137:       BIGNUM = ONE / SFMIN
138: *
139:       IF( JOB.LT.0 ) THEN
140:          IF( TOL.LE.ZERO ) THEN
141:             TOL = ABS( A( 1 ) )
142:             IF( N.GT.1 )
143:      $         TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
144:             DO 10 K = 3, N
145:                TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
146:      $               ABS( D( K-2 ) ) )
147:    10       CONTINUE
148:             TOL = TOL*EPS
149:             IF( TOL.EQ.ZERO )
150:      $         TOL = EPS
151:          END IF
152:       END IF
153: *
154:       IF( ABS( JOB ).EQ.1 ) THEN
155:          DO 20 K = 2, N
156:             IF( IN( K-1 ).EQ.0 ) THEN
157:                Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
158:             ELSE
159:                TEMP = Y( K-1 )
160:                Y( K-1 ) = Y( K )
161:                Y( K ) = TEMP - C( K-1 )*Y( K )
162:             END IF
163:    20    CONTINUE
164:          IF( JOB.EQ.1 ) THEN
165:             DO 30 K = N, 1, -1
166:                IF( K.LE.N-2 ) THEN
167:                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
168:                ELSE IF( K.EQ.N-1 ) THEN
169:                   TEMP = Y( K ) - B( K )*Y( K+1 )
170:                ELSE
171:                   TEMP = Y( K )
172:                END IF
173:                AK = A( K )
174:                ABSAK = ABS( AK )
175:                IF( ABSAK.LT.ONE ) THEN
176:                   IF( ABSAK.LT.SFMIN ) THEN
177:                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
178:      $                    THEN
179:                         INFO = K
180:                         RETURN
181:                      ELSE
182:                         TEMP = TEMP*BIGNUM
183:                         AK = AK*BIGNUM
184:                      END IF
185:                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
186:                      INFO = K
187:                      RETURN
188:                   END IF
189:                END IF
190:                Y( K ) = TEMP / AK
191:    30       CONTINUE
192:          ELSE
193:             DO 50 K = N, 1, -1
194:                IF( K.LE.N-2 ) THEN
195:                   TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
196:                ELSE IF( K.EQ.N-1 ) THEN
197:                   TEMP = Y( K ) - B( K )*Y( K+1 )
198:                ELSE
199:                   TEMP = Y( K )
200:                END IF
201:                AK = A( K )
202:                PERT = SIGN( TOL, AK )
203:    40          CONTINUE
204:                ABSAK = ABS( AK )
205:                IF( ABSAK.LT.ONE ) THEN
206:                   IF( ABSAK.LT.SFMIN ) THEN
207:                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
208:      $                    THEN
209:                         AK = AK + PERT
210:                         PERT = 2*PERT
211:                         GO TO 40
212:                      ELSE
213:                         TEMP = TEMP*BIGNUM
214:                         AK = AK*BIGNUM
215:                      END IF
216:                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
217:                      AK = AK + PERT
218:                      PERT = 2*PERT
219:                      GO TO 40
220:                   END IF
221:                END IF
222:                Y( K ) = TEMP / AK
223:    50       CONTINUE
224:          END IF
225:       ELSE
226: *
227: *        Come to here if  JOB = 2 or -2
228: *
229:          IF( JOB.EQ.2 ) THEN
230:             DO 60 K = 1, N
231:                IF( K.GE.3 ) THEN
232:                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
233:                ELSE IF( K.EQ.2 ) THEN
234:                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
235:                ELSE
236:                   TEMP = Y( K )
237:                END IF
238:                AK = A( K )
239:                ABSAK = ABS( AK )
240:                IF( ABSAK.LT.ONE ) THEN
241:                   IF( ABSAK.LT.SFMIN ) THEN
242:                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
243:      $                    THEN
244:                         INFO = K
245:                         RETURN
246:                      ELSE
247:                         TEMP = TEMP*BIGNUM
248:                         AK = AK*BIGNUM
249:                      END IF
250:                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
251:                      INFO = K
252:                      RETURN
253:                   END IF
254:                END IF
255:                Y( K ) = TEMP / AK
256:    60       CONTINUE
257:          ELSE
258:             DO 80 K = 1, N
259:                IF( K.GE.3 ) THEN
260:                   TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
261:                ELSE IF( K.EQ.2 ) THEN
262:                   TEMP = Y( K ) - B( K-1 )*Y( K-1 )
263:                ELSE
264:                   TEMP = Y( K )
265:                END IF
266:                AK = A( K )
267:                PERT = SIGN( TOL, AK )
268:    70          CONTINUE
269:                ABSAK = ABS( AK )
270:                IF( ABSAK.LT.ONE ) THEN
271:                   IF( ABSAK.LT.SFMIN ) THEN
272:                      IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
273:      $                    THEN
274:                         AK = AK + PERT
275:                         PERT = 2*PERT
276:                         GO TO 70
277:                      ELSE
278:                         TEMP = TEMP*BIGNUM
279:                         AK = AK*BIGNUM
280:                      END IF
281:                   ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
282:                      AK = AK + PERT
283:                      PERT = 2*PERT
284:                      GO TO 70
285:                   END IF
286:                END IF
287:                Y( K ) = TEMP / AK
288:    80       CONTINUE
289:          END IF
290: *
291:          DO 90 K = N, 2, -1
292:             IF( IN( K-1 ).EQ.0 ) THEN
293:                Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
294:             ELSE
295:                TEMP = Y( K-1 )
296:                Y( K-1 ) = Y( K )
297:                Y( K ) = TEMP - C( K-1 )*Y( K )
298:             END IF
299:    90    CONTINUE
300:       END IF
301: *
302: *     End of SLAGTS
303: *
304:       END
305: