366
367 IMPLICIT NONE
368
369
370
371
372
373
374 INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUT, NSIZES,
375 $ NTYPES
376 DOUBLE PRECISION THRESH
377
378
379 LOGICAL DOTYPE( * )
380 INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
381 DOUBLE PRECISION A( LDA, * ), ASAV( LDA, * ), E( * ), S( * ),
382 $ SSAV( * ), U( LDU, * ), USAV( LDU, * ),
383 $ VT( LDVT, * ), VTSAV( LDVT, * ), WORK( * )
384
385
386
387
388
389 DOUBLE PRECISION ZERO, ONE, TWO, HALF
390 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
391 $ half = 0.5d0 )
392 INTEGER MAXTYP
393 parameter( maxtyp = 5 )
394
395
396 LOGICAL BADMM, BADNN
397 CHARACTER JOBQ, JOBU, JOBVT, RANGE
398 CHARACTER*3 PATH
399 INTEGER I, IINFO, IJQ, IJU, IJVT, IL,IU, IWS, IWTMP,
400 $ ITEMP, J, JSIZE, JTYPE, LSWORK, M, MINWRK,
401 $ MMAX, MNMAX, MNMIN, MTYPES, N, NFAIL,
402 $ NMAX, NS, NSI, NSV, NTEST
403 DOUBLE PRECISION ANORM, DIF, DIV, OVFL, RTUNFL, ULP,
404 $ ULPINV, UNFL, VL, VU
405
406
407 INTEGER LIWORK, LRWORK, NUMRANK
408
409
410 DOUBLE PRECISION RWORK( 2 )
411
412
413 CHARACTER CJOB( 4 ), CJOBR( 3 ), CJOBV( 2 )
414 INTEGER IOLDSD( 4 ), ISEED2( 4 )
415 DOUBLE PRECISION RESULT( 39 )
416
417
418 DOUBLE PRECISION DLAMCH, DLARND
420
421
425
426
427 INTRINSIC abs, dble, int, max, min
428
429
430 LOGICAL LERR, OK
431 CHARACTER*32 SRNAMT
432 INTEGER INFOT, NUNIT
433
434
435 COMMON / infoc / infot, nunit, ok, lerr
436 COMMON / srnamc / srnamt
437
438
439 DATA cjob / 'N', 'O', 'S', 'A' /
440 DATA cjobr / 'A', 'V', 'I' /
441 DATA cjobv / 'N', 'V' /
442
443
444
445
446
447 info = 0
448 badmm = .false.
449 badnn = .false.
450 mmax = 1
451 nmax = 1
452 mnmax = 1
453 minwrk = 1
454 DO 10 j = 1, nsizes
455 mmax = max( mmax, mm( j ) )
456 IF( mm( j ).LT.0 )
457 $ badmm = .true.
458 nmax = max( nmax, nn( j ) )
459 IF( nn( j ).LT.0 )
460 $ badnn = .true.
461 mnmax = max( mnmax, min( mm( j ), nn( j ) ) )
462 minwrk = max( minwrk, max( 3*min( mm( j ),
463 $ nn( j ) )+max( mm( j ), nn( j ) ), 5*min( mm( j ),
464 $ nn( j )-4 ) )+2*min( mm( j ), nn( j ) )**2 )
465 10 CONTINUE
466
467
468
469 IF( nsizes.LT.0 ) THEN
470 info = -1
471 ELSE IF( badmm ) THEN
472 info = -2
473 ELSE IF( badnn ) THEN
474 info = -3
475 ELSE IF( ntypes.LT.0 ) THEN
476 info = -4
477 ELSE IF( lda.LT.max( 1, mmax ) ) THEN
478 info = -10
479 ELSE IF( ldu.LT.max( 1, mmax ) ) THEN
480 info = -12
481 ELSE IF( ldvt.LT.max( 1, nmax ) ) THEN
482 info = -14
483 ELSE IF( minwrk.GT.lwork ) THEN
484 info = -21
485 END IF
486
487 IF( info.NE.0 ) THEN
488 CALL xerbla(
'DDRVBD', -info )
489 RETURN
490 END IF
491
492
493
494 path( 1: 1 ) = 'Double precision'
495 path( 2: 3 ) = 'BD'
496 nfail = 0
497 ntest = 0
498 unfl =
dlamch(
'Safe minimum' )
499 ovfl = one / unfl
501 ulp =
dlamch(
'Precision' )
502 rtunfl = sqrt( unfl )
503 ulpinv = one / ulp
504 infot = 0
505
506
507
508 DO 240 jsize = 1, nsizes
509 m = mm( jsize )
510 n = nn( jsize )
511 mnmin = min( m, n )
512
513 IF( nsizes.NE.1 ) THEN
514 mtypes = min( maxtyp, ntypes )
515 ELSE
516 mtypes = min( maxtyp+1, ntypes )
517 END IF
518
519 DO 230 jtype = 1, mtypes
520 IF( .NOT.dotype( jtype ) )
521 $ GO TO 230
522
523 DO 20 j = 1, 4
524 ioldsd( j ) = iseed( j )
525 20 CONTINUE
526
527
528
529 IF( mtypes.GT.maxtyp )
530 $ GO TO 30
531
532 IF( jtype.EQ.1 ) THEN
533
534
535
536 CALL dlaset(
'Full', m, n, zero, zero, a, lda )
537
538 ELSE IF( jtype.EQ.2 ) THEN
539
540
541
542 CALL dlaset(
'Full', m, n, zero, one, a, lda )
543
544 ELSE
545
546
547
548 IF( jtype.EQ.3 )
549 $ anorm = one
550 IF( jtype.EQ.4 )
551 $ anorm = unfl / ulp
552 IF( jtype.EQ.5 )
553 $ anorm = ovfl*ulp
554 CALL dlatms( m, n,
'U', iseed,
'N', s, 4, dble( mnmin ),
555 $ anorm, m-1, n-1, 'N', a, lda, work, iinfo )
556 IF( iinfo.NE.0 ) THEN
557 WRITE( nout, fmt = 9996 )'Generator', iinfo, m, n,
558 $ jtype, ioldsd
559 info = abs( iinfo )
560 RETURN
561 END IF
562 END IF
563
564 30 CONTINUE
565 CALL dlacpy(
'F', m, n, a, lda, asav, lda )
566
567
568
569 DO 220 iws = 1, 4
570
571 DO 40 j = 1, 32
572 result( j ) = -one
573 40 CONTINUE
574
575
576
577 iwtmp = max( 3*min( m, n )+max( m, n ), 5*min( m, n ) )
578 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
579 lswork = min( lswork, lwork )
580 lswork = max( lswork, 1 )
581 IF( iws.EQ.4 )
582 $ lswork = lwork
583
584 IF( iws.GT.1 )
585 $
CALL dlacpy(
'F', m, n, asav, lda, a, lda )
586 srnamt = 'DGESVD'
587 CALL dgesvd(
'A',
'A', m, n, a, lda, ssav, usav, ldu,
588 $ vtsav, ldvt, work, lswork, iinfo )
589 IF( iinfo.NE.0 ) THEN
590 WRITE( nout, fmt = 9995 )'GESVD', iinfo, m, n, jtype,
591 $ lswork, ioldsd
592 info = abs( iinfo )
593 RETURN
594 END IF
595
596
597
598 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
599 $ vtsav, ldvt, work, result( 1 ) )
600 IF( m.NE.0 .AND. n.NE.0 ) THEN
601 CALL dort01(
'Columns', m, m, usav, ldu, work, lwork,
602 $ result( 2 ) )
603 CALL dort01(
'Rows', n, n, vtsav, ldvt, work, lwork,
604 $ result( 3 ) )
605 END IF
606 result( 4 ) = zero
607 DO 50 i = 1, mnmin - 1
608 IF( ssav( i ).LT.ssav( i+1 ) )
609 $ result( 4 ) = ulpinv
610 IF( ssav( i ).LT.zero )
611 $ result( 4 ) = ulpinv
612 50 CONTINUE
613 IF( mnmin.GE.1 ) THEN
614 IF( ssav( mnmin ).LT.zero )
615 $ result( 4 ) = ulpinv
616 END IF
617
618
619
620 result( 5 ) = zero
621 result( 6 ) = zero
622 result( 7 ) = zero
623 DO 80 iju = 0, 3
624 DO 70 ijvt = 0, 3
625 IF( ( iju.EQ.3 .AND. ijvt.EQ.3 ) .OR.
626 $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )GO TO 70
627 jobu = cjob( iju+1 )
628 jobvt = cjob( ijvt+1 )
629 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
630 srnamt = 'DGESVD'
631 CALL dgesvd( jobu, jobvt, m, n, a, lda, s, u, ldu,
632 $ vt, ldvt, work, lswork, iinfo )
633
634
635
636 dif = zero
637 IF( m.GT.0 .AND. n.GT.0 ) THEN
638 IF( iju.EQ.1 ) THEN
639 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
640 $ ldu, a, lda, work, lwork, dif,
641 $ iinfo )
642 ELSE IF( iju.EQ.2 ) THEN
643 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
644 $ ldu, u, ldu, work, lwork, dif,
645 $ iinfo )
646 ELSE IF( iju.EQ.3 ) THEN
647 CALL dort03(
'C', m, m, m, mnmin, usav, ldu,
648 $ u, ldu, work, lwork, dif,
649 $ iinfo )
650 END IF
651 END IF
652 result( 5 ) = max( result( 5 ), dif )
653
654
655
656 dif = zero
657 IF( m.GT.0 .AND. n.GT.0 ) THEN
658 IF( ijvt.EQ.1 ) THEN
659 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
660 $ ldvt, a, lda, work, lwork, dif,
661 $ iinfo )
662 ELSE IF( ijvt.EQ.2 ) THEN
663 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
664 $ ldvt, vt, ldvt, work, lwork,
665 $ dif, iinfo )
666 ELSE IF( ijvt.EQ.3 ) THEN
667 CALL dort03(
'R', n, n, n, mnmin, vtsav,
668 $ ldvt, vt, ldvt, work, lwork,
669 $ dif, iinfo )
670 END IF
671 END IF
672 result( 6 ) = max( result( 6 ), dif )
673
674
675
676 dif = zero
677 div = max( mnmin*ulp*s( 1 ), unfl )
678 DO 60 i = 1, mnmin - 1
679 IF( ssav( i ).LT.ssav( i+1 ) )
680 $ dif = ulpinv
681 IF( ssav( i ).LT.zero )
682 $ dif = ulpinv
683 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
684 60 CONTINUE
685 result( 7 ) = max( result( 7 ), dif )
686 70 CONTINUE
687 80 CONTINUE
688
689
690
691 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
692 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
693 lswork = min( lswork, lwork )
694 lswork = max( lswork, 1 )
695 IF( iws.EQ.4 )
696 $ lswork = lwork
697
698 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
699 srnamt = 'DGESDD'
700 CALL dgesdd(
'A', m, n, a, lda, ssav, usav, ldu, vtsav,
701 $ ldvt, work, lswork, iwork, iinfo )
702 IF( iinfo.NE.0 ) THEN
703 WRITE( nout, fmt = 9995 )'GESDD', iinfo, m, n, jtype,
704 $ lswork, ioldsd
705 info = abs( iinfo )
706 RETURN
707 END IF
708
709
710
711 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
712 $ vtsav, ldvt, work, result( 8 ) )
713 IF( m.NE.0 .AND. n.NE.0 ) THEN
714 CALL dort01(
'Columns', m, m, usav, ldu, work, lwork,
715 $ result( 9 ) )
716 CALL dort01(
'Rows', n, n, vtsav, ldvt, work, lwork,
717 $ result( 10 ) )
718 END IF
719 result( 11 ) = zero
720 DO 90 i = 1, mnmin - 1
721 IF( ssav( i ).LT.ssav( i+1 ) )
722 $ result( 11 ) = ulpinv
723 IF( ssav( i ).LT.zero )
724 $ result( 11 ) = ulpinv
725 90 CONTINUE
726 IF( mnmin.GE.1 ) THEN
727 IF( ssav( mnmin ).LT.zero )
728 $ result( 11 ) = ulpinv
729 END IF
730
731
732
733 result( 12 ) = zero
734 result( 13 ) = zero
735 result( 14 ) = zero
736 DO 110 ijq = 0, 2
737 jobq = cjob( ijq+1 )
738 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
739 srnamt = 'DGESDD'
740 CALL dgesdd( jobq, m, n, a, lda, s, u, ldu, vt, ldvt,
741 $ work, lswork, iwork, iinfo )
742
743
744
745 dif = zero
746 IF( m.GT.0 .AND. n.GT.0 ) THEN
747 IF( ijq.EQ.1 ) THEN
748 IF( m.GE.n ) THEN
749 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
750 $ ldu, a, lda, work, lwork, dif,
751 $ info )
752 ELSE
753 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
754 $ ldu, u, ldu, work, lwork, dif,
755 $ info )
756 END IF
757 ELSE IF( ijq.EQ.2 ) THEN
758 CALL dort03(
'C', m, mnmin, m, mnmin, usav, ldu,
759 $ u, ldu, work, lwork, dif, info )
760 END IF
761 END IF
762 result( 12 ) = max( result( 12 ), dif )
763
764
765
766 dif = zero
767 IF( m.GT.0 .AND. n.GT.0 ) THEN
768 IF( ijq.EQ.1 ) THEN
769 IF( m.GE.n ) THEN
770 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
771 $ ldvt, vt, ldvt, work, lwork,
772 $ dif, info )
773 ELSE
774 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
775 $ ldvt, a, lda, work, lwork, dif,
776 $ info )
777 END IF
778 ELSE IF( ijq.EQ.2 ) THEN
779 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
780 $ ldvt, vt, ldvt, work, lwork, dif,
781 $ info )
782 END IF
783 END IF
784 result( 13 ) = max( result( 13 ), dif )
785
786
787
788 dif = zero
789 div = max( mnmin*ulp*s( 1 ), unfl )
790 DO 100 i = 1, mnmin - 1
791 IF( ssav( i ).LT.ssav( i+1 ) )
792 $ dif = ulpinv
793 IF( ssav( i ).LT.zero )
794 $ dif = ulpinv
795 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
796 100 CONTINUE
797 result( 14 ) = max( result( 14 ), dif )
798 110 CONTINUE
799
800
801
802
803 result( 36 ) = zero
804 result( 37 ) = zero
805 result( 38 ) = zero
806 result( 39 ) = zero
807
808 IF( m.GE.n ) THEN
809 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
810 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
811 lswork = min( lswork, lwork )
812 lswork = max( lswork, 1 )
813 IF( iws.EQ.4 )
814 $ lswork = lwork
815
816 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
817 srnamt = 'DGESVDQ'
818
819 lrwork = 2
820 liwork = max( n, 1 )
821 CALL dgesvdq(
'H',
'N',
'N',
'A',
'A',
822 $ m, n, a, lda, ssav, usav, ldu,
823 $ vtsav, ldvt, numrank, iwork, liwork,
824 $ work, lwork, rwork, lrwork, iinfo )
825
826 IF( iinfo.NE.0 ) THEN
827 WRITE( nout, fmt = 9995 )'DGESVDQ', iinfo, m, n,
828 $ jtype, lswork, ioldsd
829 info = abs( iinfo )
830 RETURN
831 END IF
832
833
834
835 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
836 $ vtsav, ldvt, work, result( 36 ) )
837 IF( m.NE.0 .AND. n.NE.0 ) THEN
838 CALL dort01(
'Columns', m, m, usav, ldu, work,
839 $ lwork, result( 37 ) )
840 CALL dort01(
'Rows', n, n, vtsav, ldvt, work,
841 $ lwork, result( 38 ) )
842 END IF
843 result( 39 ) = zero
844 DO 199 i = 1, mnmin - 1
845 IF( ssav( i ).LT.ssav( i+1 ) )
846 $ result( 39 ) = ulpinv
847 IF( ssav( i ).LT.zero )
848 $ result( 39 ) = ulpinv
849 199 CONTINUE
850 IF( mnmin.GE.1 ) THEN
851 IF( ssav( mnmin ).LT.zero )
852 $ result( 39 ) = ulpinv
853 END IF
854 END IF
855
856
857
858
859 result( 15 ) = zero
860 result( 16 ) = zero
861 result( 17 ) = zero
862 result( 18 ) = zero
863
864 IF( m.GE.n ) THEN
865 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
866 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
867 lswork = min( lswork, lwork )
868 lswork = max( lswork, 1 )
869 IF( iws.EQ.4 )
870 $ lswork = lwork
871
872 CALL dlacpy(
'F', m, n, asav, lda, usav, lda )
873 srnamt = 'DGESVJ'
874 CALL dgesvj(
'G',
'U',
'V', m, n, usav, lda, ssav,
875 & 0, a, ldvt, work, lwork, info )
876
877
878
879 DO j=1,n
880 DO i=1,n
881 vtsav(j,i) = a(i,j)
882 END DO
883 END DO
884
885 IF( iinfo.NE.0 ) THEN
886 WRITE( nout, fmt = 9995 )'GESVJ', iinfo, m, n,
887 $ jtype, lswork, ioldsd
888 info = abs( iinfo )
889 RETURN
890 END IF
891
892
893
894 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
895 $ vtsav, ldvt, work, result( 15 ) )
896 IF( m.NE.0 .AND. n.NE.0 ) THEN
897 CALL dort01(
'Columns', m, m, usav, ldu, work,
898 $ lwork, result( 16 ) )
899 CALL dort01(
'Rows', n, n, vtsav, ldvt, work,
900 $ lwork, result( 17 ) )
901 END IF
902 result( 18 ) = zero
903 DO 120 i = 1, mnmin - 1
904 IF( ssav( i ).LT.ssav( i+1 ) )
905 $ result( 18 ) = ulpinv
906 IF( ssav( i ).LT.zero )
907 $ result( 18 ) = ulpinv
908 120 CONTINUE
909 IF( mnmin.GE.1 ) THEN
910 IF( ssav( mnmin ).LT.zero )
911 $ result( 18 ) = ulpinv
912 END IF
913 END IF
914
915
916
917
918 result( 19 ) = zero
919 result( 20 ) = zero
920 result( 21 ) = zero
921 result( 22 ) = zero
922 IF( m.GE.n ) THEN
923 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
924 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
925 lswork = min( lswork, lwork )
926 lswork = max( lswork, 1 )
927 IF( iws.EQ.4 )
928 $ lswork = lwork
929
930 CALL dlacpy(
'F', m, n, asav, lda, vtsav, lda )
931 srnamt = 'DGEJSV'
932 CALL dgejsv(
'G',
'U',
'V',
'R',
'N',
'N',
933 & m, n, vtsav, lda, ssav, usav, ldu, a, ldvt,
934 & work, lwork, iwork, info )
935
936
937
938 DO 140 j=1,n
939 DO 130 i=1,n
940 vtsav(j,i) = a(i,j)
941 130 END DO
942 140 END DO
943
944 IF( iinfo.NE.0 ) THEN
945 WRITE( nout, fmt = 9995 )'GEJSV', iinfo, m, n,
946 $ jtype, lswork, ioldsd
947 info = abs( iinfo )
948 RETURN
949 END IF
950
951
952
953 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
954 $ vtsav, ldvt, work, result( 19 ) )
955 IF( m.NE.0 .AND. n.NE.0 ) THEN
956 CALL dort01(
'Columns', m, m, usav, ldu, work,
957 $ lwork, result( 20 ) )
958 CALL dort01(
'Rows', n, n, vtsav, ldvt, work,
959 $ lwork, result( 21 ) )
960 END IF
961 result( 22 ) = zero
962 DO 150 i = 1, mnmin - 1
963 IF( ssav( i ).LT.ssav( i+1 ) )
964 $ result( 22 ) = ulpinv
965 IF( ssav( i ).LT.zero )
966 $ result( 22 ) = ulpinv
967 150 CONTINUE
968 IF( mnmin.GE.1 ) THEN
969 IF( ssav( mnmin ).LT.zero )
970 $ result( 22 ) = ulpinv
971 END IF
972 END IF
973
974
975
976 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
977 CALL dgesvdx(
'V',
'V',
'A', m, n, a, lda,
978 $ vl, vu, il, iu, ns, ssav, usav, ldu,
979 $ vtsav, ldvt, work, lwork, iwork,
980 $ iinfo )
981 IF( iinfo.NE.0 ) THEN
982 WRITE( nout, fmt = 9995 )'GESVDX', iinfo, m, n,
983 $ jtype, lswork, ioldsd
984 info = abs( iinfo )
985 RETURN
986 END IF
987
988
989
990 result( 23 ) = zero
991 result( 24 ) = zero
992 result( 25 ) = zero
993 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
994 $ vtsav, ldvt, work, result( 23 ) )
995 IF( m.NE.0 .AND. n.NE.0 ) THEN
996 CALL dort01(
'Columns', m, m, usav, ldu, work, lwork,
997 $ result( 24 ) )
998 CALL dort01(
'Rows', n, n, vtsav, ldvt, work, lwork,
999 $ result( 25 ) )
1000 END IF
1001 result( 26 ) = zero
1002 DO 160 i = 1, mnmin - 1
1003 IF( ssav( i ).LT.ssav( i+1 ) )
1004 $ result( 26 ) = ulpinv
1005 IF( ssav( i ).LT.zero )
1006 $ result( 26 ) = ulpinv
1007 160 CONTINUE
1008 IF( mnmin.GE.1 ) THEN
1009 IF( ssav( mnmin ).LT.zero )
1010 $ result( 26 ) = ulpinv
1011 END IF
1012
1013
1014
1015 result( 27 ) = zero
1016 result( 28 ) = zero
1017 result( 29 ) = zero
1018 DO 180 iju = 0, 1
1019 DO 170 ijvt = 0, 1
1020 IF( ( iju.EQ.0 .AND. ijvt.EQ.0 ) .OR.
1021 $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )GO TO 170
1022 jobu = cjobv( iju+1 )
1023 jobvt = cjobv( ijvt+1 )
1024 range = cjobr( 1 )
1025 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
1026 CALL dgesvdx( jobu, jobvt, range, m, n, a, lda,
1027 $ vl, vu, il, iu, ns, s, u, ldu,
1028 $ vt, ldvt, work, lwork, iwork,
1029 $ iinfo )
1030
1031
1032
1033 dif = zero
1034 IF( m.GT.0 .AND. n.GT.0 ) THEN
1035 IF( iju.EQ.1 ) THEN
1036 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
1037 $ ldu, u, ldu, work, lwork, dif,
1038 $ iinfo )
1039 END IF
1040 END IF
1041 result( 27 ) = max( result( 27 ), dif )
1042
1043
1044
1045 dif = zero
1046 IF( m.GT.0 .AND. n.GT.0 ) THEN
1047 IF( ijvt.EQ.1 ) THEN
1048 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
1049 $ ldvt, vt, ldvt, work, lwork,
1050 $ dif, iinfo )
1051 END IF
1052 END IF
1053 result( 28 ) = max( result( 28 ), dif )
1054
1055
1056
1057 dif = zero
1058 div = max( mnmin*ulp*s( 1 ), unfl )
1059 DO 190 i = 1, mnmin - 1
1060 IF( ssav( i ).LT.ssav( i+1 ) )
1061 $ dif = ulpinv
1062 IF( ssav( i ).LT.zero )
1063 $ dif = ulpinv
1064 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
1065 190 CONTINUE
1066 result( 29 ) = max( result( 29 ), dif )
1067 170 CONTINUE
1068 180 CONTINUE
1069
1070
1071
1072 DO 200 i = 1, 4
1073 iseed2( i ) = iseed( i )
1074 200 CONTINUE
1075 IF( mnmin.LE.1 ) THEN
1076 il = 1
1077 iu = max( 1, mnmin )
1078 ELSE
1079 il = 1 + int( ( mnmin-1 )*
dlarnd( 1, iseed2 ) )
1080 iu = 1 + int( ( mnmin-1 )*
dlarnd( 1, iseed2 ) )
1081 IF( iu.LT.il ) THEN
1082 itemp = iu
1083 iu = il
1084 il = itemp
1085 END IF
1086 END IF
1087 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
1088 CALL dgesvdx(
'V',
'V',
'I', m, n, a, lda,
1089 $ vl, vu, il, iu, nsi, s, u, ldu,
1090 $ vt, ldvt, work, lwork, iwork,
1091 $ iinfo )
1092 IF( iinfo.NE.0 ) THEN
1093 WRITE( nout, fmt = 9995 )'GESVDX', iinfo, m, n,
1094 $ jtype, lswork, ioldsd
1095 info = abs( iinfo )
1096 RETURN
1097 END IF
1098
1099 result( 30 ) = zero
1100 result( 31 ) = zero
1101 result( 32 ) = zero
1102 CALL dbdt05( m, n, asav, lda, s, nsi, u, ldu,
1103 $ vt, ldvt, work, result( 30 ) )
1104 CALL dort01(
'Columns', m, nsi, u, ldu, work, lwork,
1105 $ result( 31 ) )
1106 CALL dort01(
'Rows', nsi, n, vt, ldvt, work, lwork,
1107 $ result( 32 ) )
1108
1109
1110
1111 IF( mnmin.GT.0 .AND. nsi.GT.1 ) THEN
1112 IF( il.NE.1 ) THEN
1113 vu = ssav( il ) +
1114 $ max( half*abs( ssav( il )-ssav( il-1 ) ),
1115 $ ulp*anorm, two*rtunfl )
1116 ELSE
1117 vu = ssav( 1 ) +
1118 $ max( half*abs( ssav( ns )-ssav( 1 ) ),
1119 $ ulp*anorm, two*rtunfl )
1120 END IF
1121 IF( iu.NE.ns ) THEN
1122 vl = ssav( iu ) - max( ulp*anorm, two*rtunfl,
1123 $ half*abs( ssav( iu+1 )-ssav( iu ) ) )
1124 ELSE
1125 vl = ssav( ns ) - max( ulp*anorm, two*rtunfl,
1126 $ half*abs( ssav( ns )-ssav( 1 ) ) )
1127 END IF
1128 vl = max( vl,zero )
1129 vu = max( vu,zero )
1130 IF( vl.GE.vu ) vu = max( vu*2, vu+vl+half )
1131 ELSE
1132 vl = zero
1133 vu = one
1134 END IF
1135 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
1136 CALL dgesvdx(
'V',
'V',
'V', m, n, a, lda,
1137 $ vl, vu, il, iu, nsv, s, u, ldu,
1138 $ vt, ldvt, work, lwork, iwork,
1139 $ iinfo )
1140 IF( iinfo.NE.0 ) THEN
1141 WRITE( nout, fmt = 9995 )'GESVDX', iinfo, m, n,
1142 $ jtype, lswork, ioldsd
1143 info = abs( iinfo )
1144 RETURN
1145 END IF
1146
1147 result( 33 ) = zero
1148 result( 34 ) = zero
1149 result( 35 ) = zero
1150 CALL dbdt05( m, n, asav, lda, s, nsv, u, ldu,
1151 $ vt, ldvt, work, result( 33 ) )
1152 CALL dort01(
'Columns', m, nsv, u, ldu, work, lwork,
1153 $ result( 34 ) )
1154 CALL dort01(
'Rows', nsv, n, vt, ldvt, work, lwork,
1155 $ result( 35 ) )
1156
1157
1158
1159 DO 210 j = 1, 39
1160 IF( result( j ).GE.thresh ) THEN
1161 IF( nfail.EQ.0 ) THEN
1162 WRITE( nout, fmt = 9999 )
1163 WRITE( nout, fmt = 9998 )
1164 END IF
1165 WRITE( nout, fmt = 9997 )m, n, jtype, iws, ioldsd,
1166 $ j, result( j )
1167 nfail = nfail + 1
1168 END IF
1169 210 CONTINUE
1170 ntest = ntest + 39
1171 220 CONTINUE
1172 230 CONTINUE
1173 240 CONTINUE
1174
1175
1176
1177 CALL alasvm( path, nout, nfail, ntest, 0 )
1178
1179 9999 FORMAT( ' SVD -- Real Singular Value Decomposition Driver ',
1180 $ / ' Matrix types (see DDRVBD for details):',
1181 $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
1182 $ / ' 3 = Evenly spaced singular values near 1',
1183 $ / ' 4 = Evenly spaced singular values near underflow',
1184 $ / ' 5 = Evenly spaced singular values near overflow', / /
1185 $ ' Tests performed: ( A is dense, U and V are orthogonal,',
1186 $ / 19x, ' S is an array, and Upartial, VTpartial, and',
1187 $ / 19x, ' Spartial are partially computed U, VT and S),', / )
1188 9998 FORMAT( ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1189 $ / ' 2 = | I - U**T U | / ( M ulp ) ',
1190 $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
1191 $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
1192 $ ' decreasing order, else 1/ulp',
1193 $ / ' 5 = | U - Upartial | / ( M ulp )',
1194 $ / ' 6 = | VT - VTpartial | / ( N ulp )',
1195 $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
1196 $ / ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1197 $ / ' 9 = | I - U**T U | / ( M ulp ) ',
1198 $ / '10 = | I - VT VT**T | / ( N ulp ) ',
1199 $ / '11 = 0 if S contains min(M,N) nonnegative values in',
1200 $ ' decreasing order, else 1/ulp',
1201 $ / '12 = | U - Upartial | / ( M ulp )',
1202 $ / '13 = | VT - VTpartial | / ( N ulp )',
1203 $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )',
1204 $ / '15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1205 $ / '16 = | I - U**T U | / ( M ulp ) ',
1206 $ / '17 = | I - VT VT**T | / ( N ulp ) ',
1207 $ / '18 = 0 if S contains min(M,N) nonnegative values in',
1208 $ ' decreasing order, else 1/ulp',
1209 $ / '19 = | U - Upartial | / ( M ulp )',
1210 $ / '20 = | VT - VTpartial | / ( N ulp )',
1211 $ / '21 = | S - Spartial | / ( min(M,N) ulp |S| )',
1212 $ / '22 = 0 if S contains min(M,N) nonnegative values in',
1213 $ ' decreasing order, else 1/ulp',
1214 $ / '23 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ),',
1215 $ ' DGESVDX(V,V,A) ',
1216 $ / '24 = | I - U**T U | / ( M ulp ) ',
1217 $ / '25 = | I - VT VT**T | / ( N ulp ) ',
1218 $ / '26 = 0 if S contains min(M,N) nonnegative values in',
1219 $ ' decreasing order, else 1/ulp',
1220 $ / '27 = | U - Upartial | / ( M ulp )',
1221 $ / '28 = | VT - VTpartial | / ( N ulp )',
1222 $ / '29 = | S - Spartial | / ( min(M,N) ulp |S| )',
1223 $ / '30 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
1224 $ ' DGESVDX(V,V,I) ',
1225 $ / '31 = | I - U**T U | / ( M ulp ) ',
1226 $ / '32 = | I - VT VT**T | / ( N ulp ) ',
1227 $ / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
1228 $ ' DGESVDX(V,V,V) ',
1229 $ / '34 = | I - U**T U | / ( M ulp ) ',
1230 $ / '35 = | I - VT VT**T | / ( N ulp ) ',
1231 $ ' DGESVDQ(H,N,N,A,A',
1232 $ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1233 $ / '37 = | I - U**T U | / ( M ulp ) ',
1234 $ / '38 = | I - VT VT**T | / ( N ulp ) ',
1235 $ / '39 = 0 if S contains min(M,N) nonnegative values in',
1236 $ ' decreasing order, else 1/ulp',
1237 $ / / )
1238 9997 FORMAT( ' M=', i5, ', N=', i5, ', type ', i1, ', IWS=', i1,
1239 $ ', seed=', 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1240 9996 FORMAT( ' DDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1241 $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1242 $ i5, ')' )
1243 9995 FORMAT( ' DDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1244 $ i6, ', N=', i6, ', JTYPE=', i6, ', LSWORK=', i6, / 9x,
1245 $ 'ISEED=(', 3( i5, ',' ), i5, ')' )
1246
1247 RETURN
1248
1249
1250
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine dbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
DBDT01
subroutine dort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
DORT01
subroutine dort03(RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK, RESULT, INFO)
DORT03
subroutine dbdt05(M, N, A, LDA, S, NS, U, LDU, VT, LDVT, WORK, RESID)
DBDT05
double precision function dlarnd(IDIST, ISEED)
DLARND
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
subroutine dgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
DGESVJ
subroutine dgesvdq(JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA, S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK, WORK, LWORK, RWORK, LRWORK, INFO)
DGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE...
subroutine dgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
DGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine dgejsv(JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, WORK, LWORK, IWORK, INFO)
DGEJSV
subroutine dgesvdx(JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, IL, IU, NS, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
DGESVDX computes the singular value decomposition (SVD) for GE matrices
subroutine dgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
DGESDD