262
263
264
265
266
267
268 CHARACTER UPLO
269 INTEGER INFO, KB, LDA, LDW, N, NB
270
271
272 INTEGER IPIV( * )
273 COMPLEX A( LDA, * ), W( LDW, * ), E( * )
274
275
276
277
278
279 REAL ZERO, ONE
280 parameter( zero = 0.0e+0, one = 1.0e+0 )
281 REAL EIGHT, SEVTEN
282 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
283 COMPLEX CONE, CZERO
284 parameter( cone = ( 1.0e+0, 0.0e+0 ),
285 $ czero = ( 0.0e+0, 0.0e+0 ) )
286
287
288 LOGICAL DONE
289 INTEGER IMAX, ITEMP, II, J, JB, JJ, JMAX, K, KK, KKW,
290 $ KP, KSTEP, KW, P
291 REAL ABSAKK, ALPHA, COLMAX, STEMP, R1, ROWMAX, T,
292 $ SFMIN
293 COMPLEX D11, D21, D22, Z
294
295
296 LOGICAL LSAME
297 INTEGER ICAMAX
298 REAL SLAMCH
300
301
303
304
305 INTRINSIC abs, conjg, aimag, max, min, real, sqrt
306
307
308 REAL CABS1
309
310
311 cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
312
313
314
315 info = 0
316
317
318
319 alpha = ( one+sqrt( sevten ) ) / eight
320
321
322
324
325 IF(
lsame( uplo,
'U' ) )
THEN
326
327
328
329
330
331
332
333
334 e( 1 ) = czero
335
336
337
338 k = n
339 10 CONTINUE
340
341
342
343 kw = nb + k - n
344
345
346
347 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
348 $ GO TO 30
349
350 kstep = 1
351 p = k
352
353
354
355 IF( k.GT.1 )
356 $
CALL ccopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
357 w( k, kw ) = real( a( k, k ) )
358 IF( k.LT.n ) THEN
359 CALL cgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
360 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
361 w( k, kw ) = real( w( k, kw ) )
362 END IF
363
364
365
366
367 absakk = abs( real( w( k, kw ) ) )
368
369
370
371
372
373 IF( k.GT.1 ) THEN
374 imax =
icamax( k-1, w( 1, kw ), 1 )
375 colmax = cabs1( w( imax, kw ) )
376 ELSE
377 colmax = zero
378 END IF
379
380 IF( max( absakk, colmax ).EQ.zero ) THEN
381
382
383
384 IF( info.EQ.0 )
385 $ info = k
386 kp = k
387 a( k, k ) = real( w( k, kw ) )
388 IF( k.GT.1 )
389 $
CALL ccopy( k-1, w( 1, kw ), 1, a( 1, k ), 1 )
390
391
392
393 IF( k.GT.1 )
394 $ e( k ) = czero
395
396 ELSE
397
398
399
400
401
402
403
404
405 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
406
407
408
409 kp = k
410
411 ELSE
412
413
414
415 done = .false.
416
417 12 CONTINUE
418
419
420
421
422
423
424 IF( imax.GT.1 )
425 $
CALL ccopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ),
426 $ 1 )
427 w( imax, kw-1 ) = real( a( imax, imax ) )
428
429 CALL ccopy( k-imax, a( imax, imax+1 ), lda,
430 $ w( imax+1, kw-1 ), 1 )
431 CALL clacgv( k-imax, w( imax+1, kw-1 ), 1 )
432
433 IF( k.LT.n ) THEN
434 CALL cgemv(
'No transpose', k, n-k, -cone,
435 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
436 $ cone, w( 1, kw-1 ), 1 )
437 w( imax, kw-1 ) = real( w( imax, kw-1 ) )
438 END IF
439
440
441
442
443
444 IF( imax.NE.k ) THEN
445 jmax = imax +
icamax( k-imax, w( imax+1, kw-1 ),
446 $ 1 )
447 rowmax = cabs1( w( jmax, kw-1 ) )
448 ELSE
449 rowmax = zero
450 END IF
451
452 IF( imax.GT.1 ) THEN
453 itemp =
icamax( imax-1, w( 1, kw-1 ), 1 )
454 stemp = cabs1( w( itemp, kw-1 ) )
455 IF( stemp.GT.rowmax ) THEN
456 rowmax = stemp
457 jmax = itemp
458 END IF
459 END IF
460
461
462
463
464
465
466 IF( .NOT.( abs( real( w( imax,kw-1 ) ) )
467 $ .LT.alpha*rowmax ) ) THEN
468
469
470
471
472 kp = imax
473
474
475
476 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
477
478 done = .true.
479
480
481
482
483
484 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
485 $ THEN
486
487
488
489
490 kp = imax
491 kstep = 2
492 done = .true.
493
494
495 ELSE
496
497
498
499 p = imax
500 colmax = rowmax
501 imax = jmax
502
503
504
505 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
506
507 END IF
508
509
510
511
512 IF( .NOT.done ) GOTO 12
513
514 END IF
515
516
517
518
519
520
521
522 kk = k - kstep + 1
523
524
525
526 kkw = nb + kk - n
527
528
529
530
531 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
532
533
534
535
536
537
538 a( p, p ) = real( a( k, k ) )
539 CALL ccopy( k-1-p, a( p+1, k ), 1, a( p, p+1 ),
540 $ lda )
541 CALL clacgv( k-1-p, a( p, p+1 ), lda )
542 IF( p.GT.1 )
543 $
CALL ccopy( p-1, a( 1, k ), 1, a( 1, p ), 1 )
544
545
546
547
548
549
550 IF( k.LT.n )
551 $
CALL cswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
552 $ lda )
553 CALL cswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
554 $ ldw )
555 END IF
556
557
558
559
560 IF( kp.NE.kk ) THEN
561
562
563
564
565
566
567 a( kp, kp ) = real( a( kk, kk ) )
568 CALL ccopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
569 $ lda )
570 CALL clacgv( kk-1-kp, a( kp, kp+1 ), lda )
571 IF( kp.GT.1 )
572 $
CALL ccopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
573
574
575
576
577
578
579 IF( k.LT.n )
580 $
CALL cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
581 $ lda )
582 CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
583 $ ldw )
584 END IF
585
586 IF( kstep.EQ.1 ) THEN
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604 CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
605 IF( k.GT.1 ) THEN
606
607
608
609
610
611
612
613 t = real( a( k, k ) )
614 IF( abs( t ).GE.sfmin ) THEN
615 r1 = one / t
616 CALL csscal( k-1, r1, a( 1, k ), 1 )
617 ELSE
618 DO 14 ii = 1, k-1
619 a( ii, k ) = a( ii, k ) / t
620 14 CONTINUE
621 END IF
622
623
624
625 CALL clacgv( k-1, w( 1, kw ), 1 )
626
627
628
629 e( k ) = czero
630
631 END IF
632
633 ELSE
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650 IF( k.GT.2 ) THEN
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697 d21 = w( k-1, kw )
698 d11 = w( k, kw ) / conjg( d21 )
699 d22 = w( k-1, kw-1 ) / d21
700 t = one / ( real( d11*d22 )-one )
701
702
703
704
705
706 DO 20 j = 1, k - 2
707 a( j, k-1 ) = t*( ( d11*w( j, kw-1 )-w( j, kw ) ) /
708 $ d21 )
709 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
710 $ conjg( d21 ) )
711 20 CONTINUE
712 END IF
713
714
715
716
717
718 a( k-1, k-1 ) = w( k-1, kw-1 )
719 a( k-1, k ) = czero
720 a( k, k ) = w( k, kw )
721 e( k ) = w( k-1, kw )
722 e( k-1 ) = czero
723
724
725
726 CALL clacgv( k-1, w( 1, kw ), 1 )
727 CALL clacgv( k-2, w( 1, kw-1 ), 1 )
728
729 END IF
730
731
732
733 END IF
734
735
736
737 IF( kstep.EQ.1 ) THEN
738 ipiv( k ) = kp
739 ELSE
740 ipiv( k ) = -p
741 ipiv( k-1 ) = -kp
742 END IF
743
744
745
746 k = k - kstep
747 GO TO 10
748
749 30 CONTINUE
750
751
752
753
754
755
756
757
758 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
759 jb = min( nb, k-j+1 )
760
761
762
763 DO 40 jj = j, j + jb - 1
764 a( jj, jj ) = real( a( jj, jj ) )
765 CALL cgemv(
'No transpose', jj-j+1, n-k, -cone,
766 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
767 $ a( j, jj ), 1 )
768 a( jj, jj ) = real( a( jj, jj ) )
769 40 CONTINUE
770
771
772
773 IF( j.GE.2 )
774 $
CALL cgemm(
'No transpose',
'Transpose', j-1, jb, n-k,
775 $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
776 $ cone, a( 1, j ), lda )
777 50 CONTINUE
778
779
780
781 kb = n - k
782
783 ELSE
784
785
786
787
788
789
790
791 e( n ) = czero
792
793
794
795 k = 1
796 70 CONTINUE
797
798
799
800 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
801 $ GO TO 90
802
803 kstep = 1
804 p = k
805
806
807
808 w( k, k ) = real( a( k, k ) )
809 IF( k.LT.n )
810 $
CALL ccopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
811 IF( k.GT.1 ) THEN
812 CALL cgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
813 $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
814 w( k, k ) = real( w( k, k ) )
815 END IF
816
817
818
819
820 absakk = abs( real( w( k, k ) ) )
821
822
823
824
825
826 IF( k.LT.n ) THEN
827 imax = k +
icamax( n-k, w( k+1, k ), 1 )
828 colmax = cabs1( w( imax, k ) )
829 ELSE
830 colmax = zero
831 END IF
832
833 IF( max( absakk, colmax ).EQ.zero ) THEN
834
835
836
837 IF( info.EQ.0 )
838 $ info = k
839 kp = k
840 a( k, k ) = real( w( k, k ) )
841 IF( k.LT.n )
842 $
CALL ccopy( n-k, w( k+1, k ), 1, a( k+1, k ), 1 )
843
844
845
846 IF( k.LT.n )
847 $ e( k ) = czero
848
849 ELSE
850
851
852
853
854
855
856
857
858
859 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
860
861
862
863 kp = k
864
865 ELSE
866
867 done = .false.
868
869
870
871 72 CONTINUE
872
873
874
875
876
877
878 CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
879 CALL clacgv( imax-k, w( k, k+1 ), 1 )
880 w( imax, k+1 ) = real( a( imax, imax ) )
881
882 IF( imax.LT.n )
883 $
CALL ccopy( n-imax, a( imax+1, imax ), 1,
884 $ w( imax+1, k+1 ), 1 )
885
886 IF( k.GT.1 ) THEN
887 CALL cgemv(
'No transpose', n-k+1, k-1, -cone,
888 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
889 $ cone, w( k, k+1 ), 1 )
890 w( imax, k+1 ) = real( w( imax, k+1 ) )
891 END IF
892
893
894
895
896
897 IF( imax.NE.k ) THEN
898 jmax = k - 1 +
icamax( imax-k, w( k, k+1 ), 1 )
899 rowmax = cabs1( w( jmax, k+1 ) )
900 ELSE
901 rowmax = zero
902 END IF
903
904 IF( imax.LT.n ) THEN
905 itemp = imax +
icamax( n-imax, w( imax+1, k+1 ), 1)
906 stemp = cabs1( w( itemp, k+1 ) )
907 IF( stemp.GT.rowmax ) THEN
908 rowmax = stemp
909 jmax = itemp
910 END IF
911 END IF
912
913
914
915
916
917
918 IF( .NOT.( abs( real( w( imax,k+1 ) ) )
919 $ .LT.alpha*rowmax ) ) THEN
920
921
922
923
924 kp = imax
925
926
927
928 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
929
930 done = .true.
931
932
933
934
935
936 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
937 $ THEN
938
939
940
941
942 kp = imax
943 kstep = 2
944 done = .true.
945
946
947 ELSE
948
949
950
951 p = imax
952 colmax = rowmax
953 imax = jmax
954
955
956
957 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
958
959 END IF
960
961
962
963
964 IF( .NOT.done ) GOTO 72
965
966 END IF
967
968
969
970
971
972
973
974 kk = k + kstep - 1
975
976
977
978
979 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
980
981
982
983
984
985
986 a( p, p ) = real( a( k, k ) )
987 CALL ccopy( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
988 CALL clacgv( p-k-1, a( p, k+1 ), lda )
989 IF( p.LT.n )
990 $
CALL ccopy( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
991
992
993
994
995
996
997 IF( k.GT.1 )
998 $
CALL cswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
999 CALL cswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
1000 END IF
1001
1002
1003
1004
1005 IF( kp.NE.kk ) THEN
1006
1007
1008
1009
1010
1011
1012 a( kp, kp ) = real( a( kk, kk ) )
1013 CALL ccopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
1014 $ lda )
1015 CALL clacgv( kp-kk-1, a( kp, kk+1 ), lda )
1016 IF( kp.LT.n )
1017 $
CALL ccopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
1018
1019
1020
1021
1022
1023
1024 IF( k.GT.1 )
1025 $
CALL cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
1026 CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
1027 END IF
1028
1029 IF( kstep.EQ.1 ) THEN
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047 CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
1048 IF( k.LT.n ) THEN
1049
1050
1051
1052
1053
1054
1055
1056 t = real( a( k, k ) )
1057 IF( abs( t ).GE.sfmin ) THEN
1058 r1 = one / t
1059 CALL csscal( n-k, r1, a( k+1, k ), 1 )
1060 ELSE
1061 DO 74 ii = k + 1, n
1062 a( ii, k ) = a( ii, k ) / t
1063 74 CONTINUE
1064 END IF
1065
1066
1067
1068 CALL clacgv( n-k, w( k+1, k ), 1 )
1069
1070
1071
1072 e( k ) = czero
1073
1074 END IF
1075
1076 ELSE
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093 IF( k.LT.n-1 ) THEN
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140 d21 = w( k+1, k )
1141 d11 = w( k+1, k+1 ) / d21
1142 d22 = w( k, k ) / conjg( d21 )
1143 t = one / ( real( d11*d22 )-one )
1144
1145
1146
1147
1148
1149 DO 80 j = k + 2, n
1150 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
1151 $ conjg( d21 ) )
1152 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
1153 $ d21 )
1154 80 CONTINUE
1155 END IF
1156
1157
1158
1159
1160
1161 a( k, k ) = w( k, k )
1162 a( k+1, k ) = czero
1163 a( k+1, k+1 ) = w( k+1, k+1 )
1164 e( k ) = w( k+1, k )
1165 e( k+1 ) = czero
1166
1167
1168
1169 CALL clacgv( n-k, w( k+1, k ), 1 )
1170 CALL clacgv( n-k-1, w( k+2, k+1 ), 1 )
1171
1172 END IF
1173
1174
1175
1176 END IF
1177
1178
1179
1180 IF( kstep.EQ.1 ) THEN
1181 ipiv( k ) = kp
1182 ELSE
1183 ipiv( k ) = -p
1184 ipiv( k+1 ) = -kp
1185 END IF
1186
1187
1188
1189 k = k + kstep
1190 GO TO 70
1191
1192 90 CONTINUE
1193
1194
1195
1196
1197
1198
1199
1200
1201 DO 110 j = k, n, nb
1202 jb = min( nb, n-j+1 )
1203
1204
1205
1206 DO 100 jj = j, j + jb - 1
1207 a( jj, jj ) = real( a( jj, jj ) )
1208 CALL cgemv(
'No transpose', j+jb-jj, k-1, -cone,
1209 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
1210 $ a( jj, jj ), 1 )
1211 a( jj, jj ) = real( a( jj, jj ) )
1212 100 CONTINUE
1213
1214
1215
1216 IF( j+jb.LE.n )
1217 $
CALL cgemm(
'No transpose',
'Transpose', n-j-jb+1, jb,
1218 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
1219 $ ldw, cone, a( j+jb, j ), lda )
1220 110 CONTINUE
1221
1222
1223
1224 kb = k - 1
1225
1226 END IF
1227 RETURN
1228
1229
1230
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
integer function icamax(n, cx, incx)
ICAMAX
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
real function slamch(cmach)
SLAMCH
logical function lsame(ca, cb)
LSAME
subroutine csscal(n, sa, cx, incx)
CSSCAL
subroutine cswap(n, cx, incx, cy, incy)
CSWAP