164
  165
  166
  167
  168
  169
  170      LOGICAL            LREAL, LTRAN
  171      INTEGER            INFO, LDT, N
  172      REAL               SCALE, W
  173
  174
  175      REAL               B( * ), T( LDT, * ), WORK( * ), X( * )
  176
  177
  178
  179
  180
  181      REAL               ZERO, ONE
  182      parameter( zero = 0.0e+0, one = 1.0e+0 )
  183
  184
  185      LOGICAL            NOTRAN
  186      INTEGER            I, IERR, J, J1, J2, JNEXT, K, N1, N2
  187      REAL               BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
  188     $                   SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
  189
  190
  191      REAL               D( 2, 2 ), V( 2, 2 )
  192
  193
  194      INTEGER            ISAMAX
  195      REAL               SASUM, SDOT, SLAMCH, SLANGE
  197
  198
  200
  201
  202      INTRINSIC          abs, max
  203
  204
  205
  206
  207
  208      notran = .NOT.ltran
  209      info = 0
  210
  211
  212
  213      IF( n.EQ.0 )
  214     $   RETURN
  215
  216
  217
  219      smlnum = 
slamch( 
'S' ) / eps
 
  220      bignum = one / smlnum
  221
  222      xnorm = 
slange( 
'M', n, n, t, ldt, d )
 
  223      IF( .NOT.lreal )
  224     $   xnorm = max( xnorm, abs( w ), 
slange( 
'M', n, 1, b, n, d ) )
 
  225      smin = max( smlnum, eps*xnorm )
  226
  227
  228
  229
  230      work( 1 ) = zero
  231      DO 10 j = 2, n
  232         work( j ) = 
sasum( j-1, t( 1, j ), 1 )
 
  233   10 CONTINUE
  234
  235      IF( .NOT.lreal ) THEN
  236         DO 20 i = 2, n
  237            work( i ) = work( i ) + abs( b( i ) )
  238   20    CONTINUE
  239      END IF
  240
  241      n2 = 2*n
  242      n1 = n
  243      IF( .NOT.lreal )
  244     $   n1 = n2
  246      xmax = abs( x( k ) )
  247      scale = one
  248
  249      IF( xmax.GT.bignum ) THEN
  250         scale = bignum / xmax
  251         CALL sscal( n1, scale, x, 1 )
 
  252         xmax = bignum
  253      END IF
  254
  255      IF( lreal ) THEN
  256
  257         IF( notran ) THEN
  258
  259
  260
  261            jnext = n
  262            DO 30 j = n, 1, -1
  263               IF( j.GT.jnext )
  264     $            GO TO 30
  265               j1 = j
  266               j2 = j
  267               jnext = j - 1
  268               IF( j.GT.1 ) THEN
  269                  IF( t( j, j-1 ).NE.zero ) THEN
  270                     j1 = j - 1
  271                     jnext = j - 2
  272                  END IF
  273               END IF
  274
  275               IF( j1.EQ.j2 ) THEN
  276
  277
  278
  279
  280
  281
  282                  xj = abs( x( j1 ) )
  283                  tjj = abs( t( j1, j1 ) )
  284                  tmp = t( j1, j1 )
  285                  IF( tjj.LT.smin ) THEN
  286                     tmp = smin
  287                     tjj = smin
  288                     info = 1
  289                  END IF
  290
  291                  IF( xj.EQ.zero )
  292     $               GO TO 30
  293
  294                  IF( tjj.LT.one ) THEN
  295                     IF( xj.GT.bignum*tjj ) THEN
  296                        rec = one / xj
  297                        CALL sscal( n, rec, x, 1 )
 
  298                        scale = scale*rec
  299                        xmax = xmax*rec
  300                     END IF
  301                  END IF
  302                  x( j1 ) = x( j1 ) / tmp
  303                  xj = abs( x( j1 ) )
  304
  305
  306
  307
  308                  IF( xj.GT.one ) THEN
  309                     rec = one / xj
  310                     IF( work( j1 ).GT.( bignum-xmax )*rec ) THEN
  311                        CALL sscal( n, rec, x, 1 )
 
  312                        scale = scale*rec
  313                     END IF
  314                  END IF
  315                  IF( j1.GT.1 ) THEN
  316                     CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x,
 
  317     $                           1 )
  319                     xmax = abs( x( k ) )
  320                  END IF
  321
  322               ELSE
  323
  324
  325
  326
  327
  328
  329                  d( 1, 1 ) = x( j1 )
  330                  d( 2, 1 ) = x( j2 )
  331                  CALL slaln2( .false., 2, 1, smin, one, t( j1, j1 ),
 
  332     $                         ldt, one, one, d, 2, zero, zero, v, 2,
  333     $                         scaloc, xnorm, ierr )
  334                  IF( ierr.NE.0 )
  335     $               info = 2
  336
  337                  IF( scaloc.NE.one ) THEN
  338                     CALL sscal( n, scaloc, x, 1 )
 
  339                     scale = scale*scaloc
  340                  END IF
  341                  x( j1 ) = v( 1, 1 )
  342                  x( j2 ) = v( 2, 1 )
  343
  344
  345
  346
  347                  xj = max( abs( v( 1, 1 ) ), abs( v( 2, 1 ) ) )
  348                  IF( xj.GT.one ) THEN
  349                     rec = one / xj
  350                     IF( max( work( j1 ), work( j2 ) ).GT.
  351     $                   ( bignum-xmax )*rec ) THEN
  352                        CALL sscal( n, rec, x, 1 )
 
  353                        scale = scale*rec
  354                     END IF
  355                  END IF
  356
  357
  358
  359                  IF( j1.GT.1 ) THEN
  360                     CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x,
 
  361     $                           1 )
  362                     CALL saxpy( j1-1, -x( j2 ), t( 1, j2 ), 1, x,
 
  363     $                           1 )
  365                     xmax = abs( x( k ) )
  366                  END IF
  367
  368               END IF
  369
  370   30       CONTINUE
  371
  372         ELSE
  373
  374
  375
  376            jnext = 1
  377            DO 40 j = 1, n
  378               IF( j.LT.jnext )
  379     $            GO TO 40
  380               j1 = j
  381               j2 = j
  382               jnext = j + 1
  383               IF( j.LT.n ) THEN
  384                  IF( t( j+1, j ).NE.zero ) THEN
  385                     j2 = j + 1
  386                     jnext = j + 2
  387                  END IF
  388               END IF
  389
  390               IF( j1.EQ.j2 ) THEN
  391
  392
  393
  394
  395
  396
  397                  xj = abs( x( j1 ) )
  398                  IF( xmax.GT.one ) THEN
  399                     rec = one / xmax
  400                     IF( work( j1 ).GT.( bignum-xj )*rec ) THEN
  401                        CALL sscal( n, rec, x, 1 )
 
  402                        scale = scale*rec
  403                        xmax = xmax*rec
  404                     END IF
  405                  END IF
  406
  407                  x( j1 ) = x( j1 ) - 
sdot( j1-1, t( 1, j1 ), 1, x,
 
  408     $               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
  433
  434
  435
  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
  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
  495
  496
  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
  527
  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,
 
  539     $                           1 )
  540                     CALL saxpy( j1-1, -x( n+j1 ), t( 1, j1 ), 1,
 
  541     $                           x( n+1 ), 1 )
  542
  543                     x( 1 ) = x( 1 ) + b( j1 )*x( n+j1 )
  544                     x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 )
  545
  546                     xmax = zero
  547                     DO 50 k = 1, j1 - 1
  548                        xmax = max( xmax, abs( x( k ) )+
  549     $                         abs( x( k+n ) ) )
  550   50                CONTINUE
  551                  END IF
  552
  553               ELSE
  554
  555
  556
  557                  d( 1, 1 ) = x( j1 )
  558                  d( 2, 1 ) = x( j2 )
  559                  d( 1, 2 ) = x( n+j1 )
  560                  d( 2, 2 ) = x( n+j2 )
  561                  CALL slaln2( .false., 2, 2, sminw, one, t( j1,
 
  562     $                         j1 ),
  563     $                         ldt, one, one, d, 2, zero, -w, v, 2,
  564     $                         scaloc, xnorm, ierr )
  565                  IF( ierr.NE.0 )
  566     $               info = 2
  567
  568                  IF( scaloc.NE.one ) THEN
  569                     CALL sscal( 2*n, scaloc, x, 1 )
 
  570                     scale = scaloc*scale
  571                  END IF
  572                  x( j1 ) = v( 1, 1 )
  573                  x( j2 ) = v( 2, 1 )
  574                  x( n+j1 ) = v( 1, 2 )
  575                  x( n+j2 ) = v( 2, 2 )
  576
  577
  578
  579
  580                  xj = max( abs( v( 1, 1 ) )+abs( v( 1, 2 ) ),
  581     $                 abs( v( 2, 1 ) )+abs( v( 2, 2 ) ) )
  582                  IF( xj.GT.one ) THEN
  583                     rec = one / xj
  584                     IF( max( work( j1 ), work( j2 ) ).GT.
  585     $                   ( bignum-xmax )*rec ) THEN
  586                        CALL sscal( n2, rec, x, 1 )
 
  587                        scale = scale*rec
  588                     END IF
  589                  END IF
  590
  591
  592
  593                  IF( j1.GT.1 ) THEN
  594                     CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x,
 
  595     $                           1 )
  596                     CALL saxpy( j1-1, -x( j2 ), t( 1, j2 ), 1, x,
 
  597     $                           1 )
  598
  599                     CALL saxpy( j1-1, -x( n+j1 ), t( 1, j1 ), 1,
 
  600     $                           x( n+1 ), 1 )
  601                     CALL saxpy( j1-1, -x( n+j2 ), t( 1, j2 ), 1,
 
  602     $                           x( n+1 ), 1 )
  603
  604                     x( 1 ) = x( 1 ) + b( j1 )*x( n+j1 ) +
  605     $                        b( j2 )*x( n+j2 )
  606                     x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) -
  607     $                          b( j2 )*x( j2 )
  608
  609                     xmax = zero
  610                     DO 60 k = 1, j1 - 1
  611                        xmax = max( abs( x( k ) )+abs( x( k+n ) ),
  612     $                         xmax )
  613   60                CONTINUE
  614                  END IF
  615
  616               END IF
  617   70       CONTINUE
  618
  619         ELSE
  620
  621
  622
  623            jnext = 1
  624            DO 80 j = 1, n
  625               IF( j.LT.jnext )
  626     $            GO TO 80
  627               j1 = j
  628               j2 = j
  629               jnext = j + 1
  630               IF( j.LT.n ) THEN
  631                  IF( t( j+1, j ).NE.zero ) THEN
  632                     j2 = j + 1
  633                     jnext = j + 2
  634                  END IF
  635               END IF
  636
  637               IF( j1.EQ.j2 ) THEN
  638
  639
  640
  641
  642
  643
  644                  xj = abs( x( j1 ) ) + abs( x( j1+n ) )
  645                  IF( xmax.GT.one ) THEN
  646                     rec = one / xmax
  647                     IF( work( j1 ).GT.( bignum-xj )*rec ) THEN
  648                        CALL sscal( n2, rec, x, 1 )
 
  649                        scale = scale*rec
  650                        xmax = xmax*rec
  651                     END IF
  652                  END IF
  653
  654                  x( j1 ) = x( j1 ) - 
sdot( j1-1, t( 1, j1 ), 1, x,
 
  655     $               1 )
  656                  x( n+j1 ) = x( n+j1 ) - 
sdot( j1-1, t( 1, j1 ), 1,
 
  657     $                        x( n+1 ), 1 )
  658                  IF( j1.GT.1 ) THEN
  659                     x( j1 ) = x( j1 ) - b( j1 )*x( n+1 )
  660                     x( n+j1 ) = x( n+j1 ) + b( j1 )*x( 1 )
  661                  END IF
  662                  xj = abs( x( j1 ) ) + abs( x( j1+n ) )
  663
  664                  z = w
  665                  IF( j1.EQ.1 )
  666     $               z = b( 1 )
  667
  668
  669
  670
  671                  tjj = abs( t( j1, j1 ) ) + abs( z )
  672                  tmp = t( j1, j1 )
  673                  IF( tjj.LT.sminw ) THEN
  674                     tmp = sminw
  675                     tjj = sminw
  676                     info = 1
  677                  END IF
  678
  679                  IF( tjj.LT.one ) THEN
  680                     IF( xj.GT.bignum*tjj ) THEN
  681                        rec = one / xj
  682                        CALL sscal( n2, rec, x, 1 )
 
  683                        scale = scale*rec
  684                        xmax = xmax*rec
  685                     END IF
  686                  END IF
  687                  CALL sladiv( x( j1 ), x( n+j1 ), tmp, -z, sr, si )
 
  688                  x( j1 ) = sr
  689                  x( j1+n ) = si
  690                  xmax = max( abs( x( j1 ) )+abs( x( j1+n ) ), xmax )
  691
  692               ELSE
  693
  694
  695
  696
  697
  698
  699                  xj = max( abs( x( j1 ) )+abs( x( n+j1 ) ),
  700     $                 abs( x( j2 ) )+abs( x( n+j2 ) ) )
  701                  IF( xmax.GT.one ) THEN
  702                     rec = one / xmax
  703                     IF( max( work( j1 ), work( j2 ) ).GT.
  704     $                   ( bignum-xj ) / xmax ) THEN
  705                        CALL sscal( n2, rec, x, 1 )
 
  706                        scale = scale*rec
  707                        xmax = xmax*rec
  708                     END IF
  709                  END IF
  710
  711                  d( 1, 1 ) = x( j1 ) - 
sdot( j1-1, t( 1, j1 ), 1, x,
 
  712     $                        1 )
  713                  d( 2, 1 ) = x( j2 ) - 
sdot( j1-1, t( 1, j2 ), 1, x,
 
  714     $                        1 )
  715                  d( 1, 2 ) = x( n+j1 ) - 
sdot( j1-1, t( 1, j1 ), 1,
 
  716     $                        x( n+1 ), 1 )
  717                  d( 2, 2 ) = x( n+j2 ) - 
sdot( j1-1, t( 1, j2 ), 1,
 
  718     $                        x( n+1 ), 1 )
  719                  d( 1, 1 ) = d( 1, 1 ) - b( j1 )*x( n+1 )
  720                  d( 2, 1 ) = d( 2, 1 ) - b( j2 )*x( n+1 )
  721                  d( 1, 2 ) = d( 1, 2 ) + b( j1 )*x( 1 )
  722                  d( 2, 2 ) = d( 2, 2 ) + b( j2 )*x( 1 )
  723
  724                  CALL slaln2( .true., 2, 2, sminw, one, t( j1, j1 ),
 
  725     $                         ldt, one, one, d, 2, zero, w, v, 2,
  726     $                         scaloc, xnorm, ierr )
  727                  IF( ierr.NE.0 )
  728     $               info = 2
  729
  730                  IF( scaloc.NE.one ) THEN
  731                     CALL sscal( n2, scaloc, x, 1 )
 
  732                     scale = scaloc*scale
  733                  END IF
  734                  x( j1 ) = v( 1, 1 )
  735                  x( j2 ) = v( 2, 1 )
  736                  x( n+j1 ) = v( 1, 2 )
  737                  x( n+j2 ) = v( 2, 2 )
  738                  xmax = max( abs( x( j1 ) )+abs( x( n+j1 ) ),
  739     $                   abs( x( j2 ) )+abs( x( n+j2 ) ), xmax )
  740
  741               END IF
  742
  743   80       CONTINUE
  744
  745         END IF
  746
  747      END IF
  748
  749      RETURN
  750
  751
  752
real function sasum(n, sx, incx)
SASUM
 
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
 
real function sdot(n, sx, incx, sy, incy)
SDOT
 
integer function isamax(n, sx, incx)
ISAMAX
 
subroutine sladiv(a, b, c, d, p, q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
 
subroutine slaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
 
real function slamch(cmach)
SLAMCH
 
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
 
subroutine sscal(n, sa, sx, incx)
SSCAL