001:       SUBROUTINE SLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
002:      $                   INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       LOGICAL            LREAL, LTRAN
010:       INTEGER            INFO, LDT, N
011:       REAL               SCALE, W
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               B( * ), T( LDT, * ), WORK( * ), X( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  SLAQTR solves the real quasi-triangular system
021: *
022: *               op(T)*p = scale*c,               if LREAL = .TRUE.
023: *
024: *  or the complex quasi-triangular systems
025: *
026: *             op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.
027: *
028: *  in real arithmetic, where T is upper quasi-triangular.
029: *  If LREAL = .FALSE., then the first diagonal block of T must be
030: *  1 by 1, B is the specially structured matrix
031: *
032: *                 B = [ b(1) b(2) ... b(n) ]
033: *                     [       w            ]
034: *                     [           w        ]
035: *                     [              .     ]
036: *                     [                 w  ]
037: *
038: *  op(A) = A or A', A' denotes the conjugate transpose of
039: *  matrix A.
040: *
041: *  On input, X = [ c ].  On output, X = [ p ].
042: *                [ d ]                  [ q ]
043: *
044: *  This subroutine is designed for the condition number estimation
045: *  in routine STRSNA.
046: *
047: *  Arguments
048: *  =========
049: *
050: *  LTRAN   (input) LOGICAL
051: *          On entry, LTRAN specifies the option of conjugate transpose:
052: *             = .FALSE.,    op(T+i*B) = T+i*B,
053: *             = .TRUE.,     op(T+i*B) = (T+i*B)'.
054: *
055: *  LREAL   (input) LOGICAL
056: *          On entry, LREAL specifies the input matrix structure:
057: *             = .FALSE.,    the input is complex
058: *             = .TRUE.,     the input is real
059: *
060: *  N       (input) INTEGER
061: *          On entry, N specifies the order of T+i*B. N >= 0.
062: *
063: *  T       (input) REAL array, dimension (LDT,N)
064: *          On entry, T contains a matrix in Schur canonical form.
065: *          If LREAL = .FALSE., then the first diagonal block of T must
066: *          be 1 by 1.
067: *
068: *  LDT     (input) INTEGER
069: *          The leading dimension of the matrix T. LDT >= max(1,N).
070: *
071: *  B       (input) REAL array, dimension (N)
072: *          On entry, B contains the elements to form the matrix
073: *          B as described above.
074: *          If LREAL = .TRUE., B is not referenced.
075: *
076: *  W       (input) REAL
077: *          On entry, W is the diagonal element of the matrix B.
078: *          If LREAL = .TRUE., W is not referenced.
079: *
080: *  SCALE   (output) REAL
081: *          On exit, SCALE is the scale factor.
082: *
083: *  X       (input/output) REAL array, dimension (2*N)
084: *          On entry, X contains the right hand side of the system.
085: *          On exit, X is overwritten by the solution.
086: *
087: *  WORK    (workspace) REAL array, dimension (N)
088: *
089: *  INFO    (output) INTEGER
090: *          On exit, INFO is set to
091: *             0: successful exit.
092: *               1: the some diagonal 1 by 1 block has been perturbed by
093: *                  a small number SMIN to keep nonsingularity.
094: *               2: the some diagonal 2 by 2 block has been perturbed by
095: *                  a small number in SLALN2 to keep nonsingularity.
096: *          NOTE: In the interests of speed, this routine does not
097: *                check the inputs for errors.
098: *
099: * =====================================================================
100: *
101: *     .. Parameters ..
102:       REAL               ZERO, ONE
103:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
104: *     ..
105: *     .. Local Scalars ..
106:       LOGICAL            NOTRAN
107:       INTEGER            I, IERR, J, J1, J2, JNEXT, K, N1, N2
108:       REAL               BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
109:      $                   SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
110: *     ..
111: *     .. Local Arrays ..
112:       REAL               D( 2, 2 ), V( 2, 2 )
113: *     ..
114: *     .. External Functions ..
115:       INTEGER            ISAMAX
116:       REAL               SASUM, SDOT, SLAMCH, SLANGE
117:       EXTERNAL           ISAMAX, SASUM, SDOT, SLAMCH, SLANGE
118: *     ..
119: *     .. External Subroutines ..
120:       EXTERNAL           SAXPY, SLADIV, SLALN2, SSCAL
121: *     ..
122: *     .. Intrinsic Functions ..
123:       INTRINSIC          ABS, MAX
124: *     ..
125: *     .. Executable Statements ..
126: *
127: *     Do not test the input parameters for errors
128: *
129:       NOTRAN = .NOT.LTRAN
130:       INFO = 0
131: *
132: *     Quick return if possible
133: *
134:       IF( N.EQ.0 )
135:      $   RETURN
136: *
137: *     Set constants to control overflow
138: *
139:       EPS = SLAMCH( 'P' )
140:       SMLNUM = SLAMCH( 'S' ) / EPS
141:       BIGNUM = ONE / SMLNUM
142: *
143:       XNORM = SLANGE( 'M', N, N, T, LDT, D )
144:       IF( .NOT.LREAL )
145:      $   XNORM = MAX( XNORM, ABS( W ), SLANGE( 'M', N, 1, B, N, D ) )
146:       SMIN = MAX( SMLNUM, EPS*XNORM )
147: *
148: *     Compute 1-norm of each column of strictly upper triangular
149: *     part of T to control overflow in triangular solver.
150: *
151:       WORK( 1 ) = ZERO
152:       DO 10 J = 2, N
153:          WORK( J ) = SASUM( J-1, T( 1, J ), 1 )
154:    10 CONTINUE
155: *
156:       IF( .NOT.LREAL ) THEN
157:          DO 20 I = 2, N
158:             WORK( I ) = WORK( I ) + ABS( B( I ) )
159:    20    CONTINUE
160:       END IF
161: *
162:       N2 = 2*N
163:       N1 = N
164:       IF( .NOT.LREAL )
165:      $   N1 = N2
166:       K = ISAMAX( N1, X, 1 )
167:       XMAX = ABS( X( K ) )
168:       SCALE = ONE
169: *
170:       IF( XMAX.GT.BIGNUM ) THEN
171:          SCALE = BIGNUM / XMAX
172:          CALL SSCAL( N1, SCALE, X, 1 )
173:          XMAX = BIGNUM
174:       END IF
175: *
176:       IF( LREAL ) THEN
177: *
178:          IF( NOTRAN ) THEN
179: *
180: *           Solve T*p = scale*c
181: *
182:             JNEXT = N
183:             DO 30 J = N, 1, -1
184:                IF( J.GT.JNEXT )
185:      $            GO TO 30
186:                J1 = J
187:                J2 = J
188:                JNEXT = J - 1
189:                IF( J.GT.1 ) THEN
190:                   IF( T( J, J-1 ).NE.ZERO ) THEN
191:                      J1 = J - 1
192:                      JNEXT = J - 2
193:                   END IF
194:                END IF
195: *
196:                IF( J1.EQ.J2 ) THEN
197: *
198: *                 Meet 1 by 1 diagonal block
199: *
200: *                 Scale to avoid overflow when computing
201: *                     x(j) = b(j)/T(j,j)
202: *
203:                   XJ = ABS( X( J1 ) )
204:                   TJJ = ABS( T( J1, J1 ) )
205:                   TMP = T( J1, J1 )
206:                   IF( TJJ.LT.SMIN ) THEN
207:                      TMP = SMIN
208:                      TJJ = SMIN
209:                      INFO = 1
210:                   END IF
211: *
212:                   IF( XJ.EQ.ZERO )
213:      $               GO TO 30
214: *
215:                   IF( TJJ.LT.ONE ) THEN
216:                      IF( XJ.GT.BIGNUM*TJJ ) THEN
217:                         REC = ONE / XJ
218:                         CALL SSCAL( N, REC, X, 1 )
219:                         SCALE = SCALE*REC
220:                         XMAX = XMAX*REC
221:                      END IF
222:                   END IF
223:                   X( J1 ) = X( J1 ) / TMP
224:                   XJ = ABS( X( J1 ) )
225: *
226: *                 Scale x if necessary to avoid overflow when adding a
227: *                 multiple of column j1 of T.
228: *
229:                   IF( XJ.GT.ONE ) THEN
230:                      REC = ONE / XJ
231:                      IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
232:                         CALL SSCAL( N, REC, X, 1 )
233:                         SCALE = SCALE*REC
234:                      END IF
235:                   END IF
236:                   IF( J1.GT.1 ) THEN
237:                      CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
238:                      K = ISAMAX( J1-1, X, 1 )
239:                      XMAX = ABS( X( K ) )
240:                   END IF
241: *
242:                ELSE
243: *
244: *                 Meet 2 by 2 diagonal block
245: *
246: *                 Call 2 by 2 linear system solve, to take
247: *                 care of possible overflow by scaling factor.
248: *
249:                   D( 1, 1 ) = X( J1 )
250:                   D( 2, 1 ) = X( J2 )
251:                   CALL SLALN2( .FALSE., 2, 1, SMIN, ONE, T( J1, J1 ),
252:      $                         LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
253:      $                         SCALOC, XNORM, IERR )
254:                   IF( IERR.NE.0 )
255:      $               INFO = 2
256: *
257:                   IF( SCALOC.NE.ONE ) THEN
258:                      CALL SSCAL( N, SCALOC, X, 1 )
259:                      SCALE = SCALE*SCALOC
260:                   END IF
261:                   X( J1 ) = V( 1, 1 )
262:                   X( J2 ) = V( 2, 1 )
263: *
264: *                 Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2))
265: *                 to avoid overflow in updating right-hand side.
266: *
267:                   XJ = MAX( ABS( V( 1, 1 ) ), ABS( V( 2, 1 ) ) )
268:                   IF( XJ.GT.ONE ) THEN
269:                      REC = ONE / XJ
270:                      IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
271:      $                   ( BIGNUM-XMAX )*REC ) THEN
272:                         CALL SSCAL( N, REC, X, 1 )
273:                         SCALE = SCALE*REC
274:                      END IF
275:                   END IF
276: *
277: *                 Update right-hand side
278: *
279:                   IF( J1.GT.1 ) THEN
280:                      CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
281:                      CALL SAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
282:                      K = ISAMAX( J1-1, X, 1 )
283:                      XMAX = ABS( X( K ) )
284:                   END IF
285: *
286:                END IF
287: *
288:    30       CONTINUE
289: *
290:          ELSE
291: *
292: *           Solve T'*p = scale*c
293: *
294:             JNEXT = 1
295:             DO 40 J = 1, N
296:                IF( J.LT.JNEXT )
297:      $            GO TO 40
298:                J1 = J
299:                J2 = J
300:                JNEXT = J + 1
301:                IF( J.LT.N ) THEN
302:                   IF( T( J+1, J ).NE.ZERO ) THEN
303:                      J2 = J + 1
304:                      JNEXT = J + 2
305:                   END IF
306:                END IF
307: *
308:                IF( J1.EQ.J2 ) THEN
309: *
310: *                 1 by 1 diagonal block
311: *
312: *                 Scale if necessary to avoid overflow in forming the
313: *                 right-hand side element by inner product.
314: *
315:                   XJ = ABS( X( J1 ) )
316:                   IF( XMAX.GT.ONE ) THEN
317:                      REC = ONE / XMAX
318:                      IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
319:                         CALL SSCAL( N, REC, X, 1 )
320:                         SCALE = SCALE*REC
321:                         XMAX = XMAX*REC
322:                      END IF
323:                   END IF
324: *
325:                   X( J1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X, 1 )
326: *
327:                   XJ = ABS( X( J1 ) )
328:                   TJJ = ABS( T( J1, J1 ) )
329:                   TMP = T( J1, J1 )
330:                   IF( TJJ.LT.SMIN ) THEN
331:                      TMP = SMIN
332:                      TJJ = SMIN
333:                      INFO = 1
334:                   END IF
335: *
336:                   IF( TJJ.LT.ONE ) THEN
337:                      IF( XJ.GT.BIGNUM*TJJ ) THEN
338:                         REC = ONE / XJ
339:                         CALL SSCAL( N, REC, X, 1 )
340:                         SCALE = SCALE*REC
341:                         XMAX = XMAX*REC
342:                      END IF
343:                   END IF
344:                   X( J1 ) = X( J1 ) / TMP
345:                   XMAX = MAX( XMAX, ABS( X( J1 ) ) )
346: *
347:                ELSE
348: *
349: *                 2 by 2 diagonal block
350: *
351: *                 Scale if necessary to avoid overflow in forming the
352: *                 right-hand side elements by inner product.
353: *
354:                   XJ = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ) )
355:                   IF( XMAX.GT.ONE ) THEN
356:                      REC = ONE / XMAX
357:                      IF( MAX( WORK( J2 ), WORK( J1 ) ).GT.( BIGNUM-XJ )*
358:      $                   REC ) THEN
359:                         CALL SSCAL( N, REC, X, 1 )
360:                         SCALE = SCALE*REC
361:                         XMAX = XMAX*REC
362:                      END IF
363:                   END IF
364: *
365:                   D( 1, 1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X,
366:      $                        1 )
367:                   D( 2, 1 ) = X( J2 ) - SDOT( J1-1, T( 1, J2 ), 1, X,
368:      $                        1 )
369: *
370:                   CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J1, J1 ),
371:      $                         LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
372:      $                         SCALOC, XNORM, IERR )
373:                   IF( IERR.NE.0 )
374:      $               INFO = 2
375: *
376:                   IF( SCALOC.NE.ONE ) THEN
377:                      CALL SSCAL( N, SCALOC, X, 1 )
378:                      SCALE = SCALE*SCALOC
379:                   END IF
380:                   X( J1 ) = V( 1, 1 )
381:                   X( J2 ) = V( 2, 1 )
382:                   XMAX = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ), XMAX )
383: *
384:                END IF
385:    40       CONTINUE
386:          END IF
387: *
388:       ELSE
389: *
390:          SMINW = MAX( EPS*ABS( W ), SMIN )
391:          IF( NOTRAN ) THEN
392: *
393: *           Solve (T + iB)*(p+iq) = c+id
394: *
395:             JNEXT = N
396:             DO 70 J = N, 1, -1
397:                IF( J.GT.JNEXT )
398:      $            GO TO 70
399:                J1 = J
400:                J2 = J
401:                JNEXT = J - 1
402:                IF( J.GT.1 ) THEN
403:                   IF( T( J, J-1 ).NE.ZERO ) THEN
404:                      J1 = J - 1
405:                      JNEXT = J - 2
406:                   END IF
407:                END IF
408: *
409:                IF( J1.EQ.J2 ) THEN
410: *
411: *                 1 by 1 diagonal block
412: *
413: *                 Scale if necessary to avoid overflow in division
414: *
415:                   Z = W
416:                   IF( J1.EQ.1 )
417:      $               Z = B( 1 )
418:                   XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
419:                   TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
420:                   TMP = T( J1, J1 )
421:                   IF( TJJ.LT.SMINW ) THEN
422:                      TMP = SMINW
423:                      TJJ = SMINW
424:                      INFO = 1
425:                   END IF
426: *
427:                   IF( XJ.EQ.ZERO )
428:      $               GO TO 70
429: *
430:                   IF( TJJ.LT.ONE ) THEN
431:                      IF( XJ.GT.BIGNUM*TJJ ) THEN
432:                         REC = ONE / XJ
433:                         CALL SSCAL( N2, REC, X, 1 )
434:                         SCALE = SCALE*REC
435:                         XMAX = XMAX*REC
436:                      END IF
437:                   END IF
438:                   CALL SLADIV( X( J1 ), X( N+J1 ), TMP, Z, SR, SI )
439:                   X( J1 ) = SR
440:                   X( N+J1 ) = SI
441:                   XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
442: *
443: *                 Scale x if necessary to avoid overflow when adding a
444: *                 multiple of column j1 of T.
445: *
446:                   IF( XJ.GT.ONE ) THEN
447:                      REC = ONE / XJ
448:                      IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
449:                         CALL SSCAL( N2, REC, X, 1 )
450:                         SCALE = SCALE*REC
451:                      END IF
452:                   END IF
453: *
454:                   IF( J1.GT.1 ) THEN
455:                      CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
456:                      CALL SAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
457:      $                           X( N+1 ), 1 )
458: *
459:                      X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 )
460:                      X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 )
461: *
462:                      XMAX = ZERO
463:                      DO 50 K = 1, J1 - 1
464:                         XMAX = MAX( XMAX, ABS( X( K ) )+
465:      $                         ABS( X( K+N ) ) )
466:    50                CONTINUE
467:                   END IF
468: *
469:                ELSE
470: *
471: *                 Meet 2 by 2 diagonal block
472: *
473:                   D( 1, 1 ) = X( J1 )
474:                   D( 2, 1 ) = X( J2 )
475:                   D( 1, 2 ) = X( N+J1 )
476:                   D( 2, 2 ) = X( N+J2 )
477:                   CALL SLALN2( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ),
478:      $                         LDT, ONE, ONE, D, 2, ZERO, -W, V, 2,
479:      $                         SCALOC, XNORM, IERR )
480:                   IF( IERR.NE.0 )
481:      $               INFO = 2
482: *
483:                   IF( SCALOC.NE.ONE ) THEN
484:                      CALL SSCAL( 2*N, SCALOC, X, 1 )
485:                      SCALE = SCALOC*SCALE
486:                   END IF
487:                   X( J1 ) = V( 1, 1 )
488:                   X( J2 ) = V( 2, 1 )
489:                   X( N+J1 ) = V( 1, 2 )
490:                   X( N+J2 ) = V( 2, 2 )
491: *
492: *                 Scale X(J1), .... to avoid overflow in
493: *                 updating right hand side.
494: *
495:                   XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ),
496:      $                 ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) )
497:                   IF( XJ.GT.ONE ) THEN
498:                      REC = ONE / XJ
499:                      IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
500:      $                   ( BIGNUM-XMAX )*REC ) THEN
501:                         CALL SSCAL( N2, REC, X, 1 )
502:                         SCALE = SCALE*REC
503:                      END IF
504:                   END IF
505: *
506: *                 Update the right-hand side.
507: *
508:                   IF( J1.GT.1 ) THEN
509:                      CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
510:                      CALL SAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
511: *
512:                      CALL SAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
513:      $                           X( N+1 ), 1 )
514:                      CALL SAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1,
515:      $                           X( N+1 ), 1 )
516: *
517:                      X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) +
518:      $                        B( J2 )*X( N+J2 )
519:                      X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) -
520:      $                          B( J2 )*X( J2 )
521: *
522:                      XMAX = ZERO
523:                      DO 60 K = 1, J1 - 1
524:                         XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ),
525:      $                         XMAX )
526:    60                CONTINUE
527:                   END IF
528: *
529:                END IF
530:    70       CONTINUE
531: *
532:          ELSE
533: *
534: *           Solve (T + iB)'*(p+iq) = c+id
535: *
536:             JNEXT = 1
537:             DO 80 J = 1, N
538:                IF( J.LT.JNEXT )
539:      $            GO TO 80
540:                J1 = J
541:                J2 = J
542:                JNEXT = J + 1
543:                IF( J.LT.N ) THEN
544:                   IF( T( J+1, J ).NE.ZERO ) THEN
545:                      J2 = J + 1
546:                      JNEXT = J + 2
547:                   END IF
548:                END IF
549: *
550:                IF( J1.EQ.J2 ) THEN
551: *
552: *                 1 by 1 diagonal block
553: *
554: *                 Scale if necessary to avoid overflow in forming the
555: *                 right-hand side element by inner product.
556: *
557:                   XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
558:                   IF( XMAX.GT.ONE ) THEN
559:                      REC = ONE / XMAX
560:                      IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
561:                         CALL SSCAL( N2, REC, X, 1 )
562:                         SCALE = SCALE*REC
563:                         XMAX = XMAX*REC
564:                      END IF
565:                   END IF
566: *
567:                   X( J1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X, 1 )
568:                   X( N+J1 ) = X( N+J1 ) - SDOT( J1-1, T( 1, J1 ), 1,
569:      $                        X( N+1 ), 1 )
570:                   IF( J1.GT.1 ) THEN
571:                      X( J1 ) = X( J1 ) - B( J1 )*X( N+1 )
572:                      X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 )
573:                   END IF
574:                   XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
575: *
576:                   Z = W
577:                   IF( J1.EQ.1 )
578:      $               Z = B( 1 )
579: *
580: *                 Scale if necessary to avoid overflow in
581: *                 complex division
582: *
583:                   TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
584:                   TMP = T( J1, J1 )
585:                   IF( TJJ.LT.SMINW ) THEN
586:                      TMP = SMINW
587:                      TJJ = SMINW
588:                      INFO = 1
589:                   END IF
590: *
591:                   IF( TJJ.LT.ONE ) THEN
592:                      IF( XJ.GT.BIGNUM*TJJ ) THEN
593:                         REC = ONE / XJ
594:                         CALL SSCAL( N2, REC, X, 1 )
595:                         SCALE = SCALE*REC
596:                         XMAX = XMAX*REC
597:                      END IF
598:                   END IF
599:                   CALL SLADIV( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI )
600:                   X( J1 ) = SR
601:                   X( J1+N ) = SI
602:                   XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX )
603: *
604:                ELSE
605: *
606: *                 2 by 2 diagonal block
607: *
608: *                 Scale if necessary to avoid overflow in forming the
609: *                 right-hand side element by inner product.
610: *
611:                   XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
612:      $                 ABS( X( J2 ) )+ABS( X( N+J2 ) ) )
613:                   IF( XMAX.GT.ONE ) THEN
614:                      REC = ONE / XMAX
615:                      IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
616:      $                   ( BIGNUM-XJ ) / XMAX ) THEN
617:                         CALL SSCAL( N2, REC, X, 1 )
618:                         SCALE = SCALE*REC
619:                         XMAX = XMAX*REC
620:                      END IF
621:                   END IF
622: *
623:                   D( 1, 1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X,
624:      $                        1 )
625:                   D( 2, 1 ) = X( J2 ) - SDOT( J1-1, T( 1, J2 ), 1, X,
626:      $                        1 )
627:                   D( 1, 2 ) = X( N+J1 ) - SDOT( J1-1, T( 1, J1 ), 1,
628:      $                        X( N+1 ), 1 )
629:                   D( 2, 2 ) = X( N+J2 ) - SDOT( J1-1, T( 1, J2 ), 1,
630:      $                        X( N+1 ), 1 )
631:                   D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 )
632:                   D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 )
633:                   D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 )
634:                   D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 )
635: *
636:                   CALL SLALN2( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ),
637:      $                         LDT, ONE, ONE, D, 2, ZERO, W, V, 2,
638:      $                         SCALOC, XNORM, IERR )
639:                   IF( IERR.NE.0 )
640:      $               INFO = 2
641: *
642:                   IF( SCALOC.NE.ONE ) THEN
643:                      CALL SSCAL( N2, SCALOC, X, 1 )
644:                      SCALE = SCALOC*SCALE
645:                   END IF
646:                   X( J1 ) = V( 1, 1 )
647:                   X( J2 ) = V( 2, 1 )
648:                   X( N+J1 ) = V( 1, 2 )
649:                   X( N+J2 ) = V( 2, 2 )
650:                   XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
651:      $                   ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX )
652: *
653:                END IF
654: *
655:    80       CONTINUE
656: *
657:          END IF
658: *
659:       END IF
660: *
661:       RETURN
662: *
663: *     End of SLAQTR
664: *
665:       END
666: