401
402
403
404
405
406 IMPLICIT NONE
407
408
409 INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES,
410 $ NTYPES
411 REAL THRESH
412
413
414 LOGICAL DOTYPE( * )
415 INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
416 REAL E( * ), RWORK( * ), S( * ), SSAV( * )
417 COMPLEX A( LDA, * ), ASAV( LDA, * ), U( LDU, * ),
418 $ USAV( LDU, * ), VT( LDVT, * ),
419 $ VTSAV( LDVT, * ), WORK( * )
420
421
422
423
424
425 REAL ZERO, ONE, TWO, HALF
426 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
427 $ half = 0.5e0 )
428 COMPLEX CZERO, CONE
429 parameter( czero = ( 0.0e+0, 0.0e+0 ),
430 $ cone = ( 1.0e+0, 0.0e+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 REAL 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 REAL RESULT( 39 )
452
453
454 REAL SLAMCH, SLARND
456
457
461
462
463 INTRINSIC abs, real, 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(
'CDRVBD', -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 310 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 300 jtype = 1, mtypes
560 IF( .NOT.dotype( jtype ) )
561 $ GO TO 300
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 claset(
'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 claset(
'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 clatms( m, n,
'U', iseed,
'N', s, 4, real( 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 clacpy(
'F', m, n, a, lda, asav, lda )
613
614
615
616 DO 290 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 clacpy(
'F', m, n, asav, lda, a, lda )
635 srnamt = 'CGESVD'
636 CALL cgesvd(
'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 cbdt01( 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 cunt01(
'Columns', mnmin, m, usav, ldu, work,
651 $ lwork, rwork, result( 2 ) )
652 CALL cunt01(
'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 clacpy(
'F', m, n, asav, lda, a, lda )
679 srnamt = 'CGESVD'
680 CALL cgesvd( 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 cunt03(
'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 cunt03(
'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 cunt03(
'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 cunt03(
'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 cunt03(
'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 cunt03(
'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( real( mnmin )*ulp*s( 1 ),
727 $
slamch(
'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 clacpy(
'F', m, n, asav, lda, a, lda )
751 srnamt = 'CGESDD'
752 CALL cgesdd(
'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 cbdt01( 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 cunt01(
'Columns', mnmin, m, usav, ldu, work,
767 $ lwork, rwork, result( 9 ) )
768 CALL cunt01(
'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 clacpy(
'F', m, n, asav, lda, a, lda )
791 srnamt = 'CGESDD'
792 CALL cgesdd( 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 cunt03(
'C', m, mnmin, m, mnmin, usav,
802 $ ldu, a, lda, work, lwork, rwork,
803 $ dif, iinfo )
804 ELSE
805 CALL cunt03(
'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 cunt03(
'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 cunt03(
'R', n, mnmin, n, mnmin, vtsav,
824 $ ldvt, vt, ldvt, work, lwork,
825 $ rwork, dif, iinfo )
826 ELSE
827 CALL cunt03(
'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 cunt03(
'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( real( mnmin )*ulp*s( 1 ),
843 $
slamch(
'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
858 result( 36 ) = zero
859 result( 37 ) = zero
860 result( 38 ) = zero
861 result( 39 ) = zero
862
863 IF( m.GE.n ) THEN
864 iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
865 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
866 lswork = min( lswork, lwork )
867 lswork = max( lswork, 1 )
868 IF( iwspc.EQ.4 )
869 $ lswork = lwork
870
871 CALL clacpy(
'F', m, n, asav, lda, a, lda )
872 srnamt = 'CGESVDQ'
873
874 lrwork = max(2, m, 5*n)
875 liwork = max( n, 1 )
876 CALL cgesvdq(
'H',
'N',
'N',
'A',
'A',
877 $ m, n, a, lda, ssav, usav, ldu,
878 $ vtsav, ldvt, numrank, iwork, liwork,
879 $ work, lwork, rwork, lrwork, iinfo )
880
881 IF( iinfo.NE.0 ) THEN
882 WRITE( nounit, fmt = 9995 )'CGESVDQ', iinfo, m, n,
883 $ jtype, lswork, ioldsd
884 info = abs( iinfo )
885 RETURN
886 END IF
887
888
889
890 CALL cbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
891 $ vtsav, ldvt, work, rwork, result( 36 ) )
892 IF( m.NE.0 .AND. n.NE.0 ) THEN
893 CALL cunt01(
'Columns', m, m, usav, ldu, work,
894 $ lwork, rwork, result( 37 ) )
895 CALL cunt01(
'Rows', n, n, vtsav, ldvt, work,
896 $ lwork, rwork, result( 38 ) )
897 END IF
898 result( 39 ) = zero
899 DO 199 i = 1, mnmin - 1
900 IF( ssav( i ).LT.ssav( i+1 ) )
901 $ result( 39 ) = ulpinv
902 IF( ssav( i ).LT.zero )
903 $ result( 39 ) = ulpinv
904 199 CONTINUE
905 IF( mnmin.GE.1 ) THEN
906 IF( ssav( mnmin ).LT.zero )
907 $ result( 39 ) = ulpinv
908 END IF
909 END IF
910
911
912
913
914 result( 15 ) = zero
915 result( 16 ) = zero
916 result( 17 ) = zero
917 result( 18 ) = zero
918
919 IF( m.GE.n ) THEN
920 iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
921 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
922 lswork = min( lswork, lwork )
923 lswork = max( lswork, 1 )
924 lrwork = max(6,n)
925 IF( iwspc.EQ.4 )
926 $ lswork = lwork
927
928 CALL clacpy(
'F', m, n, asav, lda, usav, lda )
929 srnamt = 'CGESVJ'
930 CALL cgesvj(
'G',
'U',
'V', m, n, usav, lda, ssav,
931 & 0, a, ldvt, work, lwork, rwork,
932 & lrwork, iinfo )
933
934
935
936 DO j=1,n
937 DO i=1,n
938 vtsav(j,i) = conjg(a(i,j))
939 END DO
940 END DO
941
942 IF( iinfo.NE.0 ) THEN
943 WRITE( nounit, fmt = 9995 )'GESVJ', iinfo, m, n,
944 $ jtype, lswork, ioldsd
945 info = abs( iinfo )
946 RETURN
947 END IF
948
949
950
951 CALL cbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
952 $ vtsav, ldvt, work, rwork, result( 15 ) )
953 IF( m.NE.0 .AND. n.NE.0 ) THEN
954 CALL cunt01(
'Columns', m, m, usav, ldu, work,
955 $ lwork, rwork, result( 16 ) )
956 CALL cunt01(
'Rows', n, n, vtsav, ldvt, work,
957 $ lwork, rwork, result( 17 ) )
958 END IF
959 result( 18 ) = zero
960 DO 131 i = 1, mnmin - 1
961 IF( ssav( i ).LT.ssav( i+1 ) )
962 $ result( 18 ) = ulpinv
963 IF( ssav( i ).LT.zero )
964 $ result( 18 ) = ulpinv
965 131 CONTINUE
966 IF( mnmin.GE.1 ) THEN
967 IF( ssav( mnmin ).LT.zero )
968 $ result( 18 ) = ulpinv
969 END IF
970 END IF
971
972
973
974
975 result( 19 ) = zero
976 result( 20 ) = zero
977 result( 21 ) = zero
978 result( 22 ) = zero
979 IF( m.GE.n ) THEN
980 iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
981 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
982 lswork = min( lswork, lwork )
983 lswork = max( lswork, 1 )
984 IF( iwspc.EQ.4 )
985 $ lswork = lwork
986 lrwork = max( 7, n + 2*m)
987
988 CALL clacpy(
'F', m, n, asav, lda, vtsav, lda )
989 srnamt = 'CGEJSV'
990 CALL cgejsv(
'G',
'U',
'V',
'R',
'N',
'N',
991 & m, n, vtsav, lda, ssav, usav, ldu, a, ldvt,
992 & work, lwork, rwork,
993 & lrwork, iwork, iinfo )
994
995
996
997 DO 133 j=1,n
998 DO 132 i=1,n
999 vtsav(j,i) = conjg(a(i,j))
1000 132 END DO
1001 133 END DO
1002
1003 IF( iinfo.NE.0 ) THEN
1004 WRITE( nounit, fmt = 9995 )'GEJSV', iinfo, m, n,
1005 $ jtype, lswork, ioldsd
1006 info = abs( iinfo )
1007 RETURN
1008 END IF
1009
1010
1011
1012 CALL cbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
1013 $ vtsav, ldvt, work, rwork, result( 19 ) )
1014 IF( m.NE.0 .AND. n.NE.0 ) THEN
1015 CALL cunt01(
'Columns', m, m, usav, ldu, work,
1016 $ lwork, rwork, result( 20 ) )
1017 CALL cunt01(
'Rows', n, n, vtsav, ldvt, work,
1018 $ lwork, rwork, result( 21 ) )
1019 END IF
1020 result( 22 ) = zero
1021 DO 134 i = 1, mnmin - 1
1022 IF( ssav( i ).LT.ssav( i+1 ) )
1023 $ result( 22 ) = ulpinv
1024 IF( ssav( i ).LT.zero )
1025 $ result( 22 ) = ulpinv
1026 134 CONTINUE
1027 IF( mnmin.GE.1 ) THEN
1028 IF( ssav( mnmin ).LT.zero )
1029 $ result( 22 ) = ulpinv
1030 END IF
1031 END IF
1032
1033
1034
1035
1036
1037 CALL clacpy(
'F', m, n, asav, lda, a, lda )
1038 srnamt = 'CGESVDX'
1039 CALL cgesvdx(
'V',
'V',
'A', m, n, a, lda,
1040 $ vl, vu, il, iu, ns, ssav, usav, ldu,
1041 $ vtsav, ldvt, work, lwork, rwork,
1042 $ iwork, iinfo )
1043 IF( iinfo.NE.0 ) THEN
1044 WRITE( nounit, fmt = 9995 )'GESVDX', iinfo, m, n,
1045 $ jtype, lswork, ioldsd
1046 info = abs( iinfo )
1047 RETURN
1048 END IF
1049
1050
1051
1052 result( 23 ) = zero
1053 result( 24 ) = zero
1054 result( 25 ) = zero
1055 CALL cbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
1056 $ vtsav, ldvt, work, rwork, result( 23 ) )
1057 IF( m.NE.0 .AND. n.NE.0 ) THEN
1058 CALL cunt01(
'Columns', mnmin, m, usav, ldu, work,
1059 $ lwork, rwork, result( 24 ) )
1060 CALL cunt01(
'Rows', mnmin, n, vtsav, ldvt, work,
1061 $ lwork, rwork, result( 25 ) )
1062 END IF
1063 result( 26 ) = zero
1064 DO 140 i = 1, mnmin - 1
1065 IF( ssav( i ).LT.ssav( i+1 ) )
1066 $ result( 26 ) = ulpinv
1067 IF( ssav( i ).LT.zero )
1068 $ result( 26 ) = ulpinv
1069 140 CONTINUE
1070 IF( mnmin.GE.1 ) THEN
1071 IF( ssav( mnmin ).LT.zero )
1072 $ result( 26 ) = ulpinv
1073 END IF
1074
1075
1076
1077 result( 27 ) = zero
1078 result( 28 ) = zero
1079 result( 29 ) = zero
1080 DO 170 iju = 0, 1
1081 DO 160 ijvt = 0, 1
1082 IF( ( iju.EQ.0 .AND. ijvt.EQ.0 ) .OR.
1083 $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) ) GO TO 160
1084 jobu = cjobv( iju+1 )
1085 jobvt = cjobv( ijvt+1 )
1086 range = cjobr( 1 )
1087 CALL clacpy(
'F', m, n, asav, lda, a, lda )
1088 srnamt = 'CGESVDX'
1089 CALL cgesvdx( jobu, jobvt,
'A', m, n, a, lda,
1090 $ vl, vu, il, iu, ns, ssav, u, ldu,
1091 $ vt, ldvt, work, lwork, rwork,
1092 $ iwork, iinfo )
1093
1094
1095
1096 dif = zero
1097 IF( m.GT.0 .AND. n.GT.0 ) THEN
1098 IF( iju.EQ.1 ) THEN
1099 CALL cunt03(
'C', m, mnmin, m, mnmin, usav,
1100 $ ldu, u, ldu, work, lwork, rwork,
1101 $ dif, iinfo )
1102 END IF
1103 END IF
1104 result( 27 ) = max( result( 27 ), dif )
1105
1106
1107
1108 dif = zero
1109 IF( m.GT.0 .AND. n.GT.0 ) THEN
1110 IF( ijvt.EQ.1 ) THEN
1111 CALL cunt03(
'R', n, mnmin, n, mnmin, vtsav,
1112 $ ldvt, vt, ldvt, work, lwork,
1113 $ rwork, dif, iinfo )
1114 END IF
1115 END IF
1116 result( 28 ) = max( result( 28 ), dif )
1117
1118
1119
1120 dif = zero
1121 div = max( real( mnmin )*ulp*s( 1 ),
1122 $
slamch(
'Safe minimum' ) )
1123 DO 150 i = 1, mnmin - 1
1124 IF( ssav( i ).LT.ssav( i+1 ) )
1125 $ dif = ulpinv
1126 IF( ssav( i ).LT.zero )
1127 $ dif = ulpinv
1128 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
1129 150 CONTINUE
1130 result( 29) = max( result( 29 ), dif )
1131 160 CONTINUE
1132 170 CONTINUE
1133
1134
1135
1136 DO 180 i = 1, 4
1137 iseed2( i ) = iseed( i )
1138 180 CONTINUE
1139 IF( mnmin.LE.1 ) THEN
1140 il = 1
1141 iu = max( 1, mnmin )
1142 ELSE
1143 il = 1 + int( ( mnmin-1 )*
slarnd( 1, iseed2 ) )
1144 iu = 1 + int( ( mnmin-1 )*
slarnd( 1, iseed2 ) )
1145 IF( iu.LT.il ) THEN
1146 itemp = iu
1147 iu = il
1148 il = itemp
1149 END IF
1150 END IF
1151 CALL clacpy(
'F', m, n, asav, lda, a, lda )
1152 srnamt = 'CGESVDX'
1153 CALL cgesvdx(
'V',
'V',
'I', m, n, a, lda,
1154 $ vl, vu, il, iu, nsi, s, u, ldu,
1155 $ vt, ldvt, work, lwork, rwork,
1156 $ iwork, iinfo )
1157 IF( iinfo.NE.0 ) THEN
1158 WRITE( nounit, fmt = 9995 )'GESVDX', iinfo, m, n,
1159 $ jtype, lswork, ioldsd
1160 info = abs( iinfo )
1161 RETURN
1162 END IF
1163
1164 result( 30 ) = zero
1165 result( 31 ) = zero
1166 result( 32 ) = zero
1167 CALL cbdt05( m, n, asav, lda, s, nsi, u, ldu,
1168 $ vt, ldvt, work, result( 30 ) )
1169 IF( m.NE.0 .AND. n.NE.0 ) THEN
1170 CALL cunt01(
'Columns', m, nsi, u, ldu, work,
1171 $ lwork, rwork, result( 31 ) )
1172 CALL cunt01(
'Rows', nsi, n, vt, ldvt, work,
1173 $ lwork, rwork, result( 32 ) )
1174 END IF
1175
1176
1177
1178 IF( mnmin.GT.0 .AND. nsi.GT.1 ) THEN
1179 IF( il.NE.1 ) THEN
1180 vu = ssav( il ) +
1181 $ max( half*abs( ssav( il )-ssav( il-1 ) ),
1182 $ ulp*anorm, two*rtunfl )
1183 ELSE
1184 vu = ssav( 1 ) +
1185 $ max( half*abs( ssav( ns )-ssav( 1 ) ),
1186 $ ulp*anorm, two*rtunfl )
1187 END IF
1188 IF( iu.NE.ns ) THEN
1189 vl = ssav( iu ) - max( ulp*anorm, two*rtunfl,
1190 $ half*abs( ssav( iu+1 )-ssav( iu ) ) )
1191 ELSE
1192 vl = ssav( ns ) - max( ulp*anorm, two*rtunfl,
1193 $ half*abs( ssav( ns )-ssav( 1 ) ) )
1194 END IF
1195 vl = max( vl,zero )
1196 vu = max( vu,zero )
1197 IF( vl.GE.vu ) vu = max( vu*2, vu+vl+half )
1198 ELSE
1199 vl = zero
1200 vu = one
1201 END IF
1202 CALL clacpy(
'F', m, n, asav, lda, a, lda )
1203 srnamt = 'CGESVDX'
1204 CALL cgesvdx(
'V',
'V',
'V', m, n, a, lda,
1205 $ vl, vu, il, iu, nsv, s, u, ldu,
1206 $ vt, ldvt, work, lwork, rwork,
1207 $ iwork, iinfo )
1208 IF( iinfo.NE.0 ) THEN
1209 WRITE( nounit, fmt = 9995 )'GESVDX', iinfo, m, n,
1210 $ jtype, lswork, ioldsd
1211 info = abs( iinfo )
1212 RETURN
1213 END IF
1214
1215 result( 33 ) = zero
1216 result( 34 ) = zero
1217 result( 35 ) = zero
1218 CALL cbdt05( m, n, asav, lda, s, nsv, u, ldu,
1219 $ vt, ldvt, work, result( 33 ) )
1220 IF( m.NE.0 .AND. n.NE.0 ) THEN
1221 CALL cunt01(
'Columns', m, nsv, u, ldu, work,
1222 $ lwork, rwork, result( 34 ) )
1223 CALL cunt01(
'Rows', nsv, n, vt, ldvt, work,
1224 $ lwork, rwork, result( 35 ) )
1225 END IF
1226
1227
1228
1229 ntest = 0
1230 nfail = 0
1231 DO 190 j = 1, 39
1232 IF( result( j ).GE.zero )
1233 $ ntest = ntest + 1
1234 IF( result( j ).GE.thresh )
1235 $ nfail = nfail + 1
1236 190 CONTINUE
1237
1238 IF( nfail.GT.0 )
1239 $ ntestf = ntestf + 1
1240 IF( ntestf.EQ.1 ) THEN
1241 WRITE( nounit, fmt = 9999 )
1242 WRITE( nounit, fmt = 9998 )thresh
1243 ntestf = 2
1244 END IF
1245
1246 DO 200 j = 1, 39
1247 IF( result( j ).GE.thresh ) THEN
1248 WRITE( nounit, fmt = 9997 )m, n, jtype, iwspc,
1249 $ ioldsd, j, result( j )
1250 END IF
1251 200 CONTINUE
1252
1253 nerrs = nerrs + nfail
1254 ntestt = ntestt + ntest
1255
1256 290 CONTINUE
1257
1258 300 CONTINUE
1259 310 CONTINUE
1260
1261
1262
1263 CALL alasvm(
'CBD', nounit, nerrs, ntestt, 0 )
1264
1265 9999 FORMAT( ' SVD -- Complex Singular Value Decomposition Driver ',
1266 $ / ' Matrix types (see CDRVBD for details):',
1267 $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
1268 $ / ' 3 = Evenly spaced singular values near 1',
1269 $ / ' 4 = Evenly spaced singular values near underflow',
1270 $ / ' 5 = Evenly spaced singular values near overflow',
1271 $ / / ' Tests performed: ( A is dense, U and V are unitary,',
1272 $ / 19x, ' S is an array, and Upartial, VTpartial, and',
1273 $ / 19x, ' Spartial are partially computed U, VT and S),', / )
1274 9998 FORMAT( ' Tests performed with Test Threshold = ', f8.2,
1275 $ / ' CGESVD: ', /
1276 $ ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1277 $ / ' 2 = | I - U**T U | / ( M ulp ) ',
1278 $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
1279 $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
1280 $ ' decreasing order, else 1/ulp',
1281 $ / ' 5 = | U - Upartial | / ( M ulp )',
1282 $ / ' 6 = | VT - VTpartial | / ( N ulp )',
1283 $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
1284 $ / ' CGESDD: ', /
1285 $ ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1286 $ / ' 9 = | I - U**T U | / ( M ulp ) ',
1287 $ / '10 = | I - VT VT**T | / ( N ulp ) ',
1288 $ / '11 = 0 if S contains min(M,N) nonnegative values in',
1289 $ ' decreasing order, else 1/ulp',
1290 $ / '12 = | U - Upartial | / ( M ulp )',
1291 $ / '13 = | VT - VTpartial | / ( N ulp )',
1292 $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )',
1293 $ / ' CGESVJ: ', /
1294 $ / '15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1295 $ / '16 = | I - U**T U | / ( M ulp ) ',
1296 $ / '17 = | I - VT VT**T | / ( N ulp ) ',
1297 $ / '18 = 0 if S contains min(M,N) nonnegative values in',
1298 $ ' decreasing order, else 1/ulp',
1299 $ / ' CGESJV: ', /
1300 $ / '19 = | A - U diag(S) VT | / ( |A| max(M,N) ulp )',
1301 $ / '20 = | I - U**T U | / ( M ulp ) ',
1302 $ / '21 = | I - VT VT**T | / ( N ulp ) ',
1303 $ / '22 = 0 if S contains min(M,N) nonnegative values in',
1304 $ ' decreasing order, else 1/ulp',
1305 $ / ' CGESVDX(V,V,A): ', /
1306 $ '23 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1307 $ / '24 = | I - U**T U | / ( M ulp ) ',
1308 $ / '25 = | I - VT VT**T | / ( N ulp ) ',
1309 $ / '26 = 0 if S contains min(M,N) nonnegative values in',
1310 $ ' decreasing order, else 1/ulp',
1311 $ / '27 = | U - Upartial | / ( M ulp )',
1312 $ / '28 = | VT - VTpartial | / ( N ulp )',
1313 $ / '29 = | S - Spartial | / ( min(M,N) ulp |S| )',
1314 $ / ' CGESVDX(V,V,I): ',
1315 $ / '30 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp )',
1316 $ / '31 = | I - U**T U | / ( M ulp ) ',
1317 $ / '32 = | I - VT VT**T | / ( N ulp ) ',
1318 $ / ' CGESVDX(V,V,V) ',
1319 $ / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp )',
1320 $ / '34 = | I - U**T U | / ( M ulp ) ',
1321 $ / '35 = | I - VT VT**T | / ( N ulp ) ',
1322 $ ' CGESVDQ(H,N,N,A,A',
1323 $ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1324 $ / '37 = | I - U**T U | / ( M ulp ) ',
1325 $ / '38 = | I - VT VT**T | / ( N ulp ) ',
1326 $ / '39 = 0 if S contains min(M,N) nonnegative values in',
1327 $ ' decreasing order, else 1/ulp',
1328 $ / / )
1329 9997 FORMAT( ' M=', i5, ', N=', i5, ', type ', i1, ', IWS=', i1,
1330 $ ', seed=', 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1331 9996 FORMAT( ' CDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1332 $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1333 $ i5, ')' )
1334 9995 FORMAT( ' CDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1335 $ i6, ', N=', i6, ', JTYPE=', i6, ', LSWORK=', i6, / 9x,
1336 $ 'ISEED=(', 3( i5, ',' ), i5, ')' )
1337
1338 RETURN
1339
1340
1341
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
subroutine cbdt01(m, n, kd, a, lda, q, ldq, d, e, pt, ldpt, work, rwork, resid)
CBDT01
subroutine cbdt05(m, n, a, lda, s, ns, u, ldu, vt, ldvt, work, resid)
CBDT05
subroutine xerbla(srname, info)
subroutine clatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
CLATMS
subroutine cunt01(rowcol, m, n, u, ldu, work, lwork, rwork, resid)
CUNT01
subroutine cunt03(rc, mu, mv, n, k, u, ldu, v, ldv, work, lwork, rwork, result, info)
CUNT03
subroutine cgejsv(joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, cwork, lwork, rwork, lrwork, iwork, info)
CGEJSV
subroutine cgesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, iwork, info)
CGESDD
subroutine cgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, info)
CGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine cgesvdq(joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, cwork, lcwork, rwork, lrwork, info)
CGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE...
subroutine cgesvdx(jobu, jobvt, range, m, n, a, lda, vl, vu, il, iu, ns, s, u, ldu, vt, ldvt, work, lwork, rwork, iwork, info)
CGESVDX computes the singular value decomposition (SVD) for GE matrices
subroutine cgesvj(joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, cwork, lwork, rwork, lrwork, info)
CGESVJ
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
real function slarnd(idist, iseed)
SLARND