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