182
  183
  184
  185
  186
  187
  188      CHARACTER          UPLO
  189      INTEGER            INFO, KB, LDA, LDW, N, NB
  190
  191
  192      INTEGER            IPIV( * )
  193      REAL               A( LDA, * ), W( LDW, * )
  194
  195
  196
  197
  198
  199      REAL               ZERO, ONE
  200      parameter( zero = 0.0e+0, one = 1.0e+0 )
  201      REAL               EIGHT, SEVTEN
  202      parameter( eight = 8.0e+0, sevten = 17.0e+0 )
  203
  204
  205      LOGICAL            DONE
  206      INTEGER            IMAX, ITEMP, J, JJ, JMAX, JP1, JP2, K, KK,
  207     $                   KW, KKW, KP, KSTEP, P, II
  208 
  209      REAL               ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
  210     $                   STEMP, R1, ROWMAX, T, SFMIN
  211
  212
  213      LOGICAL            LSAME
  214      INTEGER            ISAMAX
  215      REAL               SLAMCH
  217
  218
  220
  221
  222      INTRINSIC          abs, max, min, sqrt
  223
  224
  225
  226      info = 0
  227
  228
  229
  230      alpha = ( one+sqrt( sevten ) ) / eight
  231
  232
  233
  235
  236      IF( 
lsame( uplo, 
'U' ) ) 
THEN 
  237
  238
  239
  240
  241
  242
  243
  244         k = n
  245   10    CONTINUE
  246
  247
  248
  249         kw = nb + k - n
  250
  251
  252
  253         IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
  254     $      GO TO 30
  255
  256         kstep = 1
  257         p = k
  258
  259
  260
  261         CALL scopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
 
  262         IF( k.LT.n )
  263     $      
CALL sgemv( 
'No transpose', k, n-k, -one, a( 1, k+1 ),
 
  264     $                  lda, w( k, kw+1 ), ldw, one, w( 1, kw ), 1 )
  265
  266
  267
  268
  269         absakk = abs( w( k, kw ) )
  270
  271
  272
  273
  274
  275         IF( k.GT.1 ) THEN
  276            imax = 
isamax( k-1, w( 1, kw ), 1 )
 
  277            colmax = abs( w( imax, kw ) )
  278         ELSE
  279            colmax = zero
  280         END IF
  281
  282         IF( max( absakk, colmax ).EQ.zero ) THEN
  283
  284
  285
  286            IF( info.EQ.0 )
  287     $         info = k
  288            kp = k
  289            CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
 
  290         ELSE
  291
  292
  293
  294
  295
  296
  297
  298
  299            IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
  300
  301
  302
  303               kp = k
  304
  305            ELSE
  306
  307               done = .false.
  308
  309
  310
  311   12          CONTINUE
  312
  313
  314
  315
  316
  317
  318                  CALL scopy( imax, a( 1, imax ), 1, w( 1, kw-1 ),
 
  319     $                        1 )
  320                  CALL scopy( k-imax, a( imax, imax+1 ), lda,
 
  321     $                        w( imax+1, kw-1 ), 1 )
  322
  323                  IF( k.LT.n )
  324     $               
CALL sgemv( 
'No transpose', k, n-k, -one,
 
  325     $                           a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
  326     $                           one, w( 1, kw-1 ), 1 )
  327
  328
  329
  330
  331
  332                  IF( imax.NE.k ) THEN
  333                     jmax = imax + 
isamax( k-imax, w( imax+1, kw-1 ),
 
  334     $                                     1 )
  335                     rowmax = abs( w( jmax, kw-1 ) )
  336                  ELSE
  337                     rowmax = zero
  338                  END IF
  339
  340                  IF( imax.GT.1 ) THEN
  341                     itemp = 
isamax( imax-1, w( 1, kw-1 ), 1 )
 
  342                     stemp = abs( w( itemp, kw-1 ) )
  343                     IF( stemp.GT.rowmax ) THEN
  344                        rowmax = stemp
  345                        jmax = itemp
  346                     END IF
  347                  END IF
  348
  349
  350
  351
  352
  353                  IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
  354     $            THEN
  355
  356
  357
  358
  359                     kp = imax
  360
  361
  362
  363                     CALL scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
 
  364
  365                     done = .true.
  366
  367
  368
  369
  370                  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
  371     $            THEN
  372
  373
  374
  375
  376                     kp = imax
  377                     kstep = 2
  378                     done = .true.
  379                  ELSE
  380
  381
  382
  383                     p = imax
  384                     colmax = rowmax
  385                     imax = jmax
  386
  387
  388
  389                     CALL scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
 
  390
  391                  END IF
  392
  393
  394
  395               IF( .NOT. done ) GOTO 12
  396
  397            END IF
  398
  399
  400
  401            kk = k - kstep + 1
  402
  403
  404
  405            kkw = nb + kk - n
  406
  407            IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
  408
  409
  410
  411               CALL scopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
 
  412               CALL scopy( p, a( 1, k ), 1, a( 1, p ), 1 )
 
  413
  414
  415
  416
  417               CALL sswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
 
  418               CALL sswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
 
  419     $                     ldw )
  420            END IF
  421
  422
  423
  424            IF( kp.NE.kk ) THEN
  425
  426
  427
  428               a( kp, k ) = a( kk, k )
  429               CALL scopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
 
  430     $                     lda )
  431               CALL scopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
 
  432
  433
  434
  435
  436               CALL sswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ),
 
  437     $                     lda )
  438               CALL sswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
 
  439     $                     ldw )
  440            END IF
  441
  442            IF( kstep.EQ.1 ) THEN
  443
  444
  445
  446
  447
  448
  449
  450
  451
  452               CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
 
  453               IF( k.GT.1 ) THEN
  454                  IF( abs( a( k, k ) ).GE.sfmin ) THEN
  455                     r1 = one / a( k, k )
  456                     CALL sscal( k-1, r1, a( 1, k ), 1 )
 
  457                  ELSE IF( a( k, k ).NE.zero ) THEN
  458                     DO 14 ii = 1, k - 1
  459                        a( ii, k ) = a( ii, k ) / a( k, k )
  460   14                CONTINUE
  461                  END IF
  462               END IF
  463
  464            ELSE
  465
  466
  467
  468
  469
  470
  471
  472
  473
  474               IF( k.GT.2 ) THEN
  475
  476
  477
  478                  d12 = w( k-1, kw )
  479                  d11 = w( k, kw ) / d12
  480                  d22 = w( k-1, kw-1 ) / d12
  481                  t = one / ( d11*d22-one )
  482                  DO 20 j = 1, k - 2
  483                     a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
  484     $                             d12 )
  485                     a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
  486     $                           d12 )
  487   20             CONTINUE
  488               END IF
  489
  490
  491
  492               a( k-1, k-1 ) = w( k-1, kw-1 )
  493               a( k-1, k ) = w( k-1, kw )
  494               a( k, k ) = w( k, kw )
  495            END IF
  496         END IF
  497
  498
  499
  500         IF( kstep.EQ.1 ) THEN
  501            ipiv( k ) = kp
  502         ELSE
  503            ipiv( k ) = -p
  504            ipiv( k-1 ) = -kp
  505         END IF
  506
  507
  508
  509         k = k - kstep
  510         GO TO 10
  511
  512   30    CONTINUE
  513
  514
  515
  516
  517
  518         CALL sgemmtr( 
'Upper', 
'No transpose', 
'Transpose', k, n-k,
 
  519     $                 -one, a( 1, k+1 ), lda, w( 1, kw+1 ), ldw,
  520     $                 one, a( 1, 1 ), lda )
  521
  522
  523
  524
  525         j = k + 1
  526   60    CONTINUE
  527
  528            kstep = 1
  529            jp1 = 1
  530            jj = j
  531            jp2 = ipiv( j )
  532            IF( jp2.LT.0 ) THEN
  533               jp2 = -jp2
  534               j = j + 1
  535               jp1 = -ipiv( j )
  536               kstep = 2
  537            END IF
  538
  539            j = j + 1
  540            IF( jp2.NE.jj .AND. j.LE.n )
  541     $         
CALL sswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
 
  542            jj = j - 1
  543            IF( jp1.NE.jj .AND. kstep.EQ.2 )
  544     $         
CALL sswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
 
  545         IF( j.LE.n )
  546     $      GO TO 60
  547
  548
  549
  550         kb = n - k
  551
  552      ELSE
  553
  554
  555
  556
  557
  558
  559
  560         k = 1
  561   70   CONTINUE
  562
  563
  564
  565         IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
  566     $      GO TO 90
  567
  568         kstep = 1
  569         p = k
  570
  571
  572
  573         CALL scopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
 
  574         IF( k.GT.1 )
  575     $      
CALL sgemv( 
'No transpose', n-k+1, k-1, -one, a( k, 1 ),
 
  576     $                  lda, w( k, 1 ), ldw, one, w( k, k ), 1 )
  577
  578
  579
  580
  581         absakk = abs( w( k, k ) )
  582
  583
  584
  585
  586
  587         IF( k.LT.n ) THEN
  588            imax = k + 
isamax( n-k, w( k+1, k ), 1 )
 
  589            colmax = abs( w( imax, k ) )
  590         ELSE
  591            colmax = zero
  592         END IF
  593
  594         IF( max( absakk, colmax ).EQ.zero ) THEN
  595
  596
  597
  598            IF( info.EQ.0 )
  599     $         info = k
  600            kp = k
  601            CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
 
  602         ELSE
  603
  604
  605
  606
  607
  608
  609
  610
  611            IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
  612
  613
  614
  615               kp = k
  616
  617            ELSE
  618
  619               done = .false.
  620
  621
  622
  623   72          CONTINUE
  624
  625
  626
  627
  628
  629
  630                  CALL scopy( imax-k, a( imax, k ), lda, w( k, k+1 ),
 
  631     $                        1)
  632                  CALL scopy( n-imax+1, a( imax, imax ), 1,
 
  633     $                        w( imax, k+1 ), 1 )
  634                  IF( k.GT.1 )
  635     $               
CALL sgemv( 
'No transpose', n-k+1, k-1, -one,
 
  636     $                           a( k, 1 ), lda, w( imax, 1 ), ldw,
  637     $                           one, w( k, k+1 ), 1 )
  638
  639
  640
  641
  642
  643                  IF( imax.NE.k ) THEN
  644                     jmax = k - 1 + 
isamax( imax-k, w( k, k+1 ), 1 )
 
  645                     rowmax = abs( w( jmax, k+1 ) )
  646                  ELSE
  647                     rowmax = zero
  648                  END IF
  649
  650                  IF( imax.LT.n ) THEN
  651                     itemp = imax + 
isamax( n-imax, w( imax+1, k+1 ),
 
  652     $                                      1)
  653                     stemp = abs( w( itemp, k+1 ) )
  654                     IF( stemp.GT.rowmax ) THEN
  655                        rowmax = stemp
  656                        jmax = itemp
  657                     END IF
  658                  END IF
  659
  660
  661
  662
  663
  664                  IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
  665     $            THEN
  666
  667
  668
  669
  670                     kp = imax
  671
  672
  673
  674                     CALL scopy( n-k+1, w( k, k+1 ), 1, w( k, k ),
 
  675     $                           1 )
  676
  677                     done = .true.
  678
  679
  680
  681
  682                  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
  683     $            THEN
  684
  685
  686
  687
  688                     kp = imax
  689                     kstep = 2
  690                     done = .true.
  691                  ELSE
  692
  693
  694
  695                     p = imax
  696                     colmax = rowmax
  697                     imax = jmax
  698
  699
  700
  701                     CALL scopy( n-k+1, w( k, k+1 ), 1, w( k, k ),
 
  702     $                           1 )
  703
  704                  END IF
  705
  706
  707
  708               IF( .NOT. done ) GOTO 72
  709
  710            END IF
  711
  712
  713
  714            kk = k + kstep - 1
  715
  716            IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
  717
  718
  719
  720               CALL scopy( p-k, a( k, k ), 1, a( p, k ), lda )
 
  721               CALL scopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
 
  722
  723
  724
  725
  726               CALL sswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
 
  727               CALL sswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
 
  728            END IF
  729
  730
  731
  732            IF( kp.NE.kk ) THEN
  733
  734
  735
  736               a( kp, k ) = a( kk, k )
  737               CALL scopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ),
 
  738     $                     lda )
  739               CALL scopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
 
  740
  741
  742
  743               CALL sswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
 
  744               CALL sswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
 
  745            END IF
  746
  747            IF( kstep.EQ.1 ) THEN
  748
  749
  750
  751
  752
  753
  754
  755
  756
  757               CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
 
  758               IF( k.LT.n ) THEN
  759                  IF( abs( a( k, k ) ).GE.sfmin ) THEN
  760                     r1 = one / a( k, k )
  761                     CALL sscal( n-k, r1, a( k+1, k ), 1 )
 
  762                  ELSE IF( a( k, k ).NE.zero ) THEN
  763                     DO 74 ii = k + 1, n
  764                        a( ii, k ) = a( ii, k ) / a( k, k )
  765   74                CONTINUE
  766                  END IF
  767               END IF
  768
  769            ELSE
  770
  771
  772
  773
  774
  775
  776
  777
  778               IF( k.LT.n-1 ) THEN
  779
  780
  781
  782                  d21 = w( k+1, k )
  783                  d11 = w( k+1, k+1 ) / d21
  784                  d22 = w( k, k ) / d21
  785                  t = one / ( d11*d22-one )
  786                  DO 80 j = k + 2, n
  787                     a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
  788     $                           d21 )
  789                     a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
  790     $                             d21 )
  791   80             CONTINUE
  792               END IF
  793
  794
  795
  796               a( k, k ) = w( k, k )
  797               a( k+1, k ) = w( k+1, k )
  798               a( k+1, k+1 ) = w( k+1, k+1 )
  799            END IF
  800         END IF
  801
  802
  803
  804         IF( kstep.EQ.1 ) THEN
  805            ipiv( k ) = kp
  806         ELSE
  807            ipiv( k ) = -p
  808            ipiv( k+1 ) = -kp
  809         END IF
  810
  811
  812
  813         k = k + kstep
  814         GO TO 70
  815
  816   90    CONTINUE
  817
  818
  819
  820
  821
  822         CALL sgemmtr( 
'Lower', 
'No transpose', 
'Transpose', n-k+1,
 
  823     $                 k-1, -one, a( k, 1 ), lda, w( k, 1 ), ldw,
  824     $                 one, a( k, k ), lda )
  825
  826
  827
  828
  829         j = k - 1
  830  120    CONTINUE
  831
  832            kstep = 1
  833            jp1 = 1
  834            jj = j
  835            jp2 = ipiv( j )
  836            IF( jp2.LT.0 ) THEN
  837               jp2 = -jp2
  838               j = j - 1
  839               jp1 = -ipiv( j )
  840               kstep = 2
  841            END IF
  842
  843            j = j - 1
  844            IF( jp2.NE.jj .AND. j.GE.1 )
  845     $         
CALL sswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
 
  846            jj = j + 1
  847            IF( jp1.NE.jj .AND. kstep.EQ.2 )
  848     $         
CALL sswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
 
  849         IF( j.GE.1 )
  850     $      GO TO 120
  851
  852
  853
  854         kb = k - 1
  855
  856      END IF
  857      RETURN
  858
  859
  860
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
subroutine sgemmtr(uplo, transa, transb, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMMTR
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
integer function isamax(n, sx, incx)
ISAMAX
real function slamch(cmach)
SLAMCH
logical function lsame(ca, cb)
LSAME
subroutine sscal(n, sa, sx, incx)
SSCAL
subroutine sswap(n, sx, incx, sy, incy)
SSWAP