177
178
179
180
181
182
183 CHARACTER UPLO
184 INTEGER INFO, KB, LDA, LDW, N, NB
185
186
187 INTEGER IPIV( * )
188 COMPLEX*16 A( LDA, * ), W( LDW, * )
189
190
191
192
193
194 DOUBLE PRECISION ZERO, ONE
195 parameter( zero = 0.0d+0, one = 1.0d+0 )
196 COMPLEX*16 CONE
197 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
198 DOUBLE PRECISION EIGHT, SEVTEN
199 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
200
201
202 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
203 $ KSTEP, KW
204 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, R1, ROWMAX, T
205 COMPLEX*16 D11, D21, D22, Z
206
207
208 LOGICAL LSAME
209 INTEGER IZAMAX
211
212
214
215
216 INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
217
218
219 DOUBLE PRECISION CABS1
220
221
222 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
223
224
225
226 info = 0
227
228
229
230 alpha = ( one+sqrt( sevten ) ) / eight
231
232 IF(
lsame( uplo,
'U' ) )
THEN
233
234
235
236
237
238
239
240
241
242 k = n
243 10 CONTINUE
244 kw = nb + k - n
245
246
247
248 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
249 $ GO TO 30
250
251 kstep = 1
252
253
254
255 CALL zcopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
256 w( k, kw ) = dble( a( k, k ) )
257 IF( k.LT.n ) THEN
258 CALL zgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
259 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
260 w( k, kw ) = dble( w( k, kw ) )
261 END IF
262
263
264
265
266 absakk = abs( dble( w( k, kw ) ) )
267
268
269
270
271
272 IF( k.GT.1 ) THEN
273 imax =
izamax( k-1, w( 1, kw ), 1 )
274 colmax = cabs1( w( imax, kw ) )
275 ELSE
276 colmax = zero
277 END IF
278
279 IF( max( absakk, colmax ).EQ.zero ) THEN
280
281
282
283 IF( info.EQ.0 )
284 $ info = k
285 kp = k
286 a( k, k ) = dble( a( k, k ) )
287 ELSE
288
289
290
291
292
293
294 IF( absakk.GE.alpha*colmax ) THEN
295
296
297
298 kp = k
299 ELSE
300
301
302
303
304
305
306 CALL zcopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
307 w( imax, kw-1 ) = dble( a( imax, imax ) )
308 CALL zcopy( k-imax, a( imax, imax+1 ), lda,
309 $ w( imax+1, kw-1 ), 1 )
310 CALL zlacgv( k-imax, w( imax+1, kw-1 ), 1 )
311 IF( k.LT.n ) THEN
312 CALL zgemv(
'No transpose', k, n-k, -cone,
313 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
314 $ cone, w( 1, kw-1 ), 1 )
315 w( imax, kw-1 ) = dble( w( imax, kw-1 ) )
316 END IF
317
318
319
320
321
322 jmax = imax +
izamax( k-imax, w( imax+1, kw-1 ), 1 )
323 rowmax = cabs1( w( jmax, kw-1 ) )
324 IF( imax.GT.1 ) THEN
325 jmax =
izamax( imax-1, w( 1, kw-1 ), 1 )
326 rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
327 END IF
328
329
330 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
331
332
333
334 kp = k
335
336
337 ELSE IF( abs( dble( w( imax, kw-1 ) ) ).GE.alpha*rowmax )
338 $ THEN
339
340
341
342
343 kp = imax
344
345
346
347 CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
348
349
350 ELSE
351
352
353
354
355 kp = imax
356 kstep = 2
357 END IF
358
359
360
361
362 END IF
363
364
365
366
367
368
369
370 kk = k - kstep + 1
371
372
373
374 kkw = nb + kk - n
375
376
377
378
379 IF( kp.NE.kk ) THEN
380
381
382
383
384
385
386 a( kp, kp ) = dble( a( kk, kk ) )
387 CALL zcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
388 $ lda )
389 CALL zlacgv( kk-1-kp, a( kp, kp+1 ), lda )
390 IF( kp.GT.1 )
391 $
CALL zcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
392
393
394
395
396
397
398 IF( k.LT.n )
399 $
CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
400 $ lda )
401 CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
402 $ ldw )
403 END IF
404
405 IF( kstep.EQ.1 ) THEN
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423 CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
424 IF( k.GT.1 ) THEN
425
426
427
428
429
430 r1 = one / dble( a( k, k ) )
431 CALL zdscal( k-1, r1, a( 1, k ), 1 )
432
433
434
435 CALL zlacgv( k-1, w( 1, kw ), 1 )
436 END IF
437
438 ELSE
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455 IF( k.GT.2 ) THEN
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499 d21 = w( k-1, kw )
500 d11 = w( k, kw ) / dconjg( d21 )
501 d22 = w( k-1, kw-1 ) / d21
502 t = one / ( dble( d11*d22 )-one )
503 d21 = t / d21
504
505
506
507
508
509 DO 20 j = 1, k - 2
510 a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
511 a( j, k ) = dconjg( d21 )*
512 $ ( d22*w( j, kw )-w( j, kw-1 ) )
513 20 CONTINUE
514 END IF
515
516
517
518 a( k-1, k-1 ) = w( k-1, kw-1 )
519 a( k-1, k ) = w( k-1, kw )
520 a( k, k ) = w( k, kw )
521
522
523
524 CALL zlacgv( k-1, w( 1, kw ), 1 )
525 CALL zlacgv( k-2, w( 1, kw-1 ), 1 )
526
527 END IF
528
529 END IF
530
531
532
533 IF( kstep.EQ.1 ) THEN
534 ipiv( k ) = kp
535 ELSE
536 ipiv( k ) = -kp
537 ipiv( k-1 ) = -kp
538 END IF
539
540
541
542 k = k - kstep
543 GO TO 10
544
545 30 CONTINUE
546
547
548
549
550
551
552
553
554 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
555 jb = min( nb, k-j+1 )
556
557
558
559 DO 40 jj = j, j + jb - 1
560 a( jj, jj ) = dble( a( jj, jj ) )
561 CALL zgemv(
'No transpose', jj-j+1, n-k, -cone,
562 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
563 $ a( j, jj ), 1 )
564 a( jj, jj ) = dble( a( jj, jj ) )
565 40 CONTINUE
566
567
568
569 CALL zgemm(
'No transpose',
'Transpose', j-1, jb, n-k,
570 $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
571 $ cone, a( 1, j ), lda )
572 50 CONTINUE
573
574
575
576
577 j = k + 1
578 60 CONTINUE
579
580
581
582
583
584 jj = j
585 jp = ipiv( j )
586 IF( jp.LT.0 ) THEN
587 jp = -jp
588
589 j = j + 1
590 END IF
591
592
593 j = j + 1
594 IF( jp.NE.jj .AND. j.LE.n )
595 $
CALL zswap( n-j+1, a( jp, j ), lda, a( jj, j ), lda )
596 IF( j.LT.n )
597 $ GO TO 60
598
599
600
601 kb = n - k
602
603 ELSE
604
605
606
607
608
609
610
611 k = 1
612 70 CONTINUE
613
614
615
616 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
617 $ GO TO 90
618
619 kstep = 1
620
621
622
623 w( k, k ) = dble( a( k, k ) )
624 IF( k.LT.n )
625 $
CALL zcopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
626 CALL zgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ), lda,
627 $ w( k, 1 ), ldw, cone, w( k, k ), 1 )
628 w( k, k ) = dble( w( k, k ) )
629
630
631
632
633 absakk = abs( dble( w( k, k ) ) )
634
635
636
637
638
639 IF( k.LT.n ) THEN
640 imax = k +
izamax( n-k, w( k+1, k ), 1 )
641 colmax = cabs1( w( imax, k ) )
642 ELSE
643 colmax = zero
644 END IF
645
646 IF( max( absakk, colmax ).EQ.zero ) THEN
647
648
649
650 IF( info.EQ.0 )
651 $ info = k
652 kp = k
653 a( k, k ) = dble( a( k, k ) )
654 ELSE
655
656
657
658
659
660
661 IF( absakk.GE.alpha*colmax ) THEN
662
663
664
665 kp = k
666 ELSE
667
668
669
670
671
672
673 CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
674 CALL zlacgv( imax-k, w( k, k+1 ), 1 )
675 w( imax, k+1 ) = dble( a( imax, imax ) )
676 IF( imax.LT.n )
677 $
CALL zcopy( n-imax, a( imax+1, imax ), 1,
678 $ w( imax+1, k+1 ), 1 )
679 CALL zgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
680 $ lda, w( imax, 1 ), ldw, cone, w( k, k+1 ),
681 $ 1 )
682 w( imax, k+1 ) = dble( w( imax, k+1 ) )
683
684
685
686
687
688 jmax = k - 1 +
izamax( imax-k, w( k, k+1 ), 1 )
689 rowmax = cabs1( w( jmax, k+1 ) )
690 IF( imax.LT.n ) THEN
691 jmax = imax +
izamax( n-imax, w( imax+1, k+1 ), 1 )
692 rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
693 END IF
694
695
696 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
697
698
699
700 kp = k
701
702
703 ELSE IF( abs( dble( w( imax, k+1 ) ) ).GE.alpha*rowmax )
704 $ THEN
705
706
707
708
709 kp = imax
710
711
712
713 CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
714
715
716 ELSE
717
718
719
720
721 kp = imax
722 kstep = 2
723 END IF
724
725
726
727
728 END IF
729
730
731
732
733
734
735
736 kk = k + kstep - 1
737
738
739
740
741 IF( kp.NE.kk ) THEN
742
743
744
745
746
747
748 a( kp, kp ) = dble( a( kk, kk ) )
749 CALL zcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
750 $ lda )
751 CALL zlacgv( kp-kk-1, a( kp, kk+1 ), lda )
752 IF( kp.LT.n )
753 $
CALL zcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
754
755
756
757
758
759
760 IF( k.GT.1 )
761 $
CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
762 CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
763 END IF
764
765 IF( kstep.EQ.1 ) THEN
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783 CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
784 IF( k.LT.n ) THEN
785
786
787
788
789
790 r1 = one / dble( a( k, k ) )
791 CALL zdscal( n-k, r1, a( k+1, k ), 1 )
792
793
794
795 CALL zlacgv( n-k, w( k+1, k ), 1 )
796 END IF
797
798 ELSE
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815 IF( k.LT.n-1 ) THEN
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859 d21 = w( k+1, k )
860 d11 = w( k+1, k+1 ) / d21
861 d22 = w( k, k ) / dconjg( d21 )
862 t = one / ( dble( d11*d22 )-one )
863 d21 = t / d21
864
865
866
867
868
869 DO 80 j = k + 2, n
870 a( j, k ) = dconjg( d21 )*
871 $ ( d11*w( j, k )-w( j, k+1 ) )
872 a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
873 80 CONTINUE
874 END IF
875
876
877
878 a( k, k ) = w( k, k )
879 a( k+1, k ) = w( k+1, k )
880 a( k+1, k+1 ) = w( k+1, k+1 )
881
882
883
884 CALL zlacgv( n-k, w( k+1, k ), 1 )
885 CALL zlacgv( n-k-1, w( k+2, k+1 ), 1 )
886
887 END IF
888
889 END IF
890
891
892
893 IF( kstep.EQ.1 ) THEN
894 ipiv( k ) = kp
895 ELSE
896 ipiv( k ) = -kp
897 ipiv( k+1 ) = -kp
898 END IF
899
900
901
902 k = k + kstep
903 GO TO 70
904
905 90 CONTINUE
906
907
908
909
910
911
912
913
914 DO 110 j = k, n, nb
915 jb = min( nb, n-j+1 )
916
917
918
919 DO 100 jj = j, j + jb - 1
920 a( jj, jj ) = dble( a( jj, jj ) )
921 CALL zgemv(
'No transpose', j+jb-jj, k-1, -cone,
922 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
923 $ a( jj, jj ), 1 )
924 a( jj, jj ) = dble( a( jj, jj ) )
925 100 CONTINUE
926
927
928
929 IF( j+jb.LE.n )
930 $
CALL zgemm(
'No transpose',
'Transpose', n-j-jb+1, jb,
931 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
932 $ ldw, cone, a( j+jb, j ), lda )
933 110 CONTINUE
934
935
936
937
938 j = k - 1
939 120 CONTINUE
940
941
942
943
944
945 jj = j
946 jp = ipiv( j )
947 IF( jp.LT.0 ) THEN
948 jp = -jp
949
950 j = j - 1
951 END IF
952
953
954 j = j - 1
955 IF( jp.NE.jj .AND. j.GE.1 )
956 $
CALL zswap( j, a( jp, 1 ), lda, a( jj, 1 ), lda )
957 IF( j.GT.1 )
958 $ GO TO 120
959
960
961
962 kb = k - 1
963
964 END IF
965 RETURN
966
967
968
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
integer function izamax(n, zx, incx)
IZAMAX
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
logical function lsame(ca, cb)
LSAME
subroutine zdscal(n, da, zx, incx)
ZDSCAL
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP