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
500 ulp =
dlamch(
'Precision' )
501 rtunfl = sqrt( unfl )
502 ulpinv = one / ulp
503 infot = 0
504
505
506
507 DO 240 jsize = 1, nsizes
508 m = mm( jsize )
509 n = nn( jsize )
510 mnmin = min( m, n )
511
512 IF( nsizes.NE.1 ) THEN
513 mtypes = min( maxtyp, ntypes )
514 ELSE
515 mtypes = min( maxtyp+1, ntypes )
516 END IF
517
518 DO 230 jtype = 1, mtypes
519 IF( .NOT.dotype( jtype ) )
520 $ GO TO 230
521
522 DO 20 j = 1, 4
523 ioldsd( j ) = iseed( j )
524 20 CONTINUE
525
526
527
528 IF( mtypes.GT.maxtyp )
529 $ GO TO 30
530
531 IF( jtype.EQ.1 ) THEN
532
533
534
535 CALL dlaset(
'Full', m, n, zero, zero, a, lda )
536
537 ELSE IF( jtype.EQ.2 ) THEN
538
539
540
541 CALL dlaset(
'Full', m, n, zero, one, a, lda )
542
543 ELSE
544
545
546
547 IF( jtype.EQ.3 )
548 $ anorm = one
549 IF( jtype.EQ.4 )
550 $ anorm = unfl / ulp
551 IF( jtype.EQ.5 )
552 $ anorm = ovfl*ulp
553 CALL dlatms( m, n,
'U', iseed,
'N', s, 4, dble( mnmin ),
554 $ anorm, m-1, n-1, 'N', a, lda, work, iinfo )
555 IF( iinfo.NE.0 ) THEN
556 WRITE( nout, fmt = 9996 )'Generator', iinfo, m, n,
557 $ jtype, ioldsd
558 info = abs( iinfo )
559 RETURN
560 END IF
561 END IF
562
563 30 CONTINUE
564 CALL dlacpy(
'F', m, n, a, lda, asav, lda )
565
566
567
568 DO 220 iws = 1, 4
569
570 DO 40 j = 1, 32
571 result( j ) = -one
572 40 CONTINUE
573
574
575
576 iwtmp = max( 3*min( m, n )+max( m, n ), 5*min( m, n ) )
577 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
578 lswork = min( lswork, lwork )
579 lswork = max( lswork, 1 )
580 IF( iws.EQ.4 )
581 $ lswork = lwork
582
583 IF( iws.GT.1 )
584 $
CALL dlacpy(
'F', m, n, asav, lda, a, lda )
585 srnamt = 'DGESVD'
586 CALL dgesvd(
'A',
'A', m, n, a, lda, ssav, usav, ldu,
587 $ vtsav, ldvt, work, lswork, iinfo )
588 IF( iinfo.NE.0 ) THEN
589 WRITE( nout, fmt = 9995 )'GESVD', iinfo, m, n, jtype,
590 $ lswork, ioldsd
591 info = abs( iinfo )
592 RETURN
593 END IF
594
595
596
597 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
598 $ vtsav, ldvt, work, result( 1 ) )
599 IF( m.NE.0 .AND. n.NE.0 ) THEN
600 CALL dort01(
'Columns', m, m, usav, ldu, work, lwork,
601 $ result( 2 ) )
602 CALL dort01(
'Rows', n, n, vtsav, ldvt, work, lwork,
603 $ result( 3 ) )
604 END IF
605 result( 4 ) = zero
606 DO 50 i = 1, mnmin - 1
607 IF( ssav( i ).LT.ssav( i+1 ) )
608 $ result( 4 ) = ulpinv
609 IF( ssav( i ).LT.zero )
610 $ result( 4 ) = ulpinv
611 50 CONTINUE
612 IF( mnmin.GE.1 ) THEN
613 IF( ssav( mnmin ).LT.zero )
614 $ result( 4 ) = ulpinv
615 END IF
616
617
618
619 result( 5 ) = zero
620 result( 6 ) = zero
621 result( 7 ) = zero
622 DO 80 iju = 0, 3
623 DO 70 ijvt = 0, 3
624 IF( ( iju.EQ.3 .AND. ijvt.EQ.3 ) .OR.
625 $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )GO TO 70
626 jobu = cjob( iju+1 )
627 jobvt = cjob( ijvt+1 )
628 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
629 srnamt = 'DGESVD'
630 CALL dgesvd( jobu, jobvt, m, n, a, lda, s, u, ldu,
631 $ vt, ldvt, work, lswork, iinfo )
632
633
634
635 dif = zero
636 IF( m.GT.0 .AND. n.GT.0 ) THEN
637 IF( iju.EQ.1 ) THEN
638 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
639 $ ldu, a, lda, work, lwork, dif,
640 $ iinfo )
641 ELSE IF( iju.EQ.2 ) THEN
642 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
643 $ ldu, u, ldu, work, lwork, dif,
644 $ iinfo )
645 ELSE IF( iju.EQ.3 ) THEN
646 CALL dort03(
'C', m, m, m, mnmin, usav, ldu,
647 $ u, ldu, work, lwork, dif,
648 $ iinfo )
649 END IF
650 END IF
651 result( 5 ) = max( result( 5 ), dif )
652
653
654
655 dif = zero
656 IF( m.GT.0 .AND. n.GT.0 ) THEN
657 IF( ijvt.EQ.1 ) THEN
658 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
659 $ ldvt, a, lda, work, lwork, dif,
660 $ iinfo )
661 ELSE IF( ijvt.EQ.2 ) THEN
662 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
663 $ ldvt, vt, ldvt, work, lwork,
664 $ dif, iinfo )
665 ELSE IF( ijvt.EQ.3 ) THEN
666 CALL dort03(
'R', n, n, n, mnmin, vtsav,
667 $ ldvt, vt, ldvt, work, lwork,
668 $ dif, iinfo )
669 END IF
670 END IF
671 result( 6 ) = max( result( 6 ), dif )
672
673
674
675 dif = zero
676 div = max( mnmin*ulp*s( 1 ), unfl )
677 DO 60 i = 1, mnmin - 1
678 IF( ssav( i ).LT.ssav( i+1 ) )
679 $ dif = ulpinv
680 IF( ssav( i ).LT.zero )
681 $ dif = ulpinv
682 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
683 60 CONTINUE
684 result( 7 ) = max( result( 7 ), dif )
685 70 CONTINUE
686 80 CONTINUE
687
688
689
690 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
691 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
692 lswork = min( lswork, lwork )
693 lswork = max( lswork, 1 )
694 IF( iws.EQ.4 )
695 $ lswork = lwork
696
697 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
698 srnamt = 'DGESDD'
699 CALL dgesdd(
'A', m, n, a, lda, ssav, usav, ldu, vtsav,
700 $ ldvt, work, lswork, iwork, iinfo )
701 IF( iinfo.NE.0 ) THEN
702 WRITE( nout, fmt = 9995 )'GESDD', iinfo, m, n, jtype,
703 $ lswork, ioldsd
704 info = abs( iinfo )
705 RETURN
706 END IF
707
708
709
710 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
711 $ vtsav, ldvt, work, result( 8 ) )
712 IF( m.NE.0 .AND. n.NE.0 ) THEN
713 CALL dort01(
'Columns', m, m, usav, ldu, work, lwork,
714 $ result( 9 ) )
715 CALL dort01(
'Rows', n, n, vtsav, ldvt, work, lwork,
716 $ result( 10 ) )
717 END IF
718 result( 11 ) = zero
719 DO 90 i = 1, mnmin - 1
720 IF( ssav( i ).LT.ssav( i+1 ) )
721 $ result( 11 ) = ulpinv
722 IF( ssav( i ).LT.zero )
723 $ result( 11 ) = ulpinv
724 90 CONTINUE
725 IF( mnmin.GE.1 ) THEN
726 IF( ssav( mnmin ).LT.zero )
727 $ result( 11 ) = ulpinv
728 END IF
729
730
731
732 result( 12 ) = zero
733 result( 13 ) = zero
734 result( 14 ) = zero
735 DO 110 ijq = 0, 2
736 jobq = cjob( ijq+1 )
737 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
738 srnamt = 'DGESDD'
739 CALL dgesdd( jobq, m, n, a, lda, s, u, ldu, vt, ldvt,
740 $ work, lswork, iwork, iinfo )
741
742
743
744 dif = zero
745 IF( m.GT.0 .AND. n.GT.0 ) THEN
746 IF( ijq.EQ.1 ) THEN
747 IF( m.GE.n ) THEN
748 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
749 $ ldu, a, lda, work, lwork, dif,
750 $ info )
751 ELSE
752 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
753 $ ldu, u, ldu, work, lwork, dif,
754 $ info )
755 END IF
756 ELSE IF( ijq.EQ.2 ) THEN
757 CALL dort03(
'C', m, mnmin, m, mnmin, usav, ldu,
758 $ u, ldu, work, lwork, dif, info )
759 END IF
760 END IF
761 result( 12 ) = max( result( 12 ), dif )
762
763
764
765 dif = zero
766 IF( m.GT.0 .AND. n.GT.0 ) THEN
767 IF( ijq.EQ.1 ) THEN
768 IF( m.GE.n ) THEN
769 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
770 $ ldvt, vt, ldvt, work, lwork,
771 $ dif, info )
772 ELSE
773 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
774 $ ldvt, a, lda, work, lwork, dif,
775 $ info )
776 END IF
777 ELSE IF( ijq.EQ.2 ) THEN
778 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
779 $ ldvt, vt, ldvt, work, lwork, dif,
780 $ info )
781 END IF
782 END IF
783 result( 13 ) = max( result( 13 ), dif )
784
785
786
787 dif = zero
788 div = max( mnmin*ulp*s( 1 ), unfl )
789 DO 100 i = 1, mnmin - 1
790 IF( ssav( i ).LT.ssav( i+1 ) )
791 $ dif = ulpinv
792 IF( ssav( i ).LT.zero )
793 $ dif = ulpinv
794 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
795 100 CONTINUE
796 result( 14 ) = max( result( 14 ), dif )
797 110 CONTINUE
798
799
800
801
802 result( 36 ) = zero
803 result( 37 ) = zero
804 result( 38 ) = zero
805 result( 39 ) = zero
806
807 IF( m.GE.n ) THEN
808 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
809 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
810 lswork = min( lswork, lwork )
811 lswork = max( lswork, 1 )
812 IF( iws.EQ.4 )
813 $ lswork = lwork
814
815 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
816 srnamt = 'DGESVDQ'
817
818 lrwork = 2
819 liwork = max( n, 1 )
820 CALL dgesvdq(
'H',
'N',
'N',
'A',
'A',
821 $ m, n, a, lda, ssav, usav, ldu,
822 $ vtsav, ldvt, numrank, iwork, liwork,
823 $ work, lwork, rwork, lrwork, iinfo )
824
825 IF( iinfo.NE.0 ) THEN
826 WRITE( nout, fmt = 9995 )'DGESVDQ', iinfo, m, n,
827 $ jtype, lswork, ioldsd
828 info = abs( iinfo )
829 RETURN
830 END IF
831
832
833
834 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
835 $ vtsav, ldvt, work, result( 36 ) )
836 IF( m.NE.0 .AND. n.NE.0 ) THEN
837 CALL dort01(
'Columns', m, m, usav, ldu, work,
838 $ lwork, result( 37 ) )
839 CALL dort01(
'Rows', n, n, vtsav, ldvt, work,
840 $ lwork, result( 38 ) )
841 END IF
842 result( 39 ) = zero
843 DO 199 i = 1, mnmin - 1
844 IF( ssav( i ).LT.ssav( i+1 ) )
845 $ result( 39 ) = ulpinv
846 IF( ssav( i ).LT.zero )
847 $ result( 39 ) = ulpinv
848 199 CONTINUE
849 IF( mnmin.GE.1 ) THEN
850 IF( ssav( mnmin ).LT.zero )
851 $ result( 39 ) = ulpinv
852 END IF
853 END IF
854
855
856
857
858 result( 15 ) = zero
859 result( 16 ) = zero
860 result( 17 ) = zero
861 result( 18 ) = zero
862
863 IF( m.GE.n ) THEN
864 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
865 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
866 lswork = min( lswork, lwork )
867 lswork = max( lswork, 1 )
868 IF( iws.EQ.4 )
869 $ lswork = lwork
870
871 CALL dlacpy(
'F', m, n, asav, lda, usav, lda )
872 srnamt = 'DGESVJ'
873 CALL dgesvj(
'G',
'U',
'V', m, n, usav, lda, ssav,
874 & 0, a, ldvt, work, lwork, info )
875
876
877
878 DO j=1,n
879 DO i=1,n
880 vtsav(j,i) = a(i,j)
881 END DO
882 END DO
883
884 IF( iinfo.NE.0 ) THEN
885 WRITE( nout, fmt = 9995 )'GESVJ', iinfo, m, n,
886 $ jtype, lswork, ioldsd
887 info = abs( iinfo )
888 RETURN
889 END IF
890
891
892
893 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
894 $ vtsav, ldvt, work, result( 15 ) )
895 IF( m.NE.0 .AND. n.NE.0 ) THEN
896 CALL dort01(
'Columns', m, m, usav, ldu, work,
897 $ lwork, result( 16 ) )
898 CALL dort01(
'Rows', n, n, vtsav, ldvt, work,
899 $ lwork, result( 17 ) )
900 END IF
901 result( 18 ) = zero
902 DO 120 i = 1, mnmin - 1
903 IF( ssav( i ).LT.ssav( i+1 ) )
904 $ result( 18 ) = ulpinv
905 IF( ssav( i ).LT.zero )
906 $ result( 18 ) = ulpinv
907 120 CONTINUE
908 IF( mnmin.GE.1 ) THEN
909 IF( ssav( mnmin ).LT.zero )
910 $ result( 18 ) = ulpinv
911 END IF
912 END IF
913
914
915
916
917 result( 19 ) = zero
918 result( 20 ) = zero
919 result( 21 ) = zero
920 result( 22 ) = zero
921 IF( m.GE.n ) THEN
922 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
923 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
924 lswork = min( lswork, lwork )
925 lswork = max( lswork, 1 )
926 IF( iws.EQ.4 )
927 $ lswork = lwork
928
929 CALL dlacpy(
'F', m, n, asav, lda, vtsav, lda )
930 srnamt = 'DGEJSV'
931 CALL dgejsv(
'G',
'U',
'V',
'R',
'N',
'N',
932 & m, n, vtsav, lda, ssav, usav, ldu, a, ldvt,
933 & work, lwork, iwork, info )
934
935
936
937 DO 140 j=1,n
938 DO 130 i=1,n
939 vtsav(j,i) = a(i,j)
940 130 END DO
941 140 END DO
942
943 IF( iinfo.NE.0 ) THEN
944 WRITE( nout, fmt = 9995 )'GEJSV', iinfo, m, n,
945 $ jtype, lswork, ioldsd
946 info = abs( iinfo )
947 RETURN
948 END IF
949
950
951
952 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
953 $ vtsav, ldvt, work, result( 19 ) )
954 IF( m.NE.0 .AND. n.NE.0 ) THEN
955 CALL dort01(
'Columns', m, m, usav, ldu, work,
956 $ lwork, result( 20 ) )
957 CALL dort01(
'Rows', n, n, vtsav, ldvt, work,
958 $ lwork, result( 21 ) )
959 END IF
960 result( 22 ) = zero
961 DO 150 i = 1, mnmin - 1
962 IF( ssav( i ).LT.ssav( i+1 ) )
963 $ result( 22 ) = ulpinv
964 IF( ssav( i ).LT.zero )
965 $ result( 22 ) = ulpinv
966 150 CONTINUE
967 IF( mnmin.GE.1 ) THEN
968 IF( ssav( mnmin ).LT.zero )
969 $ result( 22 ) = ulpinv
970 END IF
971 END IF
972
973
974
975 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
976 CALL dgesvdx(
'V',
'V',
'A', m, n, a, lda,
977 $ vl, vu, il, iu, ns, ssav, usav, ldu,
978 $ vtsav, ldvt, work, lwork, iwork,
979 $ iinfo )
980 IF( iinfo.NE.0 ) THEN
981 WRITE( nout, fmt = 9995 )'GESVDX', iinfo, m, n,
982 $ jtype, lswork, ioldsd
983 info = abs( iinfo )
984 RETURN
985 END IF
986
987
988
989 result( 23 ) = zero
990 result( 24 ) = zero
991 result( 25 ) = zero
992 CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
993 $ vtsav, ldvt, work, result( 23 ) )
994 IF( m.NE.0 .AND. n.NE.0 ) THEN
995 CALL dort01(
'Columns', m, m, usav, ldu, work, lwork,
996 $ result( 24 ) )
997 CALL dort01(
'Rows', n, n, vtsav, ldvt, work, lwork,
998 $ result( 25 ) )
999 END IF
1000 result( 26 ) = zero
1001 DO 160 i = 1, mnmin - 1
1002 IF( ssav( i ).LT.ssav( i+1 ) )
1003 $ result( 26 ) = ulpinv
1004 IF( ssav( i ).LT.zero )
1005 $ result( 26 ) = ulpinv
1006 160 CONTINUE
1007 IF( mnmin.GE.1 ) THEN
1008 IF( ssav( mnmin ).LT.zero )
1009 $ result( 26 ) = ulpinv
1010 END IF
1011
1012
1013
1014 result( 27 ) = zero
1015 result( 28 ) = zero
1016 result( 29 ) = zero
1017 DO 180 iju = 0, 1
1018 DO 170 ijvt = 0, 1
1019 IF( ( iju.EQ.0 .AND. ijvt.EQ.0 ) .OR.
1020 $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )GO TO 170
1021 jobu = cjobv( iju+1 )
1022 jobvt = cjobv( ijvt+1 )
1023 range = cjobr( 1 )
1024 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
1025 CALL dgesvdx( jobu, jobvt, range, m, n, a, lda,
1026 $ vl, vu, il, iu, ns, s, u, ldu,
1027 $ vt, ldvt, work, lwork, iwork,
1028 $ iinfo )
1029
1030
1031
1032 dif = zero
1033 IF( m.GT.0 .AND. n.GT.0 ) THEN
1034 IF( iju.EQ.1 ) THEN
1035 CALL dort03(
'C', m, mnmin, m, mnmin, usav,
1036 $ ldu, u, ldu, work, lwork, dif,
1037 $ iinfo )
1038 END IF
1039 END IF
1040 result( 27 ) = max( result( 27 ), dif )
1041
1042
1043
1044 dif = zero
1045 IF( m.GT.0 .AND. n.GT.0 ) THEN
1046 IF( ijvt.EQ.1 ) THEN
1047 CALL dort03(
'R', n, mnmin, n, mnmin, vtsav,
1048 $ ldvt, vt, ldvt, work, lwork,
1049 $ dif, iinfo )
1050 END IF
1051 END IF
1052 result( 28 ) = max( result( 28 ), dif )
1053
1054
1055
1056 dif = zero
1057 div = max( mnmin*ulp*s( 1 ), unfl )
1058 DO 190 i = 1, mnmin - 1
1059 IF( ssav( i ).LT.ssav( i+1 ) )
1060 $ dif = ulpinv
1061 IF( ssav( i ).LT.zero )
1062 $ dif = ulpinv
1063 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
1064 190 CONTINUE
1065 result( 29 ) = max( result( 29 ), dif )
1066 170 CONTINUE
1067 180 CONTINUE
1068
1069
1070
1071 DO 200 i = 1, 4
1072 iseed2( i ) = iseed( i )
1073 200 CONTINUE
1074 IF( mnmin.LE.1 ) THEN
1075 il = 1
1076 iu = max( 1, mnmin )
1077 ELSE
1078 il = 1 + int( ( mnmin-1 )*
dlarnd( 1, iseed2 ) )
1079 iu = 1 + int( ( mnmin-1 )*
dlarnd( 1, iseed2 ) )
1080 IF( iu.LT.il ) THEN
1081 itemp = iu
1082 iu = il
1083 il = itemp
1084 END IF
1085 END IF
1086 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
1087 CALL dgesvdx(
'V',
'V',
'I', m, n, a, lda,
1088 $ vl, vu, il, iu, nsi, s, u, ldu,
1089 $ vt, ldvt, work, lwork, iwork,
1090 $ iinfo )
1091 IF( iinfo.NE.0 ) THEN
1092 WRITE( nout, fmt = 9995 )'GESVDX', iinfo, m, n,
1093 $ jtype, lswork, ioldsd
1094 info = abs( iinfo )
1095 RETURN
1096 END IF
1097
1098 result( 30 ) = zero
1099 result( 31 ) = zero
1100 result( 32 ) = zero
1101 CALL dbdt05( m, n, asav, lda, s, nsi, u, ldu,
1102 $ vt, ldvt, work, result( 30 ) )
1103 CALL dort01(
'Columns', m, nsi, u, ldu, work, lwork,
1104 $ result( 31 ) )
1105 CALL dort01(
'Rows', nsi, n, vt, ldvt, work, lwork,
1106 $ result( 32 ) )
1107
1108
1109
1110 IF( mnmin.GT.0 .AND. nsi.GT.1 ) THEN
1111 IF( il.NE.1 ) THEN
1112 vu = ssav( il ) +
1113 $ max( half*abs( ssav( il )-ssav( il-1 ) ),
1114 $ ulp*anorm, two*rtunfl )
1115 ELSE
1116 vu = ssav( 1 ) +
1117 $ max( half*abs( ssav( ns )-ssav( 1 ) ),
1118 $ ulp*anorm, two*rtunfl )
1119 END IF
1120 IF( iu.NE.ns ) THEN
1121 vl = ssav( iu ) - max( ulp*anorm, two*rtunfl,
1122 $ half*abs( ssav( iu+1 )-ssav( iu ) ) )
1123 ELSE
1124 vl = ssav( ns ) - max( ulp*anorm, two*rtunfl,
1125 $ half*abs( ssav( ns )-ssav( 1 ) ) )
1126 END IF
1127 vl = max( vl,zero )
1128 vu = max( vu,zero )
1129 IF( vl.GE.vu ) vu = max( vu*2, vu+vl+half )
1130 ELSE
1131 vl = zero
1132 vu = one
1133 END IF
1134 CALL dlacpy(
'F', m, n, asav, lda, a, lda )
1135 CALL dgesvdx(
'V',
'V',
'V', m, n, a, lda,
1136 $ vl, vu, il, iu, nsv, s, u, ldu,
1137 $ vt, ldvt, work, lwork, iwork,
1138 $ iinfo )
1139 IF( iinfo.NE.0 ) THEN
1140 WRITE( nout, fmt = 9995 )'GESVDX', iinfo, m, n,
1141 $ jtype, lswork, ioldsd
1142 info = abs( iinfo )
1143 RETURN
1144 END IF
1145
1146 result( 33 ) = zero
1147 result( 34 ) = zero
1148 result( 35 ) = zero
1149 CALL dbdt05( m, n, asav, lda, s, nsv, u, ldu,
1150 $ vt, ldvt, work, result( 33 ) )
1151 CALL dort01(
'Columns', m, nsv, u, ldu, work, lwork,
1152 $ result( 34 ) )
1153 CALL dort01(
'Rows', nsv, n, vt, ldvt, work, lwork,
1154 $ result( 35 ) )
1155
1156
1157
1158 DO 210 j = 1, 39
1159 IF( result( j ).GE.thresh ) THEN
1160 IF( nfail.EQ.0 ) THEN
1161 WRITE( nout, fmt = 9999 )
1162 WRITE( nout, fmt = 9998 )
1163 END IF
1164 WRITE( nout, fmt = 9997 )m, n, jtype, iws, ioldsd,
1165 $ j, result( j )
1166 nfail = nfail + 1
1167 END IF
1168 210 CONTINUE
1169 ntest = ntest + 39
1170 220 CONTINUE
1171 230 CONTINUE
1172 240 CONTINUE
1173
1174
1175
1176 CALL alasvm( path, nout, nfail, ntest, 0 )
1177
1178 9999 FORMAT( ' SVD -- Real Singular Value Decomposition Driver ',
1179 $ / ' Matrix types (see DDRVBD for details):',
1180 $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
1181 $ / ' 3 = Evenly spaced singular values near 1',
1182 $ / ' 4 = Evenly spaced singular values near underflow',
1183 $ / ' 5 = Evenly spaced singular values near overflow', / /
1184 $ ' Tests performed: ( A is dense, U and V are orthogonal,',
1185 $ / 19x, ' S is an array, and Upartial, VTpartial, and',
1186 $ / 19x, ' Spartial are partially computed U, VT and S),', / )
1187 9998 FORMAT( ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1188 $ / ' 2 = | I - U**T U | / ( M ulp ) ',
1189 $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
1190 $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
1191 $ ' decreasing order, else 1/ulp',
1192 $ / ' 5 = | U - Upartial | / ( M ulp )',
1193 $ / ' 6 = | VT - VTpartial | / ( N ulp )',
1194 $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
1195 $ / ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1196 $ / ' 9 = | I - U**T U | / ( M ulp ) ',
1197 $ / '10 = | I - VT VT**T | / ( N ulp ) ',
1198 $ / '11 = 0 if S contains min(M,N) nonnegative values in',
1199 $ ' decreasing order, else 1/ulp',
1200 $ / '12 = | U - Upartial | / ( M ulp )',
1201 $ / '13 = | VT - VTpartial | / ( N ulp )',
1202 $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )',
1203 $ / '15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1204 $ / '16 = | I - U**T U | / ( M ulp ) ',
1205 $ / '17 = | I - VT VT**T | / ( N ulp ) ',
1206 $ / '18 = 0 if S contains min(M,N) nonnegative values in',
1207 $ ' decreasing order, else 1/ulp',
1208 $ / '19 = | U - Upartial | / ( M ulp )',
1209 $ / '20 = | VT - VTpartial | / ( N ulp )',
1210 $ / '21 = | S - Spartial | / ( min(M,N) ulp |S| )',
1211 $ / '22 = 0 if S contains min(M,N) nonnegative values in',
1212 $ ' decreasing order, else 1/ulp',
1213 $ / '23 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ),',
1214 $ ' DGESVDX(V,V,A) ',
1215 $ / '24 = | I - U**T U | / ( M ulp ) ',
1216 $ / '25 = | I - VT VT**T | / ( N ulp ) ',
1217 $ / '26 = 0 if S contains min(M,N) nonnegative values in',
1218 $ ' decreasing order, else 1/ulp',
1219 $ / '27 = | U - Upartial | / ( M ulp )',
1220 $ / '28 = | VT - VTpartial | / ( N ulp )',
1221 $ / '29 = | S - Spartial | / ( min(M,N) ulp |S| )',
1222 $ / '30 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
1223 $ ' DGESVDX(V,V,I) ',
1224 $ / '31 = | I - U**T U | / ( M ulp ) ',
1225 $ / '32 = | I - VT VT**T | / ( N ulp ) ',
1226 $ / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
1227 $ ' DGESVDX(V,V,V) ',
1228 $ / '34 = | I - U**T U | / ( M ulp ) ',
1229 $ / '35 = | I - VT VT**T | / ( N ulp ) ',
1230 $ ' DGESVDQ(H,N,N,A,A',
1231 $ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1232 $ / '37 = | I - U**T U | / ( M ulp ) ',
1233 $ / '38 = | I - VT VT**T | / ( N ulp ) ',
1234 $ / '39 = 0 if S contains min(M,N) nonnegative values in',
1235 $ ' decreasing order, else 1/ulp',
1236 $ / / )
1237 9997 FORMAT( ' M=', i5, ', N=', i5, ', type ', i1, ', IWS=', i1,
1238 $ ', seed=', 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1239 9996 FORMAT( ' DDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1240 $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1241 $ i5, ')' )
1242 9995 FORMAT( ' DDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1243 $ i6, ', N=', i6, ', JTYPE=', i6, ', LSWORK=', i6, / 9x,
1244 $ 'ISEED=(', 3( i5, ',' ), i5, ')' )
1245
1246 RETURN
1247
1248
1249
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
subroutine xerbla(srname, info)
subroutine dbdt01(m, n, kd, a, lda, q, ldq, d, e, pt, ldpt, work, resid)
DBDT01
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 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 dgejsv(joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, work, lwork, iwork, info)
DGEJSV
subroutine dgesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
DGESDD
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 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 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 dgesvj(joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, work, lwork, info)
DGESVJ
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
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.