293
294
295
296
297
298
299 CHARACTER HOWMNY, SIDE
300 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
301
302
303 LOGICAL SELECT( * )
304 DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
305 $ VR( LDVR, * ), WORK( * )
306
307
308
309
310
311
312 DOUBLE PRECISION ZERO, ONE, SAFETY
313 parameter( zero = 0.0d+0, one = 1.0d+0,
314 $ safety = 1.0d+2 )
315
316
317 LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
318 $ ILBBAD, ILCOMP, ILCPLX, LSA, LSB
319 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
320 $ J, JA, JC, JE, JR, JW, NA, NW
321 DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
322 $ BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A,
323 $ CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA,
324 $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
325 $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
326 $ XSCALE
327
328
329 DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
330 $ SUMP( 2, 2 )
331
332
333 LOGICAL LSAME
334 DOUBLE PRECISION DLAMCH
336
337
340
341
342 INTRINSIC abs, max, min
343
344
345
346
347
348 IF(
lsame( howmny,
'A' ) )
THEN
349 ihwmny = 1
350 ilall = .true.
351 ilback = .false.
352 ELSE IF(
lsame( howmny,
'S' ) )
THEN
353 ihwmny = 2
354 ilall = .false.
355 ilback = .false.
356 ELSE IF(
lsame( howmny,
'B' ) )
THEN
357 ihwmny = 3
358 ilall = .true.
359 ilback = .true.
360 ELSE
361 ihwmny = -1
362 ilall = .true.
363 END IF
364
365 IF(
lsame( side,
'R' ) )
THEN
366 iside = 1
367 compl = .false.
368 compr = .true.
369 ELSE IF(
lsame( side,
'L' ) )
THEN
370 iside = 2
371 compl = .true.
372 compr = .false.
373 ELSE IF(
lsame( side,
'B' ) )
THEN
374 iside = 3
375 compl = .true.
376 compr = .true.
377 ELSE
378 iside = -1
379 END IF
380
381 info = 0
382 IF( iside.LT.0 ) THEN
383 info = -1
384 ELSE IF( ihwmny.LT.0 ) THEN
385 info = -2
386 ELSE IF( n.LT.0 ) THEN
387 info = -4
388 ELSE IF( lds.LT.max( 1, n ) ) THEN
389 info = -6
390 ELSE IF( ldp.LT.max( 1, n ) ) THEN
391 info = -8
392 END IF
393 IF( info.NE.0 ) THEN
394 CALL xerbla(
'DTGEVC', -info )
395 RETURN
396 END IF
397
398
399
400 IF( .NOT.ilall ) THEN
401 im = 0
402 ilcplx = .false.
403 DO 10 j = 1, n
404 IF( ilcplx ) THEN
405 ilcplx = .false.
406 GO TO 10
407 END IF
408 IF( j.LT.n ) THEN
409 IF( s( j+1, j ).NE.zero )
410 $ ilcplx = .true.
411 END IF
412 IF( ilcplx ) THEN
413 IF( SELECT( j ) .OR. SELECT( j+1 ) )
414 $ im = im + 2
415 ELSE
416 IF( SELECT( j ) )
417 $ im = im + 1
418 END IF
419 10 CONTINUE
420 ELSE
421 im = n
422 END IF
423
424
425
426 ilabad = .false.
427 ilbbad = .false.
428 DO 20 j = 1, n - 1
429 IF( s( j+1, j ).NE.zero ) THEN
430 IF( p( j, j ).EQ.zero .OR. p( j+1, j+1 ).EQ.zero .OR.
431 $ p( j, j+1 ).NE.zero )ilbbad = .true.
432 IF( j.LT.n-1 ) THEN
433 IF( s( j+2, j+1 ).NE.zero )
434 $ ilabad = .true.
435 END IF
436 END IF
437 20 CONTINUE
438
439 IF( ilabad ) THEN
440 info = -5
441 ELSE IF( ilbbad ) THEN
442 info = -7
443 ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
444 info = -10
445 ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
446 info = -12
447 ELSE IF( mm.LT.im ) THEN
448 info = -13
449 END IF
450 IF( info.NE.0 ) THEN
451 CALL xerbla(
'DTGEVC', -info )
452 RETURN
453 END IF
454
455
456
457 m = im
458 IF( n.EQ.0 )
459 $ RETURN
460
461
462
463 safmin =
dlamch(
'Safe minimum' )
464 big = one / safmin
466 small = safmin*n / ulp
467 big = one / small
468 bignum = one / ( safmin*n )
469
470
471
472
473
474
475 anorm = abs( s( 1, 1 ) )
476 IF( n.GT.1 )
477 $ anorm = anorm + abs( s( 2, 1 ) )
478 bnorm = abs( p( 1, 1 ) )
479 work( 1 ) = zero
480 work( n+1 ) = zero
481
482 DO 50 j = 2, n
483 temp = zero
484 temp2 = zero
485 IF( s( j, j-1 ).EQ.zero ) THEN
486 iend = j - 1
487 ELSE
488 iend = j - 2
489 END IF
490 DO 30 i = 1, iend
491 temp = temp + abs( s( i, j ) )
492 temp2 = temp2 + abs( p( i, j ) )
493 30 CONTINUE
494 work( j ) = temp
495 work( n+j ) = temp2
496 DO 40 i = iend + 1, min( j+1, n )
497 temp = temp + abs( s( i, j ) )
498 temp2 = temp2 + abs( p( i, j ) )
499 40 CONTINUE
500 anorm = max( anorm, temp )
501 bnorm = max( bnorm, temp2 )
502 50 CONTINUE
503
504 ascale = one / max( anorm, safmin )
505 bscale = one / max( bnorm, safmin )
506
507
508
509 IF( compl ) THEN
510 ieig = 0
511
512
513
514 ilcplx = .false.
515 DO 220 je = 1, n
516
517
518
519
520
521
522 IF( ilcplx ) THEN
523 ilcplx = .false.
524 GO TO 220
525 END IF
526 nw = 1
527 IF( je.LT.n ) THEN
528 IF( s( je+1, je ).NE.zero ) THEN
529 ilcplx = .true.
530 nw = 2
531 END IF
532 END IF
533 IF( ilall ) THEN
534 ilcomp = .true.
535 ELSE IF( ilcplx ) THEN
536 ilcomp = SELECT( je ) .OR. SELECT( je+1 )
537 ELSE
538 ilcomp = SELECT( je )
539 END IF
540 IF( .NOT.ilcomp )
541 $ GO TO 220
542
543
544
545
546 IF( .NOT.ilcplx ) THEN
547 IF( abs( s( je, je ) ).LE.safmin .AND.
548 $ abs( p( je, je ) ).LE.safmin ) THEN
549
550
551
552 ieig = ieig + 1
553 DO 60 jr = 1, n
554 vl( jr, ieig ) = zero
555 60 CONTINUE
556 vl( ieig, ieig ) = one
557 GO TO 220
558 END IF
559 END IF
560
561
562
563 DO 70 jr = 1, nw*n
564 work( 2*n+jr ) = zero
565 70 CONTINUE
566
567
568
569
570
571 IF( .NOT.ilcplx ) THEN
572
573
574
575 temp = one / max( abs( s( je, je ) )*ascale,
576 $ abs( p( je, je ) )*bscale, safmin )
577 salfar = ( temp*s( je, je ) )*ascale
578 sbeta = ( temp*p( je, je ) )*bscale
579 acoef = sbeta*ascale
580 bcoefr = salfar*bscale
581 bcoefi = zero
582
583
584
585 scale = one
586 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
587 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
588 $ small
589 IF( lsa )
590 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
591 IF( lsb )
592 $ scale = max( scale, ( small / abs( salfar ) )*
593 $ min( bnorm, big ) )
594 IF( lsa .OR. lsb ) THEN
595 scale = min( scale, one /
596 $ ( safmin*max( one, abs( acoef ),
597 $ abs( bcoefr ) ) ) )
598 IF( lsa ) THEN
599 acoef = ascale*( scale*sbeta )
600 ELSE
601 acoef = scale*acoef
602 END IF
603 IF( lsb ) THEN
604 bcoefr = bscale*( scale*salfar )
605 ELSE
606 bcoefr = scale*bcoefr
607 END IF
608 END IF
609 acoefa = abs( acoef )
610 bcoefa = abs( bcoefr )
611
612
613
614 work( 2*n+je ) = one
615 xmax = one
616 ELSE
617
618
619
620 CALL dlag2( s( je, je ), lds, p( je, je ), ldp,
621 $ safmin*safety, acoef, temp, bcoefr, temp2,
622 $ bcoefi )
623 bcoefi = -bcoefi
624 IF( bcoefi.EQ.zero ) THEN
625 info = je
626 RETURN
627 END IF
628
629
630
631 acoefa = abs( acoef )
632 bcoefa = abs( bcoefr ) + abs( bcoefi )
633 scale = one
634 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
635 $ scale = ( safmin / ulp ) / acoefa
636 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
637 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
638 IF( safmin*acoefa.GT.ascale )
639 $ scale = ascale / ( safmin*acoefa )
640 IF( safmin*bcoefa.GT.bscale )
641 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
642 IF( scale.NE.one ) THEN
643 acoef = scale*acoef
644 acoefa = abs( acoef )
645 bcoefr = scale*bcoefr
646 bcoefi = scale*bcoefi
647 bcoefa = abs( bcoefr ) + abs( bcoefi )
648 END IF
649
650
651
652 temp = acoef*s( je+1, je )
653 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
654 temp2i = -bcoefi*p( je, je )
655 IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) ) THEN
656 work( 2*n+je ) = one
657 work( 3*n+je ) = zero
658 work( 2*n+je+1 ) = -temp2r / temp
659 work( 3*n+je+1 ) = -temp2i / temp
660 ELSE
661 work( 2*n+je+1 ) = one
662 work( 3*n+je+1 ) = zero
663 temp = acoef*s( je, je+1 )
664 work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
665 $ s( je+1, je+1 ) ) / temp
666 work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
667 END IF
668 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
669 $ abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
670 END IF
671
672 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
673
674
675
676
677
678
679
680 il2by2 = .false.
681
682 DO 160 j = je + nw, n
683 IF( il2by2 ) THEN
684 il2by2 = .false.
685 GO TO 160
686 END IF
687
688 na = 1
689 bdiag( 1 ) = p( j, j )
690 IF( j.LT.n ) THEN
691 IF( s( j+1, j ).NE.zero ) THEN
692 il2by2 = .true.
693 bdiag( 2 ) = p( j+1, j+1 )
694 na = 2
695 END IF
696 END IF
697
698
699
700 xscale = one / max( one, xmax )
701 temp = max( work( j ), work( n+j ),
702 $ acoefa*work( j )+bcoefa*work( n+j ) )
703 IF( il2by2 )
704 $ temp = max( temp, work( j+1 ), work( n+j+1 ),
705 $ acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
706 IF( temp.GT.bignum*xscale ) THEN
707 DO 90 jw = 0, nw - 1
708 DO 80 jr = je, j - 1
709 work( ( jw+2 )*n+jr ) = xscale*
710 $ work( ( jw+2 )*n+jr )
711 80 CONTINUE
712 90 CONTINUE
713 xmax = xmax*xscale
714 END IF
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732 DO 120 jw = 1, nw
733 DO 110 ja = 1, na
734 sums( ja, jw ) = zero
735 sump( ja, jw ) = zero
736
737 DO 100 jr = je, j - 1
738 sums( ja, jw ) = sums( ja, jw ) +
739 $ s( jr, j+ja-1 )*
740 $ work( ( jw+1 )*n+jr )
741 sump( ja, jw ) = sump( ja, jw ) +
742 $ p( jr, j+ja-1 )*
743 $ work( ( jw+1 )*n+jr )
744 100 CONTINUE
745 110 CONTINUE
746 120 CONTINUE
747
748 DO 130 ja = 1, na
749 IF( ilcplx ) THEN
750 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
751 $ bcoefr*sump( ja, 1 ) -
752 $ bcoefi*sump( ja, 2 )
753 sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
754 $ bcoefr*sump( ja, 2 ) +
755 $ bcoefi*sump( ja, 1 )
756 ELSE
757 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
758 $ bcoefr*sump( ja, 1 )
759 END IF
760 130 CONTINUE
761
762
763
764
765
766 CALL dlaln2( .true., na, nw, dmin, acoef, s( j, j ),
767 $ lds,
768 $ bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
769 $ bcoefi, work( 2*n+j ), n, scale, temp,
770 $ iinfo )
771 IF( scale.LT.one ) THEN
772 DO 150 jw = 0, nw - 1
773 DO 140 jr = je, j - 1
774 work( ( jw+2 )*n+jr ) = scale*
775 $ work( ( jw+2 )*n+jr )
776 140 CONTINUE
777 150 CONTINUE
778 xmax = scale*xmax
779 END IF
780 xmax = max( xmax, temp )
781 160 CONTINUE
782
783
784
785
786 ieig = ieig + 1
787 IF( ilback ) THEN
788 DO 170 jw = 0, nw - 1
789 CALL dgemv(
'N', n, n+1-je, one, vl( 1, je ), ldvl,
790 $ work( ( jw+2 )*n+je ), 1, zero,
791 $ work( ( jw+4 )*n+1 ), 1 )
792 170 CONTINUE
793 CALL dlacpy(
' ', n, nw, work( 4*n+1 ), n, vl( 1,
794 $ je ),
795 $ ldvl )
796 ibeg = 1
797 ELSE
798 CALL dlacpy(
' ', n, nw, work( 2*n+1 ), n, vl( 1,
799 $ ieig ),
800 $ ldvl )
801 ibeg = je
802 END IF
803
804
805
806 xmax = zero
807 IF( ilcplx ) THEN
808 DO 180 j = ibeg, n
809 xmax = max( xmax, abs( vl( j, ieig ) )+
810 $ abs( vl( j, ieig+1 ) ) )
811 180 CONTINUE
812 ELSE
813 DO 190 j = ibeg, n
814 xmax = max( xmax, abs( vl( j, ieig ) ) )
815 190 CONTINUE
816 END IF
817
818 IF( xmax.GT.safmin ) THEN
819 xscale = one / xmax
820
821 DO 210 jw = 0, nw - 1
822 DO 200 jr = ibeg, n
823 vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
824 200 CONTINUE
825 210 CONTINUE
826 END IF
827 ieig = ieig + nw - 1
828
829 220 CONTINUE
830 END IF
831
832
833
834 IF( compr ) THEN
835 ieig = im + 1
836
837
838
839 ilcplx = .false.
840 DO 500 je = n, 1, -1
841
842
843
844
845
846
847
848
849
850 IF( ilcplx ) THEN
851 ilcplx = .false.
852 GO TO 500
853 END IF
854 nw = 1
855 IF( je.GT.1 ) THEN
856 IF( s( je, je-1 ).NE.zero ) THEN
857 ilcplx = .true.
858 nw = 2
859 END IF
860 END IF
861 IF( ilall ) THEN
862 ilcomp = .true.
863 ELSE IF( ilcplx ) THEN
864 ilcomp = SELECT( je ) .OR. SELECT( je-1 )
865 ELSE
866 ilcomp = SELECT( je )
867 END IF
868 IF( .NOT.ilcomp )
869 $ GO TO 500
870
871
872
873
874 IF( .NOT.ilcplx ) THEN
875 IF( abs( s( je, je ) ).LE.safmin .AND.
876 $ abs( p( je, je ) ).LE.safmin ) THEN
877
878
879
880 ieig = ieig - 1
881 DO 230 jr = 1, n
882 vr( jr, ieig ) = zero
883 230 CONTINUE
884 vr( ieig, ieig ) = one
885 GO TO 500
886 END IF
887 END IF
888
889
890
891 DO 250 jw = 0, nw - 1
892 DO 240 jr = 1, n
893 work( ( jw+2 )*n+jr ) = zero
894 240 CONTINUE
895 250 CONTINUE
896
897
898
899
900
901 IF( .NOT.ilcplx ) THEN
902
903
904
905 temp = one / max( abs( s( je, je ) )*ascale,
906 $ abs( p( je, je ) )*bscale, safmin )
907 salfar = ( temp*s( je, je ) )*ascale
908 sbeta = ( temp*p( je, je ) )*bscale
909 acoef = sbeta*ascale
910 bcoefr = salfar*bscale
911 bcoefi = zero
912
913
914
915 scale = one
916 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
917 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
918 $ small
919 IF( lsa )
920 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
921 IF( lsb )
922 $ scale = max( scale, ( small / abs( salfar ) )*
923 $ min( bnorm, big ) )
924 IF( lsa .OR. lsb ) THEN
925 scale = min( scale, one /
926 $ ( safmin*max( one, abs( acoef ),
927 $ abs( bcoefr ) ) ) )
928 IF( lsa ) THEN
929 acoef = ascale*( scale*sbeta )
930 ELSE
931 acoef = scale*acoef
932 END IF
933 IF( lsb ) THEN
934 bcoefr = bscale*( scale*salfar )
935 ELSE
936 bcoefr = scale*bcoefr
937 END IF
938 END IF
939 acoefa = abs( acoef )
940 bcoefa = abs( bcoefr )
941
942
943
944 work( 2*n+je ) = one
945 xmax = one
946
947
948
949
950 DO 260 jr = 1, je - 1
951 work( 2*n+jr ) = bcoefr*p( jr, je ) -
952 $ acoef*s( jr, je )
953 260 CONTINUE
954 ELSE
955
956
957
958 CALL dlag2( s( je-1, je-1 ), lds, p( je-1, je-1 ),
959 $ ldp,
960 $ safmin*safety, acoef, temp, bcoefr, temp2,
961 $ bcoefi )
962 IF( bcoefi.EQ.zero ) THEN
963 info = je - 1
964 RETURN
965 END IF
966
967
968
969 acoefa = abs( acoef )
970 bcoefa = abs( bcoefr ) + abs( bcoefi )
971 scale = one
972 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
973 $ scale = ( safmin / ulp ) / acoefa
974 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
975 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
976 IF( safmin*acoefa.GT.ascale )
977 $ scale = ascale / ( safmin*acoefa )
978 IF( safmin*bcoefa.GT.bscale )
979 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
980 IF( scale.NE.one ) THEN
981 acoef = scale*acoef
982 acoefa = abs( acoef )
983 bcoefr = scale*bcoefr
984 bcoefi = scale*bcoefi
985 bcoefa = abs( bcoefr ) + abs( bcoefi )
986 END IF
987
988
989
990
991 temp = acoef*s( je, je-1 )
992 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
993 temp2i = -bcoefi*p( je, je )
994 IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) ) THEN
995 work( 2*n+je ) = one
996 work( 3*n+je ) = zero
997 work( 2*n+je-1 ) = -temp2r / temp
998 work( 3*n+je-1 ) = -temp2i / temp
999 ELSE
1000 work( 2*n+je-1 ) = one
1001 work( 3*n+je-1 ) = zero
1002 temp = acoef*s( je-1, je )
1003 work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
1004 $ s( je-1, je-1 ) ) / temp
1005 work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
1006 END IF
1007
1008 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
1009 $ abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
1010
1011
1012
1013
1014 creala = acoef*work( 2*n+je-1 )
1015 cimaga = acoef*work( 3*n+je-1 )
1016 crealb = bcoefr*work( 2*n+je-1 ) -
1017 $ bcoefi*work( 3*n+je-1 )
1018 cimagb = bcoefi*work( 2*n+je-1 ) +
1019 $ bcoefr*work( 3*n+je-1 )
1020 cre2a = acoef*work( 2*n+je )
1021 cim2a = acoef*work( 3*n+je )
1022 cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
1023 cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
1024 DO 270 jr = 1, je - 2
1025 work( 2*n+jr ) = -creala*s( jr, je-1 ) +
1026 $ crealb*p( jr, je-1 ) -
1027 $ cre2a*s( jr, je ) + cre2b*p( jr, je )
1028 work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
1029 $ cimagb*p( jr, je-1 ) -
1030 $ cim2a*s( jr, je ) + cim2b*p( jr, je )
1031 270 CONTINUE
1032 END IF
1033
1034 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
1035
1036
1037
1038 il2by2 = .false.
1039 DO 370 j = je - nw, 1, -1
1040
1041
1042
1043
1044 IF( .NOT.il2by2 .AND. j.GT.1 ) THEN
1045 IF( s( j, j-1 ).NE.zero ) THEN
1046 il2by2 = .true.
1047 GO TO 370
1048 END IF
1049 END IF
1050 bdiag( 1 ) = p( j, j )
1051 IF( il2by2 ) THEN
1052 na = 2
1053 bdiag( 2 ) = p( j+1, j+1 )
1054 ELSE
1055 na = 1
1056 END IF
1057
1058
1059
1060 CALL dlaln2( .false., na, nw, dmin, acoef, s( j, j ),
1061 $ lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
1062 $ n, bcoefr, bcoefi, sum, 2, scale, temp,
1063 $ iinfo )
1064 IF( scale.LT.one ) THEN
1065
1066 DO 290 jw = 0, nw - 1
1067 DO 280 jr = 1, je
1068 work( ( jw+2 )*n+jr ) = scale*
1069 $ work( ( jw+2 )*n+jr )
1070 280 CONTINUE
1071 290 CONTINUE
1072 END IF
1073 xmax = max( scale*xmax, temp )
1074
1075 DO 310 jw = 1, nw
1076 DO 300 ja = 1, na
1077 work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
1078 300 CONTINUE
1079 310 CONTINUE
1080
1081
1082
1083 IF( j.GT.1 ) THEN
1084
1085
1086
1087 xscale = one / max( one, xmax )
1088 temp = acoefa*work( j ) + bcoefa*work( n+j )
1089 IF( il2by2 )
1090 $ temp = max( temp, acoefa*work( j+1 )+bcoefa*
1091 $ work( n+j+1 ) )
1092 temp = max( temp, acoefa, bcoefa )
1093 IF( temp.GT.bignum*xscale ) THEN
1094
1095 DO 330 jw = 0, nw - 1
1096 DO 320 jr = 1, je
1097 work( ( jw+2 )*n+jr ) = xscale*
1098 $ work( ( jw+2 )*n+jr )
1099 320 CONTINUE
1100 330 CONTINUE
1101 xmax = xmax*xscale
1102 END IF
1103
1104
1105
1106
1107
1108
1109 DO 360 ja = 1, na
1110 IF( ilcplx ) THEN
1111 creala = acoef*work( 2*n+j+ja-1 )
1112 cimaga = acoef*work( 3*n+j+ja-1 )
1113 crealb = bcoefr*work( 2*n+j+ja-1 ) -
1114 $ bcoefi*work( 3*n+j+ja-1 )
1115 cimagb = bcoefi*work( 2*n+j+ja-1 ) +
1116 $ bcoefr*work( 3*n+j+ja-1 )
1117 DO 340 jr = 1, j - 1
1118 work( 2*n+jr ) = work( 2*n+jr ) -
1119 $ creala*s( jr, j+ja-1 ) +
1120 $ crealb*p( jr, j+ja-1 )
1121 work( 3*n+jr ) = work( 3*n+jr ) -
1122 $ cimaga*s( jr, j+ja-1 ) +
1123 $ cimagb*p( jr, j+ja-1 )
1124 340 CONTINUE
1125 ELSE
1126 creala = acoef*work( 2*n+j+ja-1 )
1127 crealb = bcoefr*work( 2*n+j+ja-1 )
1128 DO 350 jr = 1, j - 1
1129 work( 2*n+jr ) = work( 2*n+jr ) -
1130 $ creala*s( jr, j+ja-1 ) +
1131 $ crealb*p( jr, j+ja-1 )
1132 350 CONTINUE
1133 END IF
1134 360 CONTINUE
1135 END IF
1136
1137 il2by2 = .false.
1138 370 CONTINUE
1139
1140
1141
1142
1143 ieig = ieig - nw
1144 IF( ilback ) THEN
1145
1146 DO 410 jw = 0, nw - 1
1147 DO 380 jr = 1, n
1148 work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
1149 $ vr( jr, 1 )
1150 380 CONTINUE
1151
1152
1153
1154
1155
1156 DO 400 jc = 2, je
1157 DO 390 jr = 1, n
1158 work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
1159 $ work( ( jw+2 )*n+jc )*vr( jr, jc )
1160 390 CONTINUE
1161 400 CONTINUE
1162 410 CONTINUE
1163
1164 DO 430 jw = 0, nw - 1
1165 DO 420 jr = 1, n
1166 vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
1167 420 CONTINUE
1168 430 CONTINUE
1169
1170 iend = n
1171 ELSE
1172 DO 450 jw = 0, nw - 1
1173 DO 440 jr = 1, n
1174 vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
1175 440 CONTINUE
1176 450 CONTINUE
1177
1178 iend = je
1179 END IF
1180
1181
1182
1183 xmax = zero
1184 IF( ilcplx ) THEN
1185 DO 460 j = 1, iend
1186 xmax = max( xmax, abs( vr( j, ieig ) )+
1187 $ abs( vr( j, ieig+1 ) ) )
1188 460 CONTINUE
1189 ELSE
1190 DO 470 j = 1, iend
1191 xmax = max( xmax, abs( vr( j, ieig ) ) )
1192 470 CONTINUE
1193 END IF
1194
1195 IF( xmax.GT.safmin ) THEN
1196 xscale = one / xmax
1197 DO 490 jw = 0, nw - 1
1198 DO 480 jr = 1, iend
1199 vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
1200 480 CONTINUE
1201 490 CONTINUE
1202 END IF
1203 500 CONTINUE
1204 END IF
1205
1206 RETURN
1207
1208
1209
subroutine xerbla(srname, info)
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
subroutine dlaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
double precision function dlamch(cmach)
DLAMCH
logical function lsame(ca, cb)
LSAME