401
402
403
404
405
406 IMPLICIT NONE
407
408
409 INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES,
410 $ NTYPES
411 DOUBLE PRECISION THRESH
412
413
414 LOGICAL DOTYPE( * )
415 INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
416 DOUBLE PRECISION E( * ), RWORK( * ), S( * ), SSAV( * )
417 COMPLEX*16 A( LDA, * ), ASAV( LDA, * ), U( LDU, * ),
418 $ USAV( LDU, * ), VT( LDVT, * ),
419 $ VTSAV( LDVT, * ), WORK( * )
420
421
422
423
424
425 DOUBLE PRECISION ZERO, ONE, TWO, HALF
426 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
427 $ half = 0.5d0 )
428 COMPLEX*16 CZERO, CONE
429 parameter( czero = ( 0.0d+0, 0.0d+0 ),
430 $ cone = ( 1.0d+0, 0.0d+0 ) )
431 INTEGER MAXTYP
432 parameter( maxtyp = 5 )
433
434
435 LOGICAL BADMM, BADNN
436 CHARACTER JOBQ, JOBU, JOBVT, RANGE
437 INTEGER I, IINFO, IJQ, IJU, IJVT, IL, IU, ITEMP,
438 $ IWSPC, IWTMP, J, JSIZE, JTYPE, LSWORK, M,
439 $ MINWRK, MMAX, MNMAX, MNMIN, MTYPES, N,
440 $ NERRS, NFAIL, NMAX, NS, NSI, NSV, NTEST,
441 $ NTESTF, NTESTT, LRWORK
442 DOUBLE PRECISION ANORM, DIF, DIV, OVFL, RTUNFL, ULP, ULPINV,
443 $ UNFL, VL, VU
444
445
446 INTEGER LIWORK, NUMRANK
447
448
449 CHARACTER CJOB( 4 ), CJOBR( 3 ), CJOBV( 2 )
450 INTEGER IOLDSD( 4 ), ISEED2( 4 )
451 DOUBLE PRECISION RESULT( 39 )
452
453
454 DOUBLE PRECISION DLAMCH, DLARND
456
457
461
462
463 INTRINSIC abs, dble, max, min
464
465
466 CHARACTER*32 SRNAMT
467
468
469 COMMON / srnamc / srnamt
470
471
472 DATA cjob / 'N', 'O', 'S', 'A' /
473 DATA cjobr / 'A', 'V', 'I' /
474 DATA cjobv / 'N', 'V' /
475
476
477
478
479
480 info = 0
481
482
483
484 nerrs = 0
485 ntestt = 0
486 ntestf = 0
487 badmm = .false.
488 badnn = .false.
489 mmax = 1
490 nmax = 1
491 mnmax = 1
492 minwrk = 1
493 DO 10 j = 1, nsizes
494 mmax = max( mmax, mm( j ) )
495 IF( mm( j ).LT.0 )
496 $ badmm = .true.
497 nmax = max( nmax, nn( j ) )
498 IF( nn( j ).LT.0 )
499 $ badnn = .true.
500 mnmax = max( mnmax, min( mm( j ), nn( j ) ) )
501 minwrk = max( minwrk, max( 3*min( mm( j ),
502 $ nn( j ) )+max( mm( j ), nn( j ) )**2, 5*min( mm( j ),
503 $ nn( j ) ), 3*max( mm( j ), nn( j ) ) ) )
504 10 CONTINUE
505
506
507
508 IF( nsizes.LT.0 ) THEN
509 info = -1
510 ELSE IF( badmm ) THEN
511 info = -2
512 ELSE IF( badnn ) THEN
513 info = -3
514 ELSE IF( ntypes.LT.0 ) THEN
515 info = -4
516 ELSE IF( lda.LT.max( 1, mmax ) ) THEN
517 info = -10
518 ELSE IF( ldu.LT.max( 1, mmax ) ) THEN
519 info = -12
520 ELSE IF( ldvt.LT.max( 1, nmax ) ) THEN
521 info = -14
522 ELSE IF( minwrk.GT.lwork ) THEN
523 info = -21
524 END IF
525
526 IF( info.NE.0 ) THEN
527 CALL xerbla(
'ZDRVBD', -info )
528 RETURN
529 END IF
530
531
532
533 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
534 $ RETURN
535
536
537
539 ovfl = one / unfl
541 ulpinv = one / ulp
542 rtunfl = sqrt( unfl )
543
544
545
546 nerrs = 0
547
548 DO 230 jsize = 1, nsizes
549 m = mm( jsize )
550 n = nn( jsize )
551 mnmin = min( m, n )
552
553 IF( nsizes.NE.1 ) THEN
554 mtypes = min( maxtyp, ntypes )
555 ELSE
556 mtypes = min( maxtyp+1, ntypes )
557 END IF
558
559 DO 220 jtype = 1, mtypes
560 IF( .NOT.dotype( jtype ) )
561 $ GO TO 220
562 ntest = 0
563
564 DO 20 j = 1, 4
565 ioldsd( j ) = iseed( j )
566 20 CONTINUE
567
568
569
570 IF( mtypes.GT.maxtyp )
571 $ GO TO 50
572
573 IF( jtype.EQ.1 ) THEN
574
575
576
577 CALL zlaset(
'Full', m, n, czero, czero, a, lda )
578 DO 30 i = 1, min( m, n )
579 s( i ) = zero
580 30 CONTINUE
581
582 ELSE IF( jtype.EQ.2 ) THEN
583
584
585
586 CALL zlaset(
'Full', m, n, czero, cone, a, lda )
587 DO 40 i = 1, min( m, n )
588 s( i ) = one
589 40 CONTINUE
590
591 ELSE
592
593
594
595 IF( jtype.EQ.3 )
596 $ anorm = one
597 IF( jtype.EQ.4 )
598 $ anorm = unfl / ulp
599 IF( jtype.EQ.5 )
600 $ anorm = ovfl*ulp
601 CALL zlatms( m, n,
'U', iseed,
'N', s, 4, dble( mnmin ),
602 $ anorm, m-1, n-1, 'N', a, lda, work, iinfo )
603 IF( iinfo.NE.0 ) THEN
604 WRITE( nounit, fmt = 9996 )'Generator', iinfo, m, n,
605 $ jtype, ioldsd
606 info = abs( iinfo )
607 RETURN
608 END IF
609 END IF
610
611 50 CONTINUE
612 CALL zlacpy(
'F', m, n, a, lda, asav, lda )
613
614
615
616 DO 210 iwspc = 1, 4
617
618
619
620 iwtmp = 2*min( m, n )+max( m, n )
621 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
622 lswork = min( lswork, lwork )
623 lswork = max( lswork, 1 )
624 IF( iwspc.EQ.4 )
625 $ lswork = lwork
626
627 DO 60 j = 1, 35
628 result( j ) = -one
629 60 CONTINUE
630
631
632
633 IF( iwspc.GT.1 )
634 $
CALL zlacpy(
'F', m, n, asav, lda, a, lda )
635 srnamt = 'ZGESVD'
636 CALL zgesvd(
'A',
'A', m, n, a, lda, ssav, usav, ldu,
637 $ vtsav, ldvt, work, lswork, rwork, iinfo )
638 IF( iinfo.NE.0 ) THEN
639 WRITE( nounit, fmt = 9995 )'GESVD', iinfo, m, n,
640 $ jtype, lswork, ioldsd
641 info = abs( iinfo )
642 RETURN
643 END IF
644
645
646
647 CALL zbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
648 $ vtsav, ldvt, work, rwork, result( 1 ) )
649 IF( m.NE.0 .AND. n.NE.0 ) THEN
650 CALL zunt01(
'Columns', mnmin, m, usav, ldu, work,
651 $ lwork, rwork, result( 2 ) )
652 CALL zunt01(
'Rows', mnmin, n, vtsav, ldvt, work,
653 $ lwork, rwork, result( 3 ) )
654 END IF
655 result( 4 ) = 0
656 DO 70 i = 1, mnmin - 1
657 IF( ssav( i ).LT.ssav( i+1 ) )
658 $ result( 4 ) = ulpinv
659 IF( ssav( i ).LT.zero )
660 $ result( 4 ) = ulpinv
661 70 CONTINUE
662 IF( mnmin.GE.1 ) THEN
663 IF( ssav( mnmin ).LT.zero )
664 $ result( 4 ) = ulpinv
665 END IF
666
667
668
669 result( 5 ) = zero
670 result( 6 ) = zero
671 result( 7 ) = zero
672 DO 100 iju = 0, 3
673 DO 90 ijvt = 0, 3
674 IF( ( iju.EQ.3 .AND. ijvt.EQ.3 ) .OR.
675 $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )GO TO 90
676 jobu = cjob( iju+1 )
677 jobvt = cjob( ijvt+1 )
678 CALL zlacpy(
'F', m, n, asav, lda, a, lda )
679 srnamt = 'ZGESVD'
680 CALL zgesvd( jobu, jobvt, m, n, a, lda, s, u, ldu,
681 $ vt, ldvt, work, lswork, rwork, iinfo )
682
683
684
685 dif = zero
686 IF( m.GT.0 .AND. n.GT.0 ) THEN
687 IF( iju.EQ.1 ) THEN
688 CALL zunt03(
'C', m, mnmin, m, mnmin, usav,
689 $ ldu, a, lda, work, lwork, rwork,
690 $ dif, iinfo )
691 ELSE IF( iju.EQ.2 ) THEN
692 CALL zunt03(
'C', m, mnmin, m, mnmin, usav,
693 $ ldu, u, ldu, work, lwork, rwork,
694 $ dif, iinfo )
695 ELSE IF( iju.EQ.3 ) THEN
696 CALL zunt03(
'C', m, m, m, mnmin, usav, ldu,
697 $ u, ldu, work, lwork, rwork, dif,
698 $ iinfo )
699 END IF
700 END IF
701 result( 5 ) = max( result( 5 ), dif )
702
703
704
705 dif = zero
706 IF( m.GT.0 .AND. n.GT.0 ) THEN
707 IF( ijvt.EQ.1 ) THEN
708 CALL zunt03(
'R', n, mnmin, n, mnmin, vtsav,
709 $ ldvt, a, lda, work, lwork,
710 $ rwork, dif, iinfo )
711 ELSE IF( ijvt.EQ.2 ) THEN
712 CALL zunt03(
'R', n, mnmin, n, mnmin, vtsav,
713 $ ldvt, vt, ldvt, work, lwork,
714 $ rwork, dif, iinfo )
715 ELSE IF( ijvt.EQ.3 ) THEN
716 CALL zunt03(
'R', n, n, n, mnmin, vtsav,
717 $ ldvt, vt, ldvt, work, lwork,
718 $ rwork, dif, iinfo )
719 END IF
720 END IF
721 result( 6 ) = max( result( 6 ), dif )
722
723
724
725 dif = zero
726 div = max( dble( mnmin )*ulp*s( 1 ),
727 $
dlamch(
'Safe minimum' ) )
728 DO 80 i = 1, mnmin - 1
729 IF( ssav( i ).LT.ssav( i+1 ) )
730 $ dif = ulpinv
731 IF( ssav( i ).LT.zero )
732 $ dif = ulpinv
733 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
734 80 CONTINUE
735 result( 7 ) = max( result( 7 ), dif )
736 90 CONTINUE
737 100 CONTINUE
738
739
740
741 iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
742 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
743 lswork = min( lswork, lwork )
744 lswork = max( lswork, 1 )
745 IF( iwspc.EQ.4 )
746 $ lswork = lwork
747
748
749
750 CALL zlacpy(
'F', m, n, asav, lda, a, lda )
751 srnamt = 'ZGESDD'
752 CALL zgesdd(
'A', m, n, a, lda, ssav, usav, ldu, vtsav,
753 $ ldvt, work, lswork, rwork, iwork, iinfo )
754 IF( iinfo.NE.0 ) THEN
755 WRITE( nounit, fmt = 9995 )'GESDD', iinfo, m, n,
756 $ jtype, lswork, ioldsd
757 info = abs( iinfo )
758 RETURN
759 END IF
760
761
762
763 CALL zbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
764 $ vtsav, ldvt, work, rwork, result( 8 ) )
765 IF( m.NE.0 .AND. n.NE.0 ) THEN
766 CALL zunt01(
'Columns', mnmin, m, usav, ldu, work,
767 $ lwork, rwork, result( 9 ) )
768 CALL zunt01(
'Rows', mnmin, n, vtsav, ldvt, work,
769 $ lwork, rwork, result( 10 ) )
770 END IF
771 result( 11 ) = 0
772 DO 110 i = 1, mnmin - 1
773 IF( ssav( i ).LT.ssav( i+1 ) )
774 $ result( 11 ) = ulpinv
775 IF( ssav( i ).LT.zero )
776 $ result( 11 ) = ulpinv
777 110 CONTINUE
778 IF( mnmin.GE.1 ) THEN
779 IF( ssav( mnmin ).LT.zero )
780 $ result( 11 ) = ulpinv
781 END IF
782
783
784
785 result( 12 ) = zero
786 result( 13 ) = zero
787 result( 14 ) = zero
788 DO 130 ijq = 0, 2
789 jobq = cjob( ijq+1 )
790 CALL zlacpy(
'F', m, n, asav, lda, a, lda )
791 srnamt = 'ZGESDD'
792 CALL zgesdd( jobq, m, n, a, lda, s, u, ldu, vt, ldvt,
793 $ work, lswork, rwork, iwork, iinfo )
794
795
796
797 dif = zero
798 IF( m.GT.0 .AND. n.GT.0 ) THEN
799 IF( ijq.EQ.1 ) THEN
800 IF( m.GE.n ) THEN
801 CALL zunt03(
'C', m, mnmin, m, mnmin, usav,
802 $ ldu, a, lda, work, lwork, rwork,
803 $ dif, iinfo )
804 ELSE
805 CALL zunt03(
'C', m, mnmin, m, mnmin, usav,
806 $ ldu, u, ldu, work, lwork, rwork,
807 $ dif, iinfo )
808 END IF
809 ELSE IF( ijq.EQ.2 ) THEN
810 CALL zunt03(
'C', m, mnmin, m, mnmin, usav, ldu,
811 $ u, ldu, work, lwork, rwork, dif,
812 $ iinfo )
813 END IF
814 END IF
815 result( 12 ) = max( result( 12 ), dif )
816
817
818
819 dif = zero
820 IF( m.GT.0 .AND. n.GT.0 ) THEN
821 IF( ijq.EQ.1 ) THEN
822 IF( m.GE.n ) THEN
823 CALL zunt03(
'R', n, mnmin, n, mnmin, vtsav,
824 $ ldvt, vt, ldvt, work, lwork,
825 $ rwork, dif, iinfo )
826 ELSE
827 CALL zunt03(
'R', n, mnmin, n, mnmin, vtsav,
828 $ ldvt, a, lda, work, lwork,
829 $ rwork, dif, iinfo )
830 END IF
831 ELSE IF( ijq.EQ.2 ) THEN
832 CALL zunt03(
'R', n, mnmin, n, mnmin, vtsav,
833 $ ldvt, vt, ldvt, work, lwork, rwork,
834 $ dif, iinfo )
835 END IF
836 END IF
837 result( 13 ) = max( result( 13 ), dif )
838
839
840
841 dif = zero
842 div = max( dble( mnmin )*ulp*s( 1 ),
843 $
dlamch(
'Safe minimum' ) )
844 DO 120 i = 1, mnmin - 1
845 IF( ssav( i ).LT.ssav( i+1 ) )
846 $ dif = ulpinv
847 IF( ssav( i ).LT.zero )
848 $ dif = ulpinv
849 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
850 120 CONTINUE
851 result( 14 ) = max( result( 14 ), dif )
852 130 CONTINUE
853
854
855
856
857 result( 36 ) = zero
858 result( 37 ) = zero
859 result( 38 ) = zero
860 result( 39 ) = zero
861
862 IF( m.GE.n ) THEN
863 iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
864 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
865 lswork = min( lswork, lwork )
866 lswork = max( lswork, 1 )
867 IF( iwspc.EQ.4 )
868 $ lswork = lwork
869
870 CALL zlacpy(
'F', m, n, asav, lda, a, lda )
871 srnamt = 'ZGESVDQ'
872
873 lrwork = max(2, m, 5*n)
874 liwork = max( n, 1 )
875 CALL zgesvdq(
'H',
'N',
'N',
'A',
'A',
876 $ m, n, a, lda, ssav, usav, ldu,
877 $ vtsav, ldvt, numrank, iwork, liwork,
878 $ work, lwork, rwork, lrwork, iinfo )
879
880 IF( iinfo.NE.0 ) THEN
881 WRITE( nounit, fmt = 9995 )'ZGESVDQ', iinfo, m, n,
882 $ jtype, lswork, ioldsd
883 info = abs( iinfo )
884 RETURN
885 END IF
886
887
888
889 CALL zbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
890 $ vtsav, ldvt, work, rwork, result( 36 ) )
891 IF( m.NE.0 .AND. n.NE.0 ) THEN
892 CALL zunt01(
'Columns', m, m, usav, ldu, work,
893 $ lwork, rwork, result( 37 ) )
894 CALL zunt01(
'Rows', n, n, vtsav, ldvt, work,
895 $ lwork, rwork, result( 38 ) )
896 END IF
897 result( 39 ) = zero
898 DO 199 i = 1, mnmin - 1
899 IF( ssav( i ).LT.ssav( i+1 ) )
900 $ result( 39 ) = ulpinv
901 IF( ssav( i ).LT.zero )
902 $ result( 39 ) = ulpinv
903 199 CONTINUE
904 IF( mnmin.GE.1 ) THEN
905 IF( ssav( mnmin ).LT.zero )
906 $ result( 39 ) = ulpinv
907 END IF
908 END IF
909
910
911
912
913 result( 15 ) = zero
914 result( 16 ) = zero
915 result( 17 ) = zero
916 result( 18 ) = zero
917
918 IF( m.GE.n ) THEN
919 iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
920 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
921 lswork = min( lswork, lwork )
922 lswork = max( lswork, 1 )
923 lrwork = max(6,n)
924 IF( iwspc.EQ.4 )
925 $ lswork = lwork
926
927 CALL zlacpy(
'F', m, n, asav, lda, usav, lda )
928 srnamt = 'ZGESVJ'
929 CALL zgesvj(
'G',
'U',
'V', m, n, usav, lda, ssav,
930 & 0, a, ldvt, work, lwork, rwork,
931 & lrwork, iinfo )
932
933
934
935 DO j=1,n
936 DO i=1,n
937 vtsav(j,i) = conjg(a(i,j))
938 END DO
939 END DO
940
941 IF( iinfo.NE.0 ) THEN
942 WRITE( nounit, fmt = 9995 )'GESVJ', iinfo, m, n,
943 $ jtype, lswork, ioldsd
944 info = abs( iinfo )
945 RETURN
946 END IF
947
948
949
950 CALL zbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
951 $ vtsav, ldvt, work, rwork, result( 15 ) )
952 IF( m.NE.0 .AND. n.NE.0 ) THEN
953 CALL zunt01(
'Columns', m, m, usav, ldu, work,
954 $ lwork, rwork, result( 16 ) )
955 CALL zunt01(
'Rows', n, n, vtsav, ldvt, work,
956 $ lwork, rwork, result( 17 ) )
957 END IF
958 result( 18 ) = zero
959 DO 131 i = 1, mnmin - 1
960 IF( ssav( i ).LT.ssav( i+1 ) )
961 $ result( 18 ) = ulpinv
962 IF( ssav( i ).LT.zero )
963 $ result( 18 ) = ulpinv
964 131 CONTINUE
965 IF( mnmin.GE.1 ) THEN
966 IF( ssav( mnmin ).LT.zero )
967 $ result( 18 ) = ulpinv
968 END IF
969 END IF
970
971
972
973
974 result( 19 ) = zero
975 result( 20 ) = zero
976 result( 21 ) = zero
977 result( 22 ) = zero
978 IF( m.GE.n ) THEN
979 iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
980 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
981 lswork = min( lswork, lwork )
982 lswork = max( lswork, 1 )
983 IF( iwspc.EQ.4 )
984 $ lswork = lwork
985 lrwork = max( 7, n + 2*m)
986
987 CALL zlacpy(
'F', m, n, asav, lda, vtsav, lda )
988 srnamt = 'ZGEJSV'
989 CALL zgejsv(
'G',
'U',
'V',
'R',
'N',
'N',
990 & m, n, vtsav, lda, ssav, usav, ldu, a, ldvt,
991 & work, lwork, rwork,
992 & lrwork, iwork, iinfo )
993
994
995
996 DO 133 j=1,n
997 DO 132 i=1,n
998 vtsav(j,i) = conjg(a(i,j))
999 132 END DO
1000 133 END DO
1001
1002 IF( iinfo.NE.0 ) THEN
1003 WRITE( nounit, fmt = 9995 )'GEJSV', iinfo, m, n,
1004 $ jtype, lswork, ioldsd
1005 info = abs( iinfo )
1006 RETURN
1007 END IF
1008
1009
1010
1011 CALL zbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
1012 $ vtsav, ldvt, work, rwork, result( 19 ) )
1013 IF( m.NE.0 .AND. n.NE.0 ) THEN
1014 CALL zunt01(
'Columns', m, m, usav, ldu, work,
1015 $ lwork, rwork, result( 20 ) )
1016 CALL zunt01(
'Rows', n, n, vtsav, ldvt, work,
1017 $ lwork, rwork, result( 21 ) )
1018 END IF
1019 result( 22 ) = zero
1020 DO 134 i = 1, mnmin - 1
1021 IF( ssav( i ).LT.ssav( i+1 ) )
1022 $ result( 22 ) = ulpinv
1023 IF( ssav( i ).LT.zero )
1024 $ result( 22 ) = ulpinv
1025 134 CONTINUE
1026 IF( mnmin.GE.1 ) THEN
1027 IF( ssav( mnmin ).LT.zero )
1028 $ result( 22 ) = ulpinv
1029 END IF
1030 END IF
1031
1032
1033
1034
1035
1036 CALL zlacpy(
'F', m, n, asav, lda, a, lda )
1037 srnamt = 'ZGESVDX'
1038 CALL zgesvdx(
'V',
'V',
'A', m, n, a, lda,
1039 $ vl, vu, il, iu, ns, ssav, usav, ldu,
1040 $ vtsav, ldvt, work, lwork, rwork,
1041 $ iwork, iinfo )
1042 IF( iinfo.NE.0 ) THEN
1043 WRITE( nounit, fmt = 9995 )'GESVDX', iinfo, m, n,
1044 $ jtype, lswork, ioldsd
1045 info = abs( iinfo )
1046 RETURN
1047 END IF
1048
1049
1050
1051 result( 23 ) = zero
1052 result( 24 ) = zero
1053 result( 25 ) = zero
1054 CALL zbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
1055 $ vtsav, ldvt, work, rwork, result( 23 ) )
1056 IF( m.NE.0 .AND. n.NE.0 ) THEN
1057 CALL zunt01(
'Columns', mnmin, m, usav, ldu, work,
1058 $ lwork, rwork, result( 24 ) )
1059 CALL zunt01(
'Rows', mnmin, n, vtsav, ldvt, work,
1060 $ lwork, rwork, result( 25 ) )
1061 END IF
1062 result( 26 ) = zero
1063 DO 140 i = 1, mnmin - 1
1064 IF( ssav( i ).LT.ssav( i+1 ) )
1065 $ result( 26 ) = ulpinv
1066 IF( ssav( i ).LT.zero )
1067 $ result( 26 ) = ulpinv
1068 140 CONTINUE
1069 IF( mnmin.GE.1 ) THEN
1070 IF( ssav( mnmin ).LT.zero )
1071 $ result( 26 ) = ulpinv
1072 END IF
1073
1074
1075
1076 result( 27 ) = zero
1077 result( 28 ) = zero
1078 result( 29 ) = zero
1079 DO 170 iju = 0, 1
1080 DO 160 ijvt = 0, 1
1081 IF( ( iju.EQ.0 .AND. ijvt.EQ.0 ) .OR.
1082 $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) ) GO TO 160
1083 jobu = cjobv( iju+1 )
1084 jobvt = cjobv( ijvt+1 )
1085 range = cjobr( 1 )
1086 CALL zlacpy(
'F', m, n, asav, lda, a, lda )
1087 srnamt = 'ZGESVDX'
1088 CALL zgesvdx( jobu, jobvt,
'A', m, n, a, lda,
1089 $ vl, vu, il, iu, ns, ssav, u, ldu,
1090 $ vt, ldvt, work, lwork, rwork,
1091 $ iwork, iinfo )
1092
1093
1094
1095 dif = zero
1096 IF( m.GT.0 .AND. n.GT.0 ) THEN
1097 IF( iju.EQ.1 ) THEN
1098 CALL zunt03(
'C', m, mnmin, m, mnmin, usav,
1099 $ ldu, u, ldu, work, lwork, rwork,
1100 $ dif, iinfo )
1101 END IF
1102 END IF
1103 result( 27 ) = max( result( 27 ), dif )
1104
1105
1106
1107 dif = zero
1108 IF( m.GT.0 .AND. n.GT.0 ) THEN
1109 IF( ijvt.EQ.1 ) THEN
1110 CALL zunt03(
'R', n, mnmin, n, mnmin, vtsav,
1111 $ ldvt, vt, ldvt, work, lwork,
1112 $ rwork, dif, iinfo )
1113 END IF
1114 END IF
1115 result( 28 ) = max( result( 28 ), dif )
1116
1117
1118
1119 dif = zero
1120 div = max( dble( mnmin )*ulp*s( 1 ),
1121 $
dlamch(
'Safe minimum' ) )
1122 DO 150 i = 1, mnmin - 1
1123 IF( ssav( i ).LT.ssav( i+1 ) )
1124 $ dif = ulpinv
1125 IF( ssav( i ).LT.zero )
1126 $ dif = ulpinv
1127 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
1128 150 CONTINUE
1129 result( 29) = max( result( 29 ), dif )
1130 160 CONTINUE
1131 170 CONTINUE
1132
1133
1134
1135 DO 180 i = 1, 4
1136 iseed2( i ) = iseed( i )
1137 180 CONTINUE
1138 IF( mnmin.LE.1 ) THEN
1139 il = 1
1140 iu = max( 1, mnmin )
1141 ELSE
1142 il = 1 + int( ( mnmin-1 )*
dlarnd( 1, iseed2 ) )
1143 iu = 1 + int( ( mnmin-1 )*
dlarnd( 1, iseed2 ) )
1144 IF( iu.LT.il ) THEN
1145 itemp = iu
1146 iu = il
1147 il = itemp
1148 END IF
1149 END IF
1150 CALL zlacpy(
'F', m, n, asav, lda, a, lda )
1151 srnamt = 'ZGESVDX'
1152 CALL zgesvdx(
'V',
'V',
'I', m, n, a, lda,
1153 $ vl, vu, il, iu, nsi, s, u, ldu,
1154 $ vt, ldvt, work, lwork, rwork,
1155 $ iwork, iinfo )
1156 IF( iinfo.NE.0 ) THEN
1157 WRITE( nounit, fmt = 9995 )'GESVDX', iinfo, m, n,
1158 $ jtype, lswork, ioldsd
1159 info = abs( iinfo )
1160 RETURN
1161 END IF
1162
1163 result( 30 ) = zero
1164 result( 31 ) = zero
1165 result( 32 ) = zero
1166 CALL zbdt05( m, n, asav, lda, s, nsi, u, ldu,
1167 $ vt, ldvt, work, result( 30 ) )
1168 IF( m.NE.0 .AND. n.NE.0 ) THEN
1169 CALL zunt01(
'Columns', m, nsi, u, ldu, work,
1170 $ lwork, rwork, result( 31 ) )
1171 CALL zunt01(
'Rows', nsi, n, vt, ldvt, work,
1172 $ lwork, rwork, result( 32 ) )
1173 END IF
1174
1175
1176
1177 IF( mnmin.GT.0 .AND. nsi.GT.1 ) THEN
1178 IF( il.NE.1 ) THEN
1179 vu = ssav( il ) +
1180 $ max( half*abs( ssav( il )-ssav( il-1 ) ),
1181 $ ulp*anorm, two*rtunfl )
1182 ELSE
1183 vu = ssav( 1 ) +
1184 $ max( half*abs( ssav( ns )-ssav( 1 ) ),
1185 $ ulp*anorm, two*rtunfl )
1186 END IF
1187 IF( iu.NE.ns ) THEN
1188 vl = ssav( iu ) - max( ulp*anorm, two*rtunfl,
1189 $ half*abs( ssav( iu+1 )-ssav( iu ) ) )
1190 ELSE
1191 vl = ssav( ns ) - max( ulp*anorm, two*rtunfl,
1192 $ half*abs( ssav( ns )-ssav( 1 ) ) )
1193 END IF
1194 vl = max( vl,zero )
1195 vu = max( vu,zero )
1196 IF( vl.GE.vu ) vu = max( vu*2, vu+vl+half )
1197 ELSE
1198 vl = zero
1199 vu = one
1200 END IF
1201 CALL zlacpy(
'F', m, n, asav, lda, a, lda )
1202 srnamt = 'ZGESVDX'
1203 CALL zgesvdx(
'V',
'V',
'V', m, n, a, lda,
1204 $ vl, vu, il, iu, nsv, s, u, ldu,
1205 $ vt, ldvt, work, lwork, rwork,
1206 $ iwork, iinfo )
1207 IF( iinfo.NE.0 ) THEN
1208 WRITE( nounit, fmt = 9995 )'GESVDX', iinfo, m, n,
1209 $ jtype, lswork, ioldsd
1210 info = abs( iinfo )
1211 RETURN
1212 END IF
1213
1214 result( 33 ) = zero
1215 result( 34 ) = zero
1216 result( 35 ) = zero
1217 CALL zbdt05( m, n, asav, lda, s, nsv, u, ldu,
1218 $ vt, ldvt, work, result( 33 ) )
1219 IF( m.NE.0 .AND. n.NE.0 ) THEN
1220 CALL zunt01(
'Columns', m, nsv, u, ldu, work,
1221 $ lwork, rwork, result( 34 ) )
1222 CALL zunt01(
'Rows', nsv, n, vt, ldvt, work,
1223 $ lwork, rwork, result( 35 ) )
1224 END IF
1225
1226
1227
1228 ntest = 0
1229 nfail = 0
1230 DO 190 j = 1, 39
1231 IF( result( j ).GE.zero )
1232 $ ntest = ntest + 1
1233 IF( result( j ).GE.thresh )
1234 $ nfail = nfail + 1
1235 190 CONTINUE
1236
1237 IF( nfail.GT.0 )
1238 $ ntestf = ntestf + 1
1239 IF( ntestf.EQ.1 ) THEN
1240 WRITE( nounit, fmt = 9999 )
1241 WRITE( nounit, fmt = 9998 )thresh
1242 ntestf = 2
1243 END IF
1244
1245 DO 200 j = 1, 39
1246 IF( result( j ).GE.thresh ) THEN
1247 WRITE( nounit, fmt = 9997 )m, n, jtype, iwspc,
1248 $ ioldsd, j, result( j )
1249 END IF
1250 200 CONTINUE
1251
1252 nerrs = nerrs + nfail
1253 ntestt = ntestt + ntest
1254
1255 210 CONTINUE
1256
1257 220 CONTINUE
1258 230 CONTINUE
1259
1260
1261
1262 CALL alasvm(
'ZBD', nounit, nerrs, ntestt, 0 )
1263
1264 9999 FORMAT( ' SVD -- Complex Singular Value Decomposition Driver ',
1265 $ / ' Matrix types (see ZDRVBD for details):',
1266 $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
1267 $ / ' 3 = Evenly spaced singular values near 1',
1268 $ / ' 4 = Evenly spaced singular values near underflow',
1269 $ / ' 5 = Evenly spaced singular values near overflow',
1270 $ / / ' Tests performed: ( A is dense, U and V are unitary,',
1271 $ / 19x, ' S is an array, and Upartial, VTpartial, and',
1272 $ / 19x, ' Spartial are partially computed U, VT and S),', / )
1273 9998 FORMAT( ' Tests performed with Test Threshold = ', f8.2,
1274 $ / ' ZGESVD: ', /
1275 $ ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1276 $ / ' 2 = | I - U**T U | / ( M ulp ) ',
1277 $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
1278 $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
1279 $ ' decreasing order, else 1/ulp',
1280 $ / ' 5 = | U - Upartial | / ( M ulp )',
1281 $ / ' 6 = | VT - VTpartial | / ( N ulp )',
1282 $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
1283 $ / ' ZGESDD: ', /
1284 $ ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1285 $ / ' 9 = | I - U**T U | / ( M ulp ) ',
1286 $ / '10 = | I - VT VT**T | / ( N ulp ) ',
1287 $ / '11 = 0 if S contains min(M,N) nonnegative values in',
1288 $ ' decreasing order, else 1/ulp',
1289 $ / '12 = | U - Upartial | / ( M ulp )',
1290 $ / '13 = | VT - VTpartial | / ( N ulp )',
1291 $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )',
1292 $ / ' ZGESVJ: ', /
1293 $ / '15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1294 $ / '16 = | I - U**T U | / ( M ulp ) ',
1295 $ / '17 = | I - VT VT**T | / ( N ulp ) ',
1296 $ / '18 = 0 if S contains min(M,N) nonnegative values in',
1297 $ ' decreasing order, else 1/ulp',
1298 $ / ' ZGESJV: ', /
1299 $ / '19 = | A - U diag(S) VT | / ( |A| max(M,N) ulp )',
1300 $ / '20 = | I - U**T U | / ( M ulp ) ',
1301 $ / '21 = | I - VT VT**T | / ( N ulp ) ',
1302 $ / '22 = 0 if S contains min(M,N) nonnegative values in',
1303 $ ' decreasing order, else 1/ulp',
1304 $ / ' ZGESVDX(V,V,A): ', /
1305 $ '23 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1306 $ / '24 = | I - U**T U | / ( M ulp ) ',
1307 $ / '25 = | I - VT VT**T | / ( N ulp ) ',
1308 $ / '26 = 0 if S contains min(M,N) nonnegative values in',
1309 $ ' decreasing order, else 1/ulp',
1310 $ / '27 = | U - Upartial | / ( M ulp )',
1311 $ / '28 = | VT - VTpartial | / ( N ulp )',
1312 $ / '29 = | S - Spartial | / ( min(M,N) ulp |S| )',
1313 $ / ' ZGESVDX(V,V,I): ',
1314 $ / '30 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp )',
1315 $ / '31 = | I - U**T U | / ( M ulp ) ',
1316 $ / '32 = | I - VT VT**T | / ( N ulp ) ',
1317 $ / ' ZGESVDX(V,V,V) ',
1318 $ / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp )',
1319 $ / '34 = | I - U**T U | / ( M ulp ) ',
1320 $ / '35 = | I - VT VT**T | / ( N ulp ) ',
1321 $ ' ZGESVDQ(H,N,N,A,A',
1322 $ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1323 $ / '37 = | I - U**T U | / ( M ulp ) ',
1324 $ / '38 = | I - VT VT**T | / ( N ulp ) ',
1325 $ / '39 = 0 if S contains min(M,N) nonnegative values in',
1326 $ ' decreasing order, else 1/ulp',
1327 $ / / )
1328 9997 FORMAT( ' M=', i5, ', N=', i5, ', type ', i1, ', IWS=', i1,
1329 $ ', seed=', 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1330 9996 FORMAT( ' ZDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1331 $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1332 $ i5, ')' )
1333 9995 FORMAT( ' ZDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1334 $ i6, ', N=', i6, ', JTYPE=', i6, ', LSWORK=', i6, / 9x,
1335 $ 'ISEED=(', 3( i5, ',' ), i5, ')' )
1336
1337 RETURN
1338
1339
1340
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
subroutine xerbla(srname, info)
double precision function dlarnd(idist, iseed)
DLARND
subroutine zgejsv(joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, cwork, lwork, rwork, lrwork, iwork, info)
ZGEJSV
subroutine zgesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, iwork, info)
ZGESDD
subroutine zgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, info)
ZGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine zgesvdq(joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, cwork, lcwork, rwork, lrwork, info)
ZGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE...
subroutine zgesvdx(jobu, jobvt, range, m, n, a, lda, vl, vu, il, iu, ns, s, u, ldu, vt, ldvt, work, lwork, rwork, iwork, info)
ZGESVDX computes the singular value decomposition (SVD) for GE matrices
subroutine zgesvj(joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, cwork, lwork, rwork, lrwork, info)
ZGESVJ
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine zbdt01(m, n, kd, a, lda, q, ldq, d, e, pt, ldpt, work, rwork, resid)
ZBDT01
subroutine zbdt05(m, n, a, lda, s, ns, u, ldu, vt, ldvt, work, resid)
ZBDT05
subroutine zlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
ZLATMS
subroutine zunt01(rowcol, m, n, u, ldu, work, lwork, rwork, resid)
ZUNT01
subroutine zunt03(rc, mu, mv, n, k, u, ldu, v, ldv, work, lwork, rwork, result, info)
ZUNT03