343
344
345
346
347
348
349 LOGICAL COMP
350 INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
351 DOUBLE PRECISION RCDEIN, RCDVIN, THRESH
352
353
354 LOGICAL BWORK( * )
355 INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * )
356 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), HT( LDA, * ),
357 $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
358 $ WI( * ), WIT( * ), WITMP( * ), WORK( * ),
359 $ WR( * ), WRT( * ), WRTMP( * )
360
361
362
363
364
365 DOUBLE PRECISION ZERO, ONE
366 parameter( zero = 0.0d0, one = 1.0d0 )
367 DOUBLE PRECISION EPSIN
368 parameter( epsin = 5.9605d-8 )
369
370
371 CHARACTER SORT
372 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
373 $ RSUB, SDIM, SDIM1
374 DOUBLE PRECISION ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
375 $ SMLNUM, TMP, TOL, TOLIN, ULP, ULPINV, V, VIMIN,
376 $ VRMIN, WNORM
377
378
379 INTEGER IPNT( 20 )
380
381
382 LOGICAL SELVAL( 20 )
383 DOUBLE PRECISION SELWI( 20 ), SELWR( 20 )
384
385
386 INTEGER SELDIM, SELOPT
387
388
389 COMMON / sslct / selopt, seldim, selval, selwr, selwi
390
391
392 LOGICAL DSLECT
393 DOUBLE PRECISION DLAMCH, DLANGE
395
396
398
399
400 INTRINSIC abs, dble, max, min, sign, sqrt
401
402
403
404
405
406 info = 0
407 IF( thresh.LT.zero ) THEN
408 info = -3
409 ELSE IF( nounit.LE.0 ) THEN
410 info = -5
411 ELSE IF( n.LT.0 ) THEN
412 info = -6
413 ELSE IF( lda.LT.1 .OR. lda.LT.n ) THEN
414 info = -8
415 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n ) THEN
416 info = -18
417 ELSE IF( lwork.LT.3*n ) THEN
418 info = -26
419 END IF
420
421 IF( info.NE.0 ) THEN
422 CALL xerbla(
'DGET24', -info )
423 RETURN
424 END IF
425
426
427
428 DO 10 i = 1, 17
429 result( i ) = -one
430 10 CONTINUE
431
432 IF( n.EQ.0 )
433 $ RETURN
434
435
436
437 smlnum =
dlamch(
'Safe minimum' )
438 ulp =
dlamch(
'Precision' )
439 ulpinv = one / ulp
440
441
442
443 selopt = 0
444 liwork = n*n
445 DO 120 isort = 0, 1
446 IF( isort.EQ.0 ) THEN
447 sort = 'N'
448 rsub = 0
449 ELSE
450 sort = 'S'
451 rsub = 6
452 END IF
453
454
455
456 CALL dlacpy(
'F', n, n, a, lda, h, lda )
457 CALL dgeesx(
'V', sort,
dslect,
'N', n, h, lda, sdim, wr, wi,
458 $ vs, ldvs, rconde, rcondv, work, lwork, iwork,
459 $ liwork, bwork, iinfo )
460 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
461 result( 1+rsub ) = ulpinv
462 IF( jtype.NE.22 ) THEN
463 WRITE( nounit, fmt = 9998 )'DGEESX1', iinfo, n, jtype,
464 $ iseed
465 ELSE
466 WRITE( nounit, fmt = 9999 )'DGEESX1', iinfo, n,
467 $ iseed( 1 )
468 END IF
469 info = abs( iinfo )
470 RETURN
471 END IF
472 IF( isort.EQ.0 ) THEN
473 CALL dcopy( n, wr, 1, wrtmp, 1 )
474 CALL dcopy( n, wi, 1, witmp, 1 )
475 END IF
476
477
478
479 result( 1+rsub ) = zero
480 DO 30 j = 1, n - 2
481 DO 20 i = j + 2, n
482 IF( h( i, j ).NE.zero )
483 $ result( 1+rsub ) = ulpinv
484 20 CONTINUE
485 30 CONTINUE
486 DO 40 i = 1, n - 2
487 IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.zero )
488 $ result( 1+rsub ) = ulpinv
489 40 CONTINUE
490 DO 50 i = 1, n - 1
491 IF( h( i+1, i ).NE.zero ) THEN
492 IF( h( i, i ).NE.h( i+1, i+1 ) .OR. h( i, i+1 ).EQ.
493 $ zero .OR. sign( one, h( i+1, i ) ).EQ.
494 $ sign( one, h( i, i+1 ) ) )result( 1+rsub ) = ulpinv
495 END IF
496 50 CONTINUE
497
498
499
500
501
502 CALL dlacpy(
' ', n, n, a, lda, vs1, ldvs )
503
504
505
506 CALL dgemm(
'No transpose',
'No transpose', n, n, n, one, vs,
507 $ ldvs, h, lda, zero, ht, lda )
508
509
510
511 CALL dgemm(
'No transpose',
'Transpose', n, n, n, -one, ht,
512 $ lda, vs, ldvs, one, vs1, ldvs )
513
514 anorm = max(
dlange(
'1', n, n, a, lda, work ), smlnum )
515 wnorm =
dlange(
'1', n, n, vs1, ldvs, work )
516
517 IF( anorm.GT.wnorm ) THEN
518 result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
519 ELSE
520 IF( anorm.LT.one ) THEN
521 result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
522 $ ( n*ulp )
523 ELSE
524 result( 2+rsub ) = min( wnorm / anorm, dble( n ) ) /
525 $ ( n*ulp )
526 END IF
527 END IF
528
529
530
531 CALL dort01(
'Columns', n, n, vs, ldvs, work, lwork,
532 $ result( 3+rsub ) )
533
534
535
536 result( 4+rsub ) = zero
537 DO 60 i = 1, n
538 IF( h( i, i ).NE.wr( i ) )
539 $ result( 4+rsub ) = ulpinv
540 60 CONTINUE
541 IF( n.GT.1 ) THEN
542 IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
543 $ result( 4+rsub ) = ulpinv
544 IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
545 $ result( 4+rsub ) = ulpinv
546 END IF
547 DO 70 i = 1, n - 1
548 IF( h( i+1, i ).NE.zero ) THEN
549 tmp = sqrt( abs( h( i+1, i ) ) )*
550 $ sqrt( abs( h( i, i+1 ) ) )
551 result( 4+rsub ) = max( result( 4+rsub ),
552 $ abs( wi( i )-tmp ) /
553 $ max( ulp*tmp, smlnum ) )
554 result( 4+rsub ) = max( result( 4+rsub ),
555 $ abs( wi( i+1 )+tmp ) /
556 $ max( ulp*tmp, smlnum ) )
557 ELSE IF( i.GT.1 ) THEN
558 IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.zero .AND.
559 $ wi( i ).NE.zero )result( 4+rsub ) = ulpinv
560 END IF
561 70 CONTINUE
562
563
564
565 CALL dlacpy(
'F', n, n, a, lda, ht, lda )
566 CALL dgeesx(
'N', sort,
dslect,
'N', n, ht, lda, sdim, wrt,
567 $ wit, vs, ldvs, rconde, rcondv, work, lwork, iwork,
568 $ liwork, bwork, iinfo )
569 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
570 result( 5+rsub ) = ulpinv
571 IF( jtype.NE.22 ) THEN
572 WRITE( nounit, fmt = 9998 )'DGEESX2', iinfo, n, jtype,
573 $ iseed
574 ELSE
575 WRITE( nounit, fmt = 9999 )'DGEESX2', iinfo, n,
576 $ iseed( 1 )
577 END IF
578 info = abs( iinfo )
579 GO TO 250
580 END IF
581
582 result( 5+rsub ) = zero
583 DO 90 j = 1, n
584 DO 80 i = 1, n
585 IF( h( i, j ).NE.ht( i, j ) )
586 $ result( 5+rsub ) = ulpinv
587 80 CONTINUE
588 90 CONTINUE
589
590
591
592 result( 6+rsub ) = zero
593 DO 100 i = 1, n
594 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
595 $ result( 6+rsub ) = ulpinv
596 100 CONTINUE
597
598
599
600 IF( isort.EQ.1 ) THEN
601 result( 13 ) = zero
602 knteig = 0
603 DO 110 i = 1, n
604 IF(
dslect( wr( i ), wi( i ) ) .OR.
605 $
dslect( wr( i ), -wi( i ) ) )knteig = knteig + 1
606 IF( i.LT.n ) THEN
607 IF( (
dslect( wr( i+1 ), wi( i+1 ) ) .OR.
608 $
dslect( wr( i+1 ), -wi( i+1 ) ) ) .AND.
609 $ ( .NOT.(
dslect( wr( i ),
610 $ wi( i ) ) .OR.
dslect( wr( i ),
611 $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )result( 13 )
612 $ = ulpinv
613 END IF
614 110 CONTINUE
615 IF( sdim.NE.knteig )
616 $ result( 13 ) = ulpinv
617 END IF
618
619 120 CONTINUE
620
621
622
623
624 IF( lwork.GE.n+( n*n ) / 2 ) THEN
625
626
627
628 sort = 'S'
629 result( 14 ) = zero
630 result( 15 ) = zero
631 CALL dlacpy(
'F', n, n, a, lda, ht, lda )
632 CALL dgeesx(
'V', sort,
dslect,
'B', n, ht, lda, sdim1, wrt,
633 $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
634 $ iwork, liwork, bwork, iinfo )
635 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
636 result( 14 ) = ulpinv
637 result( 15 ) = ulpinv
638 IF( jtype.NE.22 ) THEN
639 WRITE( nounit, fmt = 9998 )'DGEESX3', iinfo, n, jtype,
640 $ iseed
641 ELSE
642 WRITE( nounit, fmt = 9999 )'DGEESX3', iinfo, n,
643 $ iseed( 1 )
644 END IF
645 info = abs( iinfo )
646 GO TO 250
647 END IF
648
649
650
651 DO 140 i = 1, n
652 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
653 $ result( 10 ) = ulpinv
654 DO 130 j = 1, n
655 IF( h( i, j ).NE.ht( i, j ) )
656 $ result( 11 ) = ulpinv
657 IF( vs( i, j ).NE.vs1( i, j ) )
658 $ result( 12 ) = ulpinv
659 130 CONTINUE
660 140 CONTINUE
661 IF( sdim.NE.sdim1 )
662 $ result( 13 ) = ulpinv
663
664
665
666 CALL dlacpy(
'F', n, n, a, lda, ht, lda )
667 CALL dgeesx(
'N', sort,
dslect,
'B', n, ht, lda, sdim1, wrt,
668 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
669 $ iwork, liwork, bwork, iinfo )
670 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
671 result( 14 ) = ulpinv
672 result( 15 ) = ulpinv
673 IF( jtype.NE.22 ) THEN
674 WRITE( nounit, fmt = 9998 )'DGEESX4', iinfo, n, jtype,
675 $ iseed
676 ELSE
677 WRITE( nounit, fmt = 9999 )'DGEESX4', iinfo, n,
678 $ iseed( 1 )
679 END IF
680 info = abs( iinfo )
681 GO TO 250
682 END IF
683
684
685
686 IF( rcnde1.NE.rconde )
687 $ result( 14 ) = ulpinv
688 IF( rcndv1.NE.rcondv )
689 $ result( 15 ) = ulpinv
690
691
692
693 DO 160 i = 1, n
694 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
695 $ result( 10 ) = ulpinv
696 DO 150 j = 1, n
697 IF( h( i, j ).NE.ht( i, j ) )
698 $ result( 11 ) = ulpinv
699 IF( vs( i, j ).NE.vs1( i, j ) )
700 $ result( 12 ) = ulpinv
701 150 CONTINUE
702 160 CONTINUE
703 IF( sdim.NE.sdim1 )
704 $ result( 13 ) = ulpinv
705
706
707
708 CALL dlacpy(
'F', n, n, a, lda, ht, lda )
709 CALL dgeesx(
'V', sort,
dslect,
'E', n, ht, lda, sdim1, wrt,
710 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
711 $ iwork, liwork, bwork, iinfo )
712 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
713 result( 14 ) = ulpinv
714 IF( jtype.NE.22 ) THEN
715 WRITE( nounit, fmt = 9998 )'DGEESX5', iinfo, n, jtype,
716 $ iseed
717 ELSE
718 WRITE( nounit, fmt = 9999 )'DGEESX5', iinfo, n,
719 $ iseed( 1 )
720 END IF
721 info = abs( iinfo )
722 GO TO 250
723 END IF
724
725
726
727 IF( rcnde1.NE.rconde )
728 $ result( 14 ) = ulpinv
729
730
731
732 DO 180 i = 1, n
733 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
734 $ result( 10 ) = ulpinv
735 DO 170 j = 1, n
736 IF( h( i, j ).NE.ht( i, j ) )
737 $ result( 11 ) = ulpinv
738 IF( vs( i, j ).NE.vs1( i, j ) )
739 $ result( 12 ) = ulpinv
740 170 CONTINUE
741 180 CONTINUE
742 IF( sdim.NE.sdim1 )
743 $ result( 13 ) = ulpinv
744
745
746
747 CALL dlacpy(
'F', n, n, a, lda, ht, lda )
748 CALL dgeesx(
'N', sort,
dslect,
'E', n, ht, lda, sdim1, wrt,
749 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
750 $ iwork, liwork, bwork, iinfo )
751 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
752 result( 14 ) = ulpinv
753 IF( jtype.NE.22 ) THEN
754 WRITE( nounit, fmt = 9998 )'DGEESX6', iinfo, n, jtype,
755 $ iseed
756 ELSE
757 WRITE( nounit, fmt = 9999 )'DGEESX6', iinfo, n,
758 $ iseed( 1 )
759 END IF
760 info = abs( iinfo )
761 GO TO 250
762 END IF
763
764
765
766 IF( rcnde1.NE.rconde )
767 $ result( 14 ) = ulpinv
768
769
770
771 DO 200 i = 1, n
772 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
773 $ result( 10 ) = ulpinv
774 DO 190 j = 1, n
775 IF( h( i, j ).NE.ht( i, j ) )
776 $ result( 11 ) = ulpinv
777 IF( vs( i, j ).NE.vs1( i, j ) )
778 $ result( 12 ) = ulpinv
779 190 CONTINUE
780 200 CONTINUE
781 IF( sdim.NE.sdim1 )
782 $ result( 13 ) = ulpinv
783
784
785
786 CALL dlacpy(
'F', n, n, a, lda, ht, lda )
787 CALL dgeesx(
'V', sort,
dslect,
'V', n, ht, lda, sdim1, wrt,
788 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
789 $ iwork, liwork, bwork, iinfo )
790 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
791 result( 15 ) = ulpinv
792 IF( jtype.NE.22 ) THEN
793 WRITE( nounit, fmt = 9998 )'DGEESX7', iinfo, n, jtype,
794 $ iseed
795 ELSE
796 WRITE( nounit, fmt = 9999 )'DGEESX7', iinfo, n,
797 $ iseed( 1 )
798 END IF
799 info = abs( iinfo )
800 GO TO 250
801 END IF
802
803
804
805 IF( rcndv1.NE.rcondv )
806 $ result( 15 ) = ulpinv
807
808
809
810 DO 220 i = 1, n
811 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
812 $ result( 10 ) = ulpinv
813 DO 210 j = 1, n
814 IF( h( i, j ).NE.ht( i, j ) )
815 $ result( 11 ) = ulpinv
816 IF( vs( i, j ).NE.vs1( i, j ) )
817 $ result( 12 ) = ulpinv
818 210 CONTINUE
819 220 CONTINUE
820 IF( sdim.NE.sdim1 )
821 $ result( 13 ) = ulpinv
822
823
824
825 CALL dlacpy(
'F', n, n, a, lda, ht, lda )
826 CALL dgeesx(
'N', sort,
dslect,
'V', n, ht, lda, sdim1, wrt,
827 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
828 $ iwork, liwork, bwork, iinfo )
829 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
830 result( 15 ) = ulpinv
831 IF( jtype.NE.22 ) THEN
832 WRITE( nounit, fmt = 9998 )'DGEESX8', iinfo, n, jtype,
833 $ iseed
834 ELSE
835 WRITE( nounit, fmt = 9999 )'DGEESX8', iinfo, n,
836 $ iseed( 1 )
837 END IF
838 info = abs( iinfo )
839 GO TO 250
840 END IF
841
842
843
844 IF( rcndv1.NE.rcondv )
845 $ result( 15 ) = ulpinv
846
847
848
849 DO 240 i = 1, n
850 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
851 $ result( 10 ) = ulpinv
852 DO 230 j = 1, n
853 IF( h( i, j ).NE.ht( i, j ) )
854 $ result( 11 ) = ulpinv
855 IF( vs( i, j ).NE.vs1( i, j ) )
856 $ result( 12 ) = ulpinv
857 230 CONTINUE
858 240 CONTINUE
859 IF( sdim.NE.sdim1 )
860 $ result( 13 ) = ulpinv
861
862 END IF
863
864 250 CONTINUE
865
866
867
868
869 IF( comp ) THEN
870
871
872
873
874
875 seldim = n
876 selopt = 1
877 eps = max( ulp, epsin )
878 DO 260 i = 1, n
879 ipnt( i ) = i
880 selval( i ) = .false.
881 selwr( i ) = wrtmp( i )
882 selwi( i ) = witmp( i )
883 260 CONTINUE
884 DO 280 i = 1, n - 1
885 kmin = i
886 vrmin = wrtmp( i )
887 vimin = witmp( i )
888 DO 270 j = i + 1, n
889 IF( wrtmp( j ).LT.vrmin ) THEN
890 kmin = j
891 vrmin = wrtmp( j )
892 vimin = witmp( j )
893 END IF
894 270 CONTINUE
895 wrtmp( kmin ) = wrtmp( i )
896 witmp( kmin ) = witmp( i )
897 wrtmp( i ) = vrmin
898 witmp( i ) = vimin
899 itmp = ipnt( i )
900 ipnt( i ) = ipnt( kmin )
901 ipnt( kmin ) = itmp
902 280 CONTINUE
903 DO 290 i = 1, nslct
904 selval( ipnt( islct( i ) ) ) = .true.
905 290 CONTINUE
906
907
908
909 CALL dlacpy(
'F', n, n, a, lda, ht, lda )
910 CALL dgeesx(
'N',
'S',
dslect,
'B', n, ht, lda, sdim1, wrt,
911 $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
912 $ iwork, liwork, bwork, iinfo )
913 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
914 result( 16 ) = ulpinv
915 result( 17 ) = ulpinv
916 WRITE( nounit, fmt = 9999 )'DGEESX9', iinfo, n, iseed( 1 )
917 info = abs( iinfo )
918 GO TO 300
919 END IF
920
921
922
923
924 anorm =
dlange(
'1', n, n, a, lda, work )
925 v = max( dble( n )*eps*anorm, smlnum )
926 IF( anorm.EQ.zero )
927 $ v = one
928 IF( v.GT.rcondv ) THEN
929 tol = one
930 ELSE
931 tol = v / rcondv
932 END IF
933 IF( v.GT.rcdvin ) THEN
934 tolin = one
935 ELSE
936 tolin = v / rcdvin
937 END IF
938 tol = max( tol, smlnum / eps )
939 tolin = max( tolin, smlnum / eps )
940 IF( eps*( rcdein-tolin ).GT.rconde+tol ) THEN
941 result( 16 ) = ulpinv
942 ELSE IF( rcdein-tolin.GT.rconde+tol ) THEN
943 result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
944 ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) ) THEN
945 result( 16 ) = ulpinv
946 ELSE IF( rcdein+tolin.LT.rconde-tol ) THEN
947 result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
948 ELSE
949 result( 16 ) = one
950 END IF
951
952
953
954
955 IF( v.GT.rcondv*rconde ) THEN
956 tol = rcondv
957 ELSE
958 tol = v / rconde
959 END IF
960 IF( v.GT.rcdvin*rcdein ) THEN
961 tolin = rcdvin
962 ELSE
963 tolin = v / rcdein
964 END IF
965 tol = max( tol, smlnum / eps )
966 tolin = max( tolin, smlnum / eps )
967 IF( eps*( rcdvin-tolin ).GT.rcondv+tol ) THEN
968 result( 17 ) = ulpinv
969 ELSE IF( rcdvin-tolin.GT.rcondv+tol ) THEN
970 result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
971 ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) ) THEN
972 result( 17 ) = ulpinv
973 ELSE IF( rcdvin+tolin.LT.rcondv-tol ) THEN
974 result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
975 ELSE
976 result( 17 ) = one
977 END IF
978
979 300 CONTINUE
980
981 END IF
982
983 9999 FORMAT( ' DGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
984 $ i6, ', INPUT EXAMPLE NUMBER = ', i4 )
985 9998 FORMAT( ' DGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
986 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
987
988 RETURN
989
990
991
subroutine xerbla(srname, info)
subroutine dort01(rowcol, m, n, u, ldu, work, lwork, resid)
DORT01
logical function dslect(zr, zi)
DSLECT
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
subroutine dgeesx(jobvs, sort, select, sense, n, a, lda, sdim, wr, wi, vs, ldvs, rconde, rcondv, work, lwork, iwork, liwork, bwork, info)
DGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...