239
240
241
242
243
244
245 CHARACTER DIAG, NORMIN, TRANS, UPLO
246 INTEGER INFO, LDA, N
247 REAL SCALE
248
249
250 REAL CNORM( * )
251 COMPLEX A( LDA, * ), X( * )
252
253
254
255
256
257 REAL ZERO, HALF, ONE, TWO
258 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
259 $ two = 2.0e+0 )
260
261
262 LOGICAL NOTRAN, NOUNIT, UPPER
263 INTEGER I, IMAX, J, JFIRST, JINC, JLAST
264 REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
265 $ XBND, XJ, XMAX
266 COMPLEX CSUMJ, TJJS, USCAL, ZDUM
267
268
269 LOGICAL LSAME
270 INTEGER ICAMAX, ISAMAX
271 REAL SCASUM, SLAMCH
272 COMPLEX CDOTC, CDOTU, CLADIV
275
276
278
279
280 INTRINSIC abs, aimag, cmplx, conjg, max, min, real
281
282
283 REAL CABS1, CABS2
284
285
286 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
287 cabs2( zdum ) = abs( real( zdum ) / 2. ) +
288 $ abs( aimag( zdum ) / 2. )
289
290
291
292 info = 0
293 upper =
lsame( uplo,
'U' )
294 notran =
lsame( trans,
'N' )
295 nounit =
lsame( diag,
'N' )
296
297
298
299 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
300 info = -1
301 ELSE IF( .NOT.notran .AND. .NOT.
lsame( trans,
'T' ) .AND. .NOT.
302 $
lsame( trans,
'C' ) )
THEN
303 info = -2
304 ELSE IF( .NOT.nounit .AND. .NOT.
lsame( diag,
'U' ) )
THEN
305 info = -3
306 ELSE IF( .NOT.
lsame( normin,
'Y' ) .AND. .NOT.
307 $
lsame( normin,
'N' ) )
THEN
308 info = -4
309 ELSE IF( n.LT.0 ) THEN
310 info = -5
311 ELSE IF( lda.LT.max( 1, n ) ) THEN
312 info = -7
313 END IF
314 IF( info.NE.0 ) THEN
315 CALL xerbla(
'CLATRS', -info )
316 RETURN
317 END IF
318
319
320
321 scale = one
322 IF( n.EQ.0 )
323 $ RETURN
324
325
326
327 smlnum =
slamch(
'Safe minimum' ) /
slamch(
'Precision' )
328 bignum = one / smlnum
329
330 IF(
lsame( normin,
'N' ) )
THEN
331
332
333
334 IF( upper ) THEN
335
336
337
338 DO 10 j = 1, n
339 cnorm( j ) =
scasum( j-1, a( 1, j ), 1 )
340 10 CONTINUE
341 ELSE
342
343
344
345 DO 20 j = 1, n - 1
346 cnorm( j ) =
scasum( n-j, a( j+1, j ), 1 )
347 20 CONTINUE
348 cnorm( n ) = zero
349 END IF
350 END IF
351
352
353
354
355 imax =
isamax( n, cnorm, 1 )
356 tmax = cnorm( imax )
357 IF( tmax.LE.bignum*half ) THEN
358 tscal = one
359 ELSE
360
361
362
363
364 IF ( tmax.LE.
slamch(
'Overflow') )
THEN
365
366 tscal = half / ( smlnum*tmax )
367 CALL sscal( n, tscal, cnorm, 1 )
368 ELSE
369
370
371
372
373
374 tmax = zero
375 IF( upper ) THEN
376
377
378
379 DO j = 2, n
380 DO i = 1, j - 1
381 tmax = max( tmax, abs( real( a( i, j ) ) ),
382 $ abs( aimag(a( i, j ) ) ) )
383 END DO
384 END DO
385 ELSE
386
387
388
389 DO j = 1, n - 1
390 DO i = j + 1, n
391 tmax = max( tmax, abs( real( a( i, j ) ) ),
392 $ abs( aimag(a( i, j ) ) ) )
393 END DO
394 END DO
395 END IF
396
397 IF( tmax.LE.
slamch(
'Overflow') )
THEN
398 tscal = one / ( smlnum*tmax )
399 DO j = 1, n
400 IF( cnorm( j ).LE.
slamch(
'Overflow') )
THEN
401 cnorm( j ) = cnorm( j )*tscal
402 ELSE
403
404
405 tscal = two * tscal
406 cnorm( j ) = zero
407 IF( upper ) THEN
408 DO i = 1, j - 1
409 cnorm( j ) = cnorm( j ) +
410 $ tscal * cabs2( a( i, j ) )
411 END DO
412 ELSE
413 DO i = j + 1, n
414 cnorm( j ) = cnorm( j ) +
415 $ tscal * cabs2( a( i, j ) )
416 END DO
417 END IF
418 tscal = tscal * half
419 END IF
420 END DO
421 ELSE
422
423
424 CALL ctrsv( uplo, trans, diag, n, a, lda, x, 1 )
425 RETURN
426 END IF
427 END IF
428 END IF
429
430
431
432
433 xmax = zero
434 DO 30 j = 1, n
435 xmax = max( xmax, cabs2( x( j ) ) )
436 30 CONTINUE
437 xbnd = xmax
438
439 IF( notran ) THEN
440
441
442
443 IF( upper ) THEN
444 jfirst = n
445 jlast = 1
446 jinc = -1
447 ELSE
448 jfirst = 1
449 jlast = n
450 jinc = 1
451 END IF
452
453 IF( tscal.NE.one ) THEN
454 grow = zero
455 GO TO 60
456 END IF
457
458 IF( nounit ) THEN
459
460
461
462
463
464
465 grow = half / max( xbnd, smlnum )
466 xbnd = grow
467 DO 40 j = jfirst, jlast, jinc
468
469
470
471 IF( grow.LE.smlnum )
472 $ GO TO 60
473
474 tjjs = a( j, j )
475 tjj = cabs1( tjjs )
476
477 IF( tjj.GE.smlnum ) THEN
478
479
480
481 xbnd = min( xbnd, min( one, tjj )*grow )
482 ELSE
483
484
485
486 xbnd = zero
487 END IF
488
489 IF( tjj+cnorm( j ).GE.smlnum ) THEN
490
491
492
493 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
494 ELSE
495
496
497
498 grow = zero
499 END IF
500 40 CONTINUE
501 grow = xbnd
502 ELSE
503
504
505
506
507
508 grow = min( one, half / max( xbnd, smlnum ) )
509 DO 50 j = jfirst, jlast, jinc
510
511
512
513 IF( grow.LE.smlnum )
514 $ GO TO 60
515
516
517
518 grow = grow*( one / ( one+cnorm( j ) ) )
519 50 CONTINUE
520 END IF
521 60 CONTINUE
522
523 ELSE
524
525
526
527 IF( upper ) THEN
528 jfirst = 1
529 jlast = n
530 jinc = 1
531 ELSE
532 jfirst = n
533 jlast = 1
534 jinc = -1
535 END IF
536
537 IF( tscal.NE.one ) THEN
538 grow = zero
539 GO TO 90
540 END IF
541
542 IF( nounit ) THEN
543
544
545
546
547
548
549 grow = half / max( xbnd, smlnum )
550 xbnd = grow
551 DO 70 j = jfirst, jlast, jinc
552
553
554
555 IF( grow.LE.smlnum )
556 $ GO TO 90
557
558
559
560 xj = one + cnorm( j )
561 grow = min( grow, xbnd / xj )
562
563 tjjs = a( j, j )
564 tjj = cabs1( tjjs )
565
566 IF( tjj.GE.smlnum ) THEN
567
568
569
570 IF( xj.GT.tjj )
571 $ xbnd = xbnd*( tjj / xj )
572 ELSE
573
574
575
576 xbnd = zero
577 END IF
578 70 CONTINUE
579 grow = min( grow, xbnd )
580 ELSE
581
582
583
584
585
586 grow = min( one, half / max( xbnd, smlnum ) )
587 DO 80 j = jfirst, jlast, jinc
588
589
590
591 IF( grow.LE.smlnum )
592 $ GO TO 90
593
594
595
596 xj = one + cnorm( j )
597 grow = grow / xj
598 80 CONTINUE
599 END IF
600 90 CONTINUE
601 END IF
602
603 IF( ( grow*tscal ).GT.smlnum ) THEN
604
605
606
607
608 CALL ctrsv( uplo, trans, diag, n, a, lda, x, 1 )
609 ELSE
610
611
612
613 IF( xmax.GT.bignum*half ) THEN
614
615
616
617
618 scale = ( bignum*half ) / xmax
619 CALL csscal( n, scale, x, 1 )
620 xmax = bignum
621 ELSE
622 xmax = xmax*two
623 END IF
624
625 IF( notran ) THEN
626
627
628
629 DO 110 j = jfirst, jlast, jinc
630
631
632
633 xj = cabs1( x( j ) )
634 IF( nounit ) THEN
635 tjjs = a( j, j )*tscal
636 ELSE
637 tjjs = tscal
638 IF( tscal.EQ.one )
639 $ GO TO 105
640 END IF
641 tjj = cabs1( tjjs )
642 IF( tjj.GT.smlnum ) THEN
643
644
645
646 IF( tjj.LT.one ) THEN
647 IF( xj.GT.tjj*bignum ) THEN
648
649
650
651 rec = one / xj
652 CALL csscal( n, rec, x, 1 )
653 scale = scale*rec
654 xmax = xmax*rec
655 END IF
656 END IF
657 x( j ) =
cladiv( x( j ), tjjs )
658 xj = cabs1( x( j ) )
659 ELSE IF( tjj.GT.zero ) THEN
660
661
662
663 IF( xj.GT.tjj*bignum ) THEN
664
665
666
667
668 rec = ( tjj*bignum ) / xj
669 IF( cnorm( j ).GT.one ) THEN
670
671
672
673
674 rec = rec / cnorm( j )
675 END IF
676 CALL csscal( n, rec, x, 1 )
677 scale = scale*rec
678 xmax = xmax*rec
679 END IF
680 x( j ) =
cladiv( x( j ), tjjs )
681 xj = cabs1( x( j ) )
682 ELSE
683
684
685
686
687 DO 100 i = 1, n
688 x( i ) = zero
689 100 CONTINUE
690 x( j ) = one
691 xj = one
692 scale = zero
693 xmax = zero
694 END IF
695 105 CONTINUE
696
697
698
699
700 IF( xj.GT.one ) THEN
701 rec = one / xj
702 IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
703
704
705
706 rec = rec*half
707 CALL csscal( n, rec, x, 1 )
708 scale = scale*rec
709 END IF
710 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
711
712
713
714 CALL csscal( n, half, x, 1 )
715 scale = scale*half
716 END IF
717
718 IF( upper ) THEN
719 IF( j.GT.1 ) THEN
720
721
722
723
724 CALL caxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
725 $ 1 )
727 xmax = cabs1( x( i ) )
728 END IF
729 ELSE
730 IF( j.LT.n ) THEN
731
732
733
734
735 CALL caxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
736 $ x( j+1 ), 1 )
737 i = j +
icamax( n-j, x( j+1 ), 1 )
738 xmax = cabs1( x( i ) )
739 END IF
740 END IF
741 110 CONTINUE
742
743 ELSE IF(
lsame( trans,
'T' ) )
THEN
744
745
746
747 DO 150 j = jfirst, jlast, jinc
748
749
750
751
752 xj = cabs1( x( j ) )
753 uscal = tscal
754 rec = one / max( xmax, one )
755 IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
756
757
758
759 rec = rec*half
760 IF( nounit ) THEN
761 tjjs = a( j, j )*tscal
762 ELSE
763 tjjs = tscal
764 END IF
765 tjj = cabs1( tjjs )
766 IF( tjj.GT.one ) THEN
767
768
769
770 rec = min( one, rec*tjj )
771 uscal =
cladiv( uscal, tjjs )
772 END IF
773 IF( rec.LT.one ) THEN
774 CALL csscal( n, rec, x, 1 )
775 scale = scale*rec
776 xmax = xmax*rec
777 END IF
778 END IF
779
780 csumj = zero
781 IF( uscal.EQ.cmplx( one ) ) THEN
782
783
784
785
786 IF( upper ) THEN
787 csumj =
cdotu( j-1, a( 1, j ), 1, x, 1 )
788 ELSE IF( j.LT.n ) THEN
789 csumj =
cdotu( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
790 END IF
791 ELSE
792
793
794
795 IF( upper ) THEN
796 DO 120 i = 1, j - 1
797 csumj = csumj + ( a( i, j )*uscal )*x( i )
798 120 CONTINUE
799 ELSE IF( j.LT.n ) THEN
800 DO 130 i = j + 1, n
801 csumj = csumj + ( a( i, j )*uscal )*x( i )
802 130 CONTINUE
803 END IF
804 END IF
805
806 IF( uscal.EQ.cmplx( tscal ) ) THEN
807
808
809
810
811 x( j ) = x( j ) - csumj
812 xj = cabs1( x( j ) )
813 IF( nounit ) THEN
814 tjjs = a( j, j )*tscal
815 ELSE
816 tjjs = tscal
817 IF( tscal.EQ.one )
818 $ GO TO 145
819 END IF
820
821
822
823 tjj = cabs1( tjjs )
824 IF( tjj.GT.smlnum ) THEN
825
826
827
828 IF( tjj.LT.one ) THEN
829 IF( xj.GT.tjj*bignum ) THEN
830
831
832
833 rec = one / xj
834 CALL csscal( n, rec, x, 1 )
835 scale = scale*rec
836 xmax = xmax*rec
837 END IF
838 END IF
839 x( j ) =
cladiv( x( j ), tjjs )
840 ELSE IF( tjj.GT.zero ) THEN
841
842
843
844 IF( xj.GT.tjj*bignum ) THEN
845
846
847
848 rec = ( tjj*bignum ) / xj
849 CALL csscal( n, rec, x, 1 )
850 scale = scale*rec
851 xmax = xmax*rec
852 END IF
853 x( j ) =
cladiv( x( j ), tjjs )
854 ELSE
855
856
857
858
859 DO 140 i = 1, n
860 x( i ) = zero
861 140 CONTINUE
862 x( j ) = one
863 scale = zero
864 xmax = zero
865 END IF
866 145 CONTINUE
867 ELSE
868
869
870
871
872 x( j ) =
cladiv( x( j ), tjjs ) - csumj
873 END IF
874 xmax = max( xmax, cabs1( x( j ) ) )
875 150 CONTINUE
876
877 ELSE
878
879
880
881 DO 190 j = jfirst, jlast, jinc
882
883
884
885
886 xj = cabs1( x( j ) )
887 uscal = tscal
888 rec = one / max( xmax, one )
889 IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
890
891
892
893 rec = rec*half
894 IF( nounit ) THEN
895 tjjs = conjg( a( j, j ) )*tscal
896 ELSE
897 tjjs = tscal
898 END IF
899 tjj = cabs1( tjjs )
900 IF( tjj.GT.one ) THEN
901
902
903
904 rec = min( one, rec*tjj )
905 uscal =
cladiv( uscal, tjjs )
906 END IF
907 IF( rec.LT.one ) THEN
908 CALL csscal( n, rec, x, 1 )
909 scale = scale*rec
910 xmax = xmax*rec
911 END IF
912 END IF
913
914 csumj = zero
915 IF( uscal.EQ.cmplx( one ) ) THEN
916
917
918
919
920 IF( upper ) THEN
921 csumj =
cdotc( j-1, a( 1, j ), 1, x, 1 )
922 ELSE IF( j.LT.n ) THEN
923 csumj =
cdotc( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
924 END IF
925 ELSE
926
927
928
929 IF( upper ) THEN
930 DO 160 i = 1, j - 1
931 csumj = csumj + ( conjg( a( i, j ) )*uscal )*
932 $ x( i )
933 160 CONTINUE
934 ELSE IF( j.LT.n ) THEN
935 DO 170 i = j + 1, n
936 csumj = csumj + ( conjg( a( i, j ) )*uscal )*
937 $ x( i )
938 170 CONTINUE
939 END IF
940 END IF
941
942 IF( uscal.EQ.cmplx( tscal ) ) THEN
943
944
945
946
947 x( j ) = x( j ) - csumj
948 xj = cabs1( x( j ) )
949 IF( nounit ) THEN
950 tjjs = conjg( a( j, j ) )*tscal
951 ELSE
952 tjjs = tscal
953 IF( tscal.EQ.one )
954 $ GO TO 185
955 END IF
956
957
958
959 tjj = cabs1( tjjs )
960 IF( tjj.GT.smlnum ) THEN
961
962
963
964 IF( tjj.LT.one ) THEN
965 IF( xj.GT.tjj*bignum ) THEN
966
967
968
969 rec = one / xj
970 CALL csscal( n, rec, x, 1 )
971 scale = scale*rec
972 xmax = xmax*rec
973 END IF
974 END IF
975 x( j ) =
cladiv( x( j ), tjjs )
976 ELSE IF( tjj.GT.zero ) THEN
977
978
979
980 IF( xj.GT.tjj*bignum ) THEN
981
982
983
984 rec = ( tjj*bignum ) / xj
985 CALL csscal( n, rec, x, 1 )
986 scale = scale*rec
987 xmax = xmax*rec
988 END IF
989 x( j ) =
cladiv( x( j ), tjjs )
990 ELSE
991
992
993
994
995 DO 180 i = 1, n
996 x( i ) = zero
997 180 CONTINUE
998 x( j ) = one
999 scale = zero
1000 xmax = zero
1001 END IF
1002 185 CONTINUE
1003 ELSE
1004
1005
1006
1007
1008 x( j ) =
cladiv( x( j ), tjjs ) - csumj
1009 END IF
1010 xmax = max( xmax, cabs1( x( j ) ) )
1011 190 CONTINUE
1012 END IF
1013 scale = scale / tscal
1014 END IF
1015
1016
1017
1018 IF( tscal.NE.one ) THEN
1019 CALL sscal( n, one / tscal, cnorm, 1 )
1020 END IF
1021
1022 RETURN
1023
1024
1025
subroutine xerbla(srname, info)
real function scasum(n, cx, incx)
SCASUM
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
complex function cdotu(n, cx, incx, cy, incy)
CDOTU
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
integer function isamax(n, sx, incx)
ISAMAX
integer function icamax(n, cx, incx)
ICAMAX
complex function cladiv(x, y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
real function slamch(cmach)
SLAMCH
logical function lsame(ca, cb)
LSAME
subroutine csscal(n, sa, cx, incx)
CSSCAL
subroutine sscal(n, sa, sx, incx)
SSCAL
subroutine ctrsv(uplo, trans, diag, n, a, lda, x, incx)
CTRSV