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