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