SRC/slaqtr.f

   1 *> \brief \b SLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of special form, in real arithmetic.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLAQTR + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqtr.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqtr.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqtr.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
  22 *                          INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       LOGICAL            LREAL, LTRAN
  26 *       INTEGER            INFO, LDT, N
  27 *       REAL               SCALE, W
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       REAL               B( * ), T( LDT, * ), WORK( * ), X( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SLAQTR solves the real quasi-triangular system
  40 *>
  41 *>              op(T)*p = scale*c,               if LREAL = .TRUE.
  42 *>
  43 *> or the complex quasi-triangular systems
  44 *>
  45 *>            op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.
  46 *>
  47 *> in real arithmetic, where T is upper quasi-triangular.
  48 *> If LREAL = .FALSE., then the first diagonal block of T must be
  49 *> 1 by 1, B is the specially structured matrix
  50 *>
  51 *>                B = [ b(1) b(2) ... b(n) ]
  52 *>                    [       w            ]
  53 *>                    [           w        ]
  54 *>                    [              .     ]
  55 *>                    [                 w  ]
  56 *>
  57 *> op(A) = A or A**T, A**T denotes the transpose of
  58 *> matrix A.
  59 *>
  60 *> On input, X = [ c ].  On output, X = [ p ].
  61 *>               [ d ]                  [ q ]
  62 *>
  63 *> This subroutine is designed for the condition number estimation
  64 *> in routine STRSNA.
  65 *> \endverbatim
  66 *
  67 *  Arguments:
  68 *  ==========
  69 *
  70 *> \param[in] LTRAN
  71 *> \verbatim
  72 *>          LTRAN is LOGICAL
  73 *>          On entry, LTRAN specifies the option of conjugate transpose:
  74 *>             = .FALSE.,    op(T+i*B) = T+i*B,
  75 *>             = .TRUE.,     op(T+i*B) = (T+i*B)**T.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] LREAL
  79 *> \verbatim
  80 *>          LREAL is LOGICAL
  81 *>          On entry, LREAL specifies the input matrix structure:
  82 *>             = .FALSE.,    the input is complex
  83 *>             = .TRUE.,     the input is real
  84 *> \endverbatim
  85 *>
  86 *> \param[in] N
  87 *> \verbatim
  88 *>          N is INTEGER
  89 *>          On entry, N specifies the order of T+i*B. N >= 0.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] T
  93 *> \verbatim
  94 *>          T is REAL array, dimension (LDT,N)
  95 *>          On entry, T contains a matrix in Schur canonical form.
  96 *>          If LREAL = .FALSE., then the first diagonal block of T must
  97 *>          be 1 by 1.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] LDT
 101 *> \verbatim
 102 *>          LDT is INTEGER
 103 *>          The leading dimension of the matrix T. LDT >= max(1,N).
 104 *> \endverbatim
 105 *>
 106 *> \param[in] B
 107 *> \verbatim
 108 *>          B is REAL array, dimension (N)
 109 *>          On entry, B contains the elements to form the matrix
 110 *>          B as described above.
 111 *>          If LREAL = .TRUE., B is not referenced.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] W
 115 *> \verbatim
 116 *>          W is REAL
 117 *>          On entry, W is the diagonal element of the matrix B.
 118 *>          If LREAL = .TRUE., W is not referenced.
 119 *> \endverbatim
 120 *>
 121 *> \param[out] SCALE
 122 *> \verbatim
 123 *>          SCALE is REAL
 124 *>          On exit, SCALE is the scale factor.
 125 *> \endverbatim
 126 *>
 127 *> \param[in,out] X
 128 *> \verbatim
 129 *>          X is REAL array, dimension (2*N)
 130 *>          On entry, X contains the right hand side of the system.
 131 *>          On exit, X is overwritten by the solution.
 132 *> \endverbatim
 133 *>
 134 *> \param[out] WORK
 135 *> \verbatim
 136 *>          WORK is REAL array, dimension (N)
 137 *> \endverbatim
 138 *>
 139 *> \param[out] INFO
 140 *> \verbatim
 141 *>          INFO is INTEGER
 142 *>          On exit, INFO is set to
 143 *>             0: successful exit.
 144 *>               1: the some diagonal 1 by 1 block has been perturbed by
 145 *>                  a small number SMIN to keep nonsingularity.
 146 *>               2: the some diagonal 2 by 2 block has been perturbed by
 147 *>                  a small number in SLALN2 to keep nonsingularity.
 148 *>          NOTE: In the interests of speed, this routine does not
 149 *>                check the inputs for errors.
 150 *> \endverbatim
 151 *
 152 *  Authors:
 153 *  ========
 154 *
 155 *> \author Univ. of Tennessee
 156 *> \author Univ. of California Berkeley
 157 *> \author Univ. of Colorado Denver
 158 *> \author NAG Ltd.
 159 *
 160 *> \date September 2012
 161 *
 162 *> \ingroup realOTHERauxiliary
 163 *
 164 *  =====================================================================
 165       SUBROUTINE SLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
 166      $                   INFO )
 167 *
 168 *  -- LAPACK auxiliary routine (version 3.4.2) --
 169 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 170 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 171 *     September 2012
 172 *
 173 *     .. Scalar Arguments ..
 174       LOGICAL            LREAL, LTRAN
 175       INTEGER            INFO, LDT, N
 176       REAL               SCALE, W
 177 *     ..
 178 *     .. Array Arguments ..
 179       REAL               B( * ), T( LDT, * ), WORK( * ), X( * )
 180 *     ..
 181 *
 182 * =====================================================================
 183 *
 184 *     .. Parameters ..
 185       REAL               ZERO, ONE
 186       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 187 *     ..
 188 *     .. Local Scalars ..
 189       LOGICAL            NOTRAN
 190       INTEGER            I, IERR, J, J1, J2, JNEXT, K, N1, N2
 191       REAL               BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
 192      $                   SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
 193 *     ..
 194 *     .. Local Arrays ..
 195       REAL               D( 2, 2 ), V( 2, 2 )
 196 *     ..
 197 *     .. External Functions ..
 198       INTEGER            ISAMAX
 199       REAL               SASUM, SDOT, SLAMCH, SLANGE
 200       EXTERNAL           ISAMAX, SASUM, SDOT, SLAMCH, SLANGE
 201 *     ..
 202 *     .. External Subroutines ..
 203       EXTERNAL           SAXPY, SLADIV, SLALN2, SSCAL
 204 *     ..
 205 *     .. Intrinsic Functions ..
 206       INTRINSIC          ABS, MAX
 207 *     ..
 208 *     .. Executable Statements ..
 209 *
 210 *     Do not test the input parameters for errors
 211 *
 212       NOTRAN = .NOT.LTRAN
 213       INFO = 0
 214 *
 215 *     Quick return if possible
 216 *
 217       IF( N.EQ.0 )
 218      $   RETURN
 219 *
 220 *     Set constants to control overflow
 221 *
 222       EPS = SLAMCH( 'P' )
 223       SMLNUM = SLAMCH( 'S' ) / EPS
 224       BIGNUM = ONE / SMLNUM
 225 *
 226       XNORM = SLANGE( 'M', N, N, T, LDT, D )
 227       IF( .NOT.LREAL )
 228      $   XNORM = MAX( XNORM, ABS( W ), SLANGE( 'M', N, 1, B, N, D ) )
 229       SMIN = MAX( SMLNUM, EPS*XNORM )
 230 *
 231 *     Compute 1-norm of each column of strictly upper triangular
 232 *     part of T to control overflow in triangular solver.
 233 *
 234       WORK( 1 ) = ZERO
 235       DO 10 J = 2, N
 236          WORK( J ) = SASUM( J-1, T( 1, J ), 1 )
 237    10 CONTINUE
 238 *
 239       IF( .NOT.LREAL ) THEN
 240          DO 20 I = 2, N
 241             WORK( I ) = WORK( I ) + ABS( B( I ) )
 242    20    CONTINUE
 243       END IF
 244 *
 245       N2 = 2*N
 246       N1 = N
 247       IF( .NOT.LREAL )
 248      $   N1 = N2
 249       K = ISAMAX( N1, X, 1 )
 250       XMAX = ABS( X( K ) )
 251       SCALE = ONE
 252 *
 253       IF( XMAX.GT.BIGNUM ) THEN
 254          SCALE = BIGNUM / XMAX
 255          CALL SSCAL( N1, SCALE, X, 1 )
 256          XMAX = BIGNUM
 257       END IF
 258 *
 259       IF( LREAL ) THEN
 260 *
 261          IF( NOTRAN ) THEN
 262 *
 263 *           Solve T*p = scale*c
 264 *
 265             JNEXT = N
 266             DO 30 J = N, 1, -1
 267                IF( J.GT.JNEXT )
 268      $            GO TO 30
 269                J1 = J
 270                J2 = J
 271                JNEXT = J - 1
 272                IF( J.GT.1 ) THEN
 273                   IF( T( J, J-1 ).NE.ZERO ) THEN
 274                      J1 = J - 1
 275                      JNEXT = J - 2
 276                   END IF
 277                END IF
 278 *
 279                IF( J1.EQ.J2 ) THEN
 280 *
 281 *                 Meet 1 by 1 diagonal block
 282 *
 283 *                 Scale to avoid overflow when computing
 284 *                     x(j) = b(j)/T(j,j)
 285 *
 286                   XJ = ABS( X( J1 ) )
 287                   TJJ = ABS( T( J1, J1 ) )
 288                   TMP = T( J1, J1 )
 289                   IF( TJJ.LT.SMIN ) THEN
 290                      TMP = SMIN
 291                      TJJ = SMIN
 292                      INFO = 1
 293                   END IF
 294 *
 295                   IF( XJ.EQ.ZERO )
 296      $               GO TO 30
 297 *
 298                   IF( TJJ.LT.ONE ) THEN
 299                      IF( XJ.GT.BIGNUM*TJJ ) THEN
 300                         REC = ONE / XJ
 301                         CALL SSCAL( N, REC, X, 1 )
 302                         SCALE = SCALE*REC
 303                         XMAX = XMAX*REC
 304                      END IF
 305                   END IF
 306                   X( J1 ) = X( J1 ) / TMP
 307                   XJ = ABS( X( J1 ) )
 308 *
 309 *                 Scale x if necessary to avoid overflow when adding a
 310 *                 multiple of column j1 of T.
 311 *
 312                   IF( XJ.GT.ONE ) THEN
 313                      REC = ONE / XJ
 314                      IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
 315                         CALL SSCAL( N, REC, X, 1 )
 316                         SCALE = SCALE*REC
 317                      END IF
 318                   END IF
 319                   IF( J1.GT.1 ) THEN
 320                      CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
 321                      K = ISAMAX( J1-1, X, 1 )
 322                      XMAX = ABS( X( K ) )
 323                   END IF
 324 *
 325                ELSE
 326 *
 327 *                 Meet 2 by 2 diagonal block
 328 *
 329 *                 Call 2 by 2 linear system solve, to take
 330 *                 care of possible overflow by scaling factor.
 331 *
 332                   D( 1, 1 ) = X( J1 )
 333                   D( 2, 1 ) = X( J2 )
 334                   CALL SLALN2( .FALSE., 2, 1, SMIN, ONE, T( J1, J1 ),
 335      $                         LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
 336      $                         SCALOC, XNORM, IERR )
 337                   IF( IERR.NE.0 )
 338      $               INFO = 2
 339 *
 340                   IF( SCALOC.NE.ONE ) THEN
 341                      CALL SSCAL( N, SCALOC, X, 1 )
 342                      SCALE = SCALE*SCALOC
 343                   END IF
 344                   X( J1 ) = V( 1, 1 )
 345                   X( J2 ) = V( 2, 1 )
 346 *
 347 *                 Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2))
 348 *                 to avoid overflow in updating right-hand side.
 349 *
 350                   XJ = MAX( ABS( V( 1, 1 ) ), ABS( V( 2, 1 ) ) )
 351                   IF( XJ.GT.ONE ) THEN
 352                      REC = ONE / XJ
 353                      IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
 354      $                   ( BIGNUM-XMAX )*REC ) THEN
 355                         CALL SSCAL( N, REC, X, 1 )
 356                         SCALE = SCALE*REC
 357                      END IF
 358                   END IF
 359 *
 360 *                 Update right-hand side
 361 *
 362                   IF( J1.GT.1 ) THEN
 363                      CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
 364                      CALL SAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
 365                      K = ISAMAX( J1-1, X, 1 )
 366                      XMAX = ABS( X( K ) )
 367                   END IF
 368 *
 369                END IF
 370 *
 371    30       CONTINUE
 372 *
 373          ELSE
 374 *
 375 *           Solve T**T*p = scale*c
 376 *
 377             JNEXT = 1
 378             DO 40 J = 1, N
 379                IF( J.LT.JNEXT )
 380      $            GO TO 40
 381                J1 = J
 382                J2 = J
 383                JNEXT = J + 1
 384                IF( J.LT.N ) THEN
 385                   IF( T( J+1, J ).NE.ZERO ) THEN
 386                      J2 = J + 1
 387                      JNEXT = J + 2
 388                   END IF
 389                END IF
 390 *
 391                IF( J1.EQ.J2 ) THEN
 392 *
 393 *                 1 by 1 diagonal block
 394 *
 395 *                 Scale if necessary to avoid overflow in forming the
 396 *                 right-hand side element by inner product.
 397 *
 398                   XJ = ABS( X( J1 ) )
 399                   IF( XMAX.GT.ONE ) THEN
 400                      REC = ONE / XMAX
 401                      IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
 402                         CALL SSCAL( N, REC, X, 1 )
 403                         SCALE = SCALE*REC
 404                         XMAX = XMAX*REC
 405                      END IF
 406                   END IF
 407 *
 408                   X( J1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X, 1 )
 409 *
 410                   XJ = ABS( X( J1 ) )
 411                   TJJ = ABS( T( J1, J1 ) )
 412                   TMP = T( J1, J1 )
 413                   IF( TJJ.LT.SMIN ) THEN
 414                      TMP = SMIN
 415                      TJJ = SMIN
 416                      INFO = 1
 417                   END IF
 418 *
 419                   IF( TJJ.LT.ONE ) THEN
 420                      IF( XJ.GT.BIGNUM*TJJ ) THEN
 421                         REC = ONE / XJ
 422                         CALL SSCAL( N, REC, X, 1 )
 423                         SCALE = SCALE*REC
 424                         XMAX = XMAX*REC
 425                      END IF
 426                   END IF
 427                   X( J1 ) = X( J1 ) / TMP
 428                   XMAX = MAX( XMAX, ABS( X( J1 ) ) )
 429 *
 430                ELSE
 431 *
 432 *                 2 by 2 diagonal block
 433 *
 434 *                 Scale if necessary to avoid overflow in forming the
 435 *                 right-hand side elements by inner product.
 436 *
 437                   XJ = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ) )
 438                   IF( XMAX.GT.ONE ) THEN
 439                      REC = ONE / XMAX
 440                      IF( MAX( WORK( J2 ), WORK( J1 ) ).GT.( BIGNUM-XJ )*
 441      $                   REC ) THEN
 442                         CALL SSCAL( N, REC, X, 1 )
 443                         SCALE = SCALE*REC
 444                         XMAX = XMAX*REC
 445                      END IF
 446                   END IF
 447 *
 448                   D( 1, 1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X,
 449      $                        1 )
 450                   D( 2, 1 ) = X( J2 ) - SDOT( J1-1, T( 1, J2 ), 1, X,
 451      $                        1 )
 452 *
 453                   CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J1, J1 ),
 454      $                         LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
 455      $                         SCALOC, XNORM, IERR )
 456                   IF( IERR.NE.0 )
 457      $               INFO = 2
 458 *
 459                   IF( SCALOC.NE.ONE ) THEN
 460                      CALL SSCAL( N, SCALOC, X, 1 )
 461                      SCALE = SCALE*SCALOC
 462                   END IF
 463                   X( J1 ) = V( 1, 1 )
 464                   X( J2 ) = V( 2, 1 )
 465                   XMAX = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ), XMAX )
 466 *
 467                END IF
 468    40       CONTINUE
 469          END IF
 470 *
 471       ELSE
 472 *
 473          SMINW = MAX( EPS*ABS( W ), SMIN )
 474          IF( NOTRAN ) THEN
 475 *
 476 *           Solve (T + iB)*(p+iq) = c+id
 477 *
 478             JNEXT = N
 479             DO 70 J = N, 1, -1
 480                IF( J.GT.JNEXT )
 481      $            GO TO 70
 482                J1 = J
 483                J2 = J
 484                JNEXT = J - 1
 485                IF( J.GT.1 ) THEN
 486                   IF( T( J, J-1 ).NE.ZERO ) THEN
 487                      J1 = J - 1
 488                      JNEXT = J - 2
 489                   END IF
 490                END IF
 491 *
 492                IF( J1.EQ.J2 ) THEN
 493 *
 494 *                 1 by 1 diagonal block
 495 *
 496 *                 Scale if necessary to avoid overflow in division
 497 *
 498                   Z = W
 499                   IF( J1.EQ.1 )
 500      $               Z = B( 1 )
 501                   XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
 502                   TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
 503                   TMP = T( J1, J1 )
 504                   IF( TJJ.LT.SMINW ) THEN
 505                      TMP = SMINW
 506                      TJJ = SMINW
 507                      INFO = 1
 508                   END IF
 509 *
 510                   IF( XJ.EQ.ZERO )
 511      $               GO TO 70
 512 *
 513                   IF( TJJ.LT.ONE ) THEN
 514                      IF( XJ.GT.BIGNUM*TJJ ) THEN
 515                         REC = ONE / XJ
 516                         CALL SSCAL( N2, REC, X, 1 )
 517                         SCALE = SCALE*REC
 518                         XMAX = XMAX*REC
 519                      END IF
 520                   END IF
 521                   CALL SLADIV( X( J1 ), X( N+J1 ), TMP, Z, SR, SI )
 522                   X( J1 ) = SR
 523                   X( N+J1 ) = SI
 524                   XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
 525 *
 526 *                 Scale x if necessary to avoid overflow when adding a
 527 *                 multiple of column j1 of T.
 528 *
 529                   IF( XJ.GT.ONE ) THEN
 530                      REC = ONE / XJ
 531                      IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
 532                         CALL SSCAL( N2, REC, X, 1 )
 533                         SCALE = SCALE*REC
 534                      END IF
 535                   END IF
 536 *
 537                   IF( J1.GT.1 ) THEN
 538                      CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
 539                      CALL SAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
 540      $                           X( N+1 ), 1 )
 541 *
 542                      X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 )
 543                      X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 )
 544 *
 545                      XMAX = ZERO
 546                      DO 50 K = 1, J1 - 1
 547                         XMAX = MAX( XMAX, ABS( X( K ) )+
 548      $                         ABS( X( K+N ) ) )
 549    50                CONTINUE
 550                   END IF
 551 *
 552                ELSE
 553 *
 554 *                 Meet 2 by 2 diagonal block
 555 *
 556                   D( 1, 1 ) = X( J1 )
 557                   D( 2, 1 ) = X( J2 )
 558                   D( 1, 2 ) = X( N+J1 )
 559                   D( 2, 2 ) = X( N+J2 )
 560                   CALL SLALN2( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ),
 561      $                         LDT, ONE, ONE, D, 2, ZERO, -W, V, 2,
 562      $                         SCALOC, XNORM, IERR )
 563                   IF( IERR.NE.0 )
 564      $               INFO = 2
 565 *
 566                   IF( SCALOC.NE.ONE ) THEN
 567                      CALL SSCAL( 2*N, SCALOC, X, 1 )
 568                      SCALE = SCALOC*SCALE
 569                   END IF
 570                   X( J1 ) = V( 1, 1 )
 571                   X( J2 ) = V( 2, 1 )
 572                   X( N+J1 ) = V( 1, 2 )
 573                   X( N+J2 ) = V( 2, 2 )
 574 *
 575 *                 Scale X(J1), .... to avoid overflow in
 576 *                 updating right hand side.
 577 *
 578                   XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ),
 579      $                 ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) )
 580                   IF( XJ.GT.ONE ) THEN
 581                      REC = ONE / XJ
 582                      IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
 583      $                   ( BIGNUM-XMAX )*REC ) THEN
 584                         CALL SSCAL( N2, REC, X, 1 )
 585                         SCALE = SCALE*REC
 586                      END IF
 587                   END IF
 588 *
 589 *                 Update the right-hand side.
 590 *
 591                   IF( J1.GT.1 ) THEN
 592                      CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
 593                      CALL SAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
 594 *
 595                      CALL SAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
 596      $                           X( N+1 ), 1 )
 597                      CALL SAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1,
 598      $                           X( N+1 ), 1 )
 599 *
 600                      X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) +
 601      $                        B( J2 )*X( N+J2 )
 602                      X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) -
 603      $                          B( J2 )*X( J2 )
 604 *
 605                      XMAX = ZERO
 606                      DO 60 K = 1, J1 - 1
 607                         XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ),
 608      $                         XMAX )
 609    60                CONTINUE
 610                   END IF
 611 *
 612                END IF
 613    70       CONTINUE
 614 *
 615          ELSE
 616 *
 617 *           Solve (T + iB)**T*(p+iq) = c+id
 618 *
 619             JNEXT = 1
 620             DO 80 J = 1, N
 621                IF( J.LT.JNEXT )
 622      $            GO TO 80
 623                J1 = J
 624                J2 = J
 625                JNEXT = J + 1
 626                IF( J.LT.N ) THEN
 627                   IF( T( J+1, J ).NE.ZERO ) THEN
 628                      J2 = J + 1
 629                      JNEXT = J + 2
 630                   END IF
 631                END IF
 632 *
 633                IF( J1.EQ.J2 ) THEN
 634 *
 635 *                 1 by 1 diagonal block
 636 *
 637 *                 Scale if necessary to avoid overflow in forming the
 638 *                 right-hand side element by inner product.
 639 *
 640                   XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
 641                   IF( XMAX.GT.ONE ) THEN
 642                      REC = ONE / XMAX
 643                      IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
 644                         CALL SSCAL( N2, REC, X, 1 )
 645                         SCALE = SCALE*REC
 646                         XMAX = XMAX*REC
 647                      END IF
 648                   END IF
 649 *
 650                   X( J1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X, 1 )
 651                   X( N+J1 ) = X( N+J1 ) - SDOT( J1-1, T( 1, J1 ), 1,
 652      $                        X( N+1 ), 1 )
 653                   IF( J1.GT.1 ) THEN
 654                      X( J1 ) = X( J1 ) - B( J1 )*X( N+1 )
 655                      X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 )
 656                   END IF
 657                   XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
 658 *
 659                   Z = W
 660                   IF( J1.EQ.1 )
 661      $               Z = B( 1 )
 662 *
 663 *                 Scale if necessary to avoid overflow in
 664 *                 complex division
 665 *
 666                   TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
 667                   TMP = T( J1, J1 )
 668                   IF( TJJ.LT.SMINW ) THEN
 669                      TMP = SMINW
 670                      TJJ = SMINW
 671                      INFO = 1
 672                   END IF
 673 *
 674                   IF( TJJ.LT.ONE ) THEN
 675                      IF( XJ.GT.BIGNUM*TJJ ) THEN
 676                         REC = ONE / XJ
 677                         CALL SSCAL( N2, REC, X, 1 )
 678                         SCALE = SCALE*REC
 679                         XMAX = XMAX*REC
 680                      END IF
 681                   END IF
 682                   CALL SLADIV( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI )
 683                   X( J1 ) = SR
 684                   X( J1+N ) = SI
 685                   XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX )
 686 *
 687                ELSE
 688 *
 689 *                 2 by 2 diagonal block
 690 *
 691 *                 Scale if necessary to avoid overflow in forming the
 692 *                 right-hand side element by inner product.
 693 *
 694                   XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
 695      $                 ABS( X( J2 ) )+ABS( X( N+J2 ) ) )
 696                   IF( XMAX.GT.ONE ) THEN
 697                      REC = ONE / XMAX
 698                      IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
 699      $                   ( BIGNUM-XJ ) / XMAX ) THEN
 700                         CALL SSCAL( N2, REC, X, 1 )
 701                         SCALE = SCALE*REC
 702                         XMAX = XMAX*REC
 703                      END IF
 704                   END IF
 705 *
 706                   D( 1, 1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X,
 707      $                        1 )
 708                   D( 2, 1 ) = X( J2 ) - SDOT( J1-1, T( 1, J2 ), 1, X,
 709      $                        1 )
 710                   D( 1, 2 ) = X( N+J1 ) - SDOT( J1-1, T( 1, J1 ), 1,
 711      $                        X( N+1 ), 1 )
 712                   D( 2, 2 ) = X( N+J2 ) - SDOT( J1-1, T( 1, J2 ), 1,
 713      $                        X( N+1 ), 1 )
 714                   D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 )
 715                   D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 )
 716                   D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 )
 717                   D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 )
 718 *
 719                   CALL SLALN2( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ),
 720      $                         LDT, ONE, ONE, D, 2, ZERO, W, V, 2,
 721      $                         SCALOC, XNORM, IERR )
 722                   IF( IERR.NE.0 )
 723      $               INFO = 2
 724 *
 725                   IF( SCALOC.NE.ONE ) THEN
 726                      CALL SSCAL( N2, SCALOC, X, 1 )
 727                      SCALE = SCALOC*SCALE
 728                   END IF
 729                   X( J1 ) = V( 1, 1 )
 730                   X( J2 ) = V( 2, 1 )
 731                   X( N+J1 ) = V( 1, 2 )
 732                   X( N+J2 ) = V( 2, 2 )
 733                   XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
 734      $                   ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX )
 735 *
 736                END IF
 737 *
 738    80       CONTINUE
 739 *
 740          END IF
 741 *
 742       END IF
 743 *
 744       RETURN
 745 *
 746 *     End of SLAQTR
 747 *
 748       END