318 SUBROUTINE ddrvbd( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
319 $ a, lda, u, ldu, vt, ldvt, asav, usav, vtsav, s,
320 $ ssav, e, work, lwork, iwork, nout, info )
328 INTEGER info, lda, ldu, ldvt, lwork, nout, nsizes,
330 DOUBLE PRECISION thresh
334 INTEGER iseed( 4 ), iwork( * ), mm( * ), nn( * )
335 DOUBLE PRECISION a( lda, * ), asav( lda, * ), e( * ), s( * ),
336 $ ssav( * ), u( ldu, * ), usav( ldu, * ),
337 $ vt( ldvt, * ), vtsav( ldvt, * ), work( * )
343 DOUBLE PRECISION zero, one
344 parameter( zero = 0.0d0, one = 1.0d0 )
346 parameter( maxtyp = 5 )
350 CHARACTER jobq, jobu, jobvt
352 INTEGER i, iinfo, ijq, iju, ijvt, iws, iwtmp, j, jsize,
353 $ jtype, lswork, m, minwrk, mmax, mnmax, mnmin,
354 $ mtypes, n, nfail, nmax, ntest
355 DOUBLE PRECISION anorm, dif, div, ovfl, ulp, ulpinv, unfl
360 DOUBLE PRECISION result( 22 )
372 INTRINSIC abs, dble, max, min
380 common / infoc / infot, nunit, ok, lerr
381 common / srnamc / srnamt
384 DATA cjob /
'N',
'O',
'S',
'A' /
398 mmax = max( mmax, mm( j ) )
401 nmax = max( nmax, nn( j ) )
404 mnmax = max( mnmax, min( mm( j ), nn( j ) ) )
405 minwrk = max( minwrk, max( 3*min( mm( j ),
406 $ nn( j ) )+max( mm( j ), nn( j ) ), 5*min( mm( j ),
407 $ nn( j )-4 ) )+2*min( mm( j ), nn( j ) )**2 )
412 IF( nsizes.LT.0 )
THEN
414 ELSE IF( badmm )
THEN
416 ELSE IF( badnn )
THEN
418 ELSE IF( ntypes.LT.0 )
THEN
420 ELSE IF( lda.LT.max( 1, mmax ) )
THEN
422 ELSE IF( ldu.LT.max( 1, mmax ) )
THEN
424 ELSE IF( ldvt.LT.max( 1, nmax ) )
THEN
426 ELSE IF( minwrk.GT.lwork )
THEN
431 CALL
xerbla(
'DDRVBD', -info )
437 path( 1: 1 ) =
'Double precision'
441 unfl =
dlamch(
'Safe minimum' )
444 ulp =
dlamch(
'Precision' )
450 DO 150 jsize = 1, nsizes
455 IF( nsizes.NE.1 )
THEN
456 mtypes = min( maxtyp, ntypes )
458 mtypes = min( maxtyp+1, ntypes )
461 DO 140 jtype = 1, mtypes
462 IF( .NOT.dotype( jtype ) )
466 ioldsd( j ) = iseed( j )
471 IF( mtypes.GT.maxtyp )
474 IF( jtype.EQ.1 )
THEN
478 CALL
dlaset(
'Full', m, n, zero, zero, a, lda )
480 ELSE IF( jtype.EQ.2 )
THEN
484 CALL
dlaset(
'Full', m, n, zero, one, a, lda )
496 CALL
dlatms( m, n,
'U', iseed,
'N', s, 4, dble( mnmin ),
497 $ anorm, m-1, n-1,
'N', a, lda, work, iinfo )
498 IF( iinfo.NE.0 )
THEN
499 WRITE( nout, fmt = 9996 )
'Generator', iinfo, m, n,
507 CALL
dlacpy(
'F', m, n, a, lda, asav, lda )
519 iwtmp = max( 3*min( m, n )+max( m, n ), 5*min( m, n ) )
520 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
521 lswork = min( lswork, lwork )
522 lswork = max( lswork, 1 )
527 $ CALL
dlacpy(
'F', m, n, asav, lda, a, lda )
529 CALL
dgesvd(
'A',
'A', m, n, a, lda, ssav, usav, ldu,
530 $ vtsav, ldvt, work, lswork, iinfo )
531 IF( iinfo.NE.0 )
THEN
532 WRITE( nout, fmt = 9995 )
'GESVD', iinfo, m, n, jtype,
540 CALL
dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
541 $ vtsav, ldvt, work, result( 1 ) )
542 IF( m.NE.0 .AND. n.NE.0 )
THEN
543 CALL
dort01(
'Columns', m, m, usav, ldu, work, lwork,
545 CALL
dort01(
'Rows', n, n, vtsav, ldvt, work, lwork,
549 DO 50 i = 1, mnmin - 1
550 IF( ssav( i ).LT.ssav( i+1 ) )
551 $ result( 4 ) = ulpinv
552 IF( ssav( i ).LT.zero )
553 $ result( 4 ) = ulpinv
555 IF( mnmin.GE.1 )
THEN
556 IF( ssav( mnmin ).LT.zero )
557 $ result( 4 ) = ulpinv
567 IF( ( iju.EQ.3 .AND. ijvt.EQ.3 ) .OR.
568 $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )go to 70
570 jobvt = cjob( ijvt+1 )
571 CALL
dlacpy(
'F', m, n, asav, lda, a, lda )
573 CALL
dgesvd( jobu, jobvt, m, n, a, lda, s, u, ldu,
574 $ vt, ldvt, work, lswork, iinfo )
579 IF( m.GT.0 .AND. n.GT.0 )
THEN
581 CALL
dort03(
'C', m, mnmin, m, mnmin, usav,
582 $ ldu, a, lda, work, lwork, dif,
584 ELSE IF( iju.EQ.2 )
THEN
585 CALL
dort03(
'C', m, mnmin, m, mnmin, usav,
586 $ ldu, u, ldu, work, lwork, dif,
588 ELSE IF( iju.EQ.3 )
THEN
589 CALL
dort03(
'C', m, m, m, mnmin, usav, ldu,
590 $ u, ldu, work, lwork, dif,
594 result( 5 ) = max( result( 5 ), dif )
599 IF( m.GT.0 .AND. n.GT.0 )
THEN
601 CALL
dort03(
'R', n, mnmin, n, mnmin, vtsav,
602 $ ldvt, a, lda, work, lwork, dif,
604 ELSE IF( ijvt.EQ.2 )
THEN
605 CALL
dort03(
'R', n, mnmin, n, mnmin, vtsav,
606 $ ldvt, vt, ldvt, work, lwork,
608 ELSE IF( ijvt.EQ.3 )
THEN
609 CALL
dort03(
'R', n, n, n, mnmin, vtsav,
610 $ ldvt, vt, ldvt, work, lwork,
614 result( 6 ) = max( result( 6 ), dif )
619 div = max( dble( mnmin )*ulp*s( 1 ), unfl )
620 DO 60 i = 1, mnmin - 1
621 IF( ssav( i ).LT.ssav( i+1 ) )
623 IF( ssav( i ).LT.zero )
625 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
627 result( 7 ) = max( result( 7 ), dif )
633 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
634 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
635 lswork = min( lswork, lwork )
636 lswork = max( lswork, 1 )
640 CALL
dlacpy(
'F', m, n, asav, lda, a, lda )
642 CALL
dgesdd(
'A', m, n, a, lda, ssav, usav, ldu, vtsav,
643 $ ldvt, work, lswork, iwork, iinfo )
644 IF( iinfo.NE.0 )
THEN
645 WRITE( nout, fmt = 9995 )
'GESDD', iinfo, m, n, jtype,
653 CALL
dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
654 $ vtsav, ldvt, work, result( 8 ) )
655 IF( m.NE.0 .AND. n.NE.0 )
THEN
656 CALL
dort01(
'Columns', m, m, usav, ldu, work, lwork,
658 CALL
dort01(
'Rows', n, n, vtsav, ldvt, work, lwork,
662 DO 90 i = 1, mnmin - 1
663 IF( ssav( i ).LT.ssav( i+1 ) )
664 $ result( 11 ) = ulpinv
665 IF( ssav( i ).LT.zero )
666 $ result( 11 ) = ulpinv
668 IF( mnmin.GE.1 )
THEN
669 IF( ssav( mnmin ).LT.zero )
670 $ result( 11 ) = ulpinv
680 CALL
dlacpy(
'F', m, n, asav, lda, a, lda )
682 CALL
dgesdd( jobq, m, n, a, lda, s, u, ldu, vt, ldvt,
683 $ work, lswork, iwork, iinfo )
688 IF( m.GT.0 .AND. n.GT.0 )
THEN
691 CALL
dort03(
'C', m, mnmin, m, mnmin, usav,
692 $ ldu, a, lda, work, lwork, dif,
695 CALL
dort03(
'C', m, mnmin, m, mnmin, usav,
696 $ ldu, u, ldu, work, lwork, dif,
699 ELSE IF( ijq.EQ.2 )
THEN
700 CALL
dort03(
'C', m, mnmin, m, mnmin, usav, ldu,
701 $ u, ldu, work, lwork, dif, info )
704 result( 12 ) = max( result( 12 ), dif )
709 IF( m.GT.0 .AND. n.GT.0 )
THEN
712 CALL
dort03(
'R', n, mnmin, n, mnmin, vtsav,
713 $ ldvt, vt, ldvt, work, lwork,
716 CALL
dort03(
'R', n, mnmin, n, mnmin, vtsav,
717 $ ldvt, a, lda, work, lwork, dif,
720 ELSE IF( ijq.EQ.2 )
THEN
721 CALL
dort03(
'R', n, mnmin, n, mnmin, vtsav,
722 $ ldvt, vt, ldvt, work, lwork, dif,
726 result( 13 ) = max( result( 13 ), dif )
731 div = max( dble( mnmin )*ulp*s( 1 ), unfl )
732 DO 100 i = 1, mnmin - 1
733 IF( ssav( i ).LT.ssav( i+1 ) )
735 IF( ssav( i ).LT.zero )
737 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
739 result( 14 ) = max( result( 14 ), dif )
751 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
752 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
753 lswork = min( lswork, lwork )
754 lswork = max( lswork, 1 )
758 CALL
dlacpy(
'F', m, n, asav, lda, usav, lda )
760 CALL
dgesvj(
'G',
'U',
'V', m, n, usav, lda, ssav,
761 & 0, a, ldvt, work, lwork, info )
772 IF( iinfo.NE.0 )
THEN
773 WRITE( nout, fmt = 9995 )
'GESVJ', iinfo, m, n,
774 $ jtype, lswork, ioldsd
781 CALL
dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
782 $ vtsav, ldvt, work, result( 15 ) )
783 IF( m.NE.0 .AND. n.NE.0 )
THEN
784 CALL
dort01(
'Columns', m, m, usav, ldu, work,
785 $ lwork, result( 16 ) )
786 CALL
dort01(
'Rows', n, n, vtsav, ldvt, work,
787 $ lwork, result( 17 ) )
790 DO 200 i = 1, mnmin - 1
791 IF( ssav( i ).LT.ssav( i+1 ) )
792 $ result( 18 ) = ulpinv
793 IF( ssav( i ).LT.zero )
794 $ result( 18 ) = ulpinv
796 IF( mnmin.GE.1 )
THEN
797 IF( ssav( mnmin ).LT.zero )
798 $ result( 18 ) = ulpinv
810 iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
811 lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
812 lswork = min( lswork, lwork )
813 lswork = max( lswork, 1 )
817 CALL
dlacpy(
'F', m, n, asav, lda, vtsav, lda )
819 CALL
dgejsv(
'G',
'U',
'V',
'R',
'N',
'N',
820 & m, n, vtsav, lda, ssav, usav, ldu, a, ldvt,
821 & work, lwork, iwork, info )
832 IF( iinfo.NE.0 )
THEN
833 WRITE( nout, fmt = 9995 )
'GESVJ', iinfo, m, n,
834 $ jtype, lswork, ioldsd
841 CALL
dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
842 $ vtsav, ldvt, work, result( 19 ) )
843 IF( m.NE.0 .AND. n.NE.0 )
THEN
844 CALL
dort01(
'Columns', m, m, usav, ldu, work,
845 $ lwork, result( 20 ) )
846 CALL
dort01(
'Rows', n, n, vtsav, ldvt, work,
847 $ lwork, result( 21 ) )
850 DO 300 i = 1, mnmin - 1
851 IF( ssav( i ).LT.ssav( i+1 ) )
852 $ result( 22 ) = ulpinv
853 IF( ssav( i ).LT.zero )
854 $ result( 22 ) = ulpinv
856 IF( mnmin.GE.1 )
THEN
857 IF( ssav( mnmin ).LT.zero )
858 $ result( 22 ) = ulpinv
865 IF( result( j ).GE.thresh )
THEN
866 IF( nfail.EQ.0 )
THEN
867 WRITE( nout, fmt = 9999 )
868 WRITE( nout, fmt = 9998 )
870 WRITE( nout, fmt = 9997 )m, n, jtype, iws, ioldsd,
883 CALL
alasvm( path, nout, nfail, ntest, 0 )
885 9999 format(
' SVD -- Real Singular Value Decomposition Driver ',
886 $ /
' Matrix types (see DDRVBD for details):',
887 $ / /
' 1 = Zero matrix', /
' 2 = Identity matrix',
888 $ /
' 3 = Evenly spaced singular values near 1',
889 $ /
' 4 = Evenly spaced singular values near underflow',
890 $ /
' 5 = Evenly spaced singular values near overflow', / /
891 $
' Tests performed: ( A is dense, U and V are orthogonal,',
892 $ / 19x,
' S is an array, and Upartial, VTpartial, and',
893 $ / 19x,
' Spartial are partially computed U, VT and S),', / )
894 9998 format(
' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
895 $ /
' 2 = | I - U**T U | / ( M ulp ) ',
896 $ /
' 3 = | I - VT VT**T | / ( N ulp ) ',
897 $ /
' 4 = 0 if S contains min(M,N) nonnegative values in',
898 $
' decreasing order, else 1/ulp',
899 $ /
' 5 = | U - Upartial | / ( M ulp )',
900 $ /
' 6 = | VT - VTpartial | / ( N ulp )',
901 $ /
' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
902 $ /
' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
903 $ /
' 9 = | I - U**T U | / ( M ulp ) ',
904 $ /
'10 = | I - VT VT**T | / ( N ulp ) ',
905 $ /
'11 = 0 if S contains min(M,N) nonnegative values in',
906 $
' decreasing order, else 1/ulp',
907 $ /
'12 = | U - Upartial | / ( M ulp )',
908 $ /
'13 = | VT - VTpartial | / ( N ulp )',
909 $ /
'14 = | S - Spartial | / ( min(M,N) ulp |S| )',
910 $ /
'15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
911 $ /
'16 = | I - U**T U | / ( M ulp ) ',
912 $ /
'17 = | I - VT VT**T | / ( N ulp ) ',
913 $ /
'18 = 0 if S contains min(M,N) nonnegative values in',
914 $
' decreasing order, else 1/ulp',
915 $ /
'19 = | U - Upartial | / ( M ulp )',
916 $ /
'20 = | VT - VTpartial | / ( N ulp )',
917 $ /
'21 = | S - Spartial | / ( min(M,N) ulp |S| )', / / )
918 9997 format(
' M=', i5,
', N=', i5,
', type ', i1,
', IWS=', i1,
919 $
', seed=', 4( i4,
',' ),
' test(', i2,
')=', g11.4 )
920 9996 format(
' DDRVBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
921 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
923 9995 format(
' DDRVBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
924 $ i6,
', N=', i6,
', JTYPE=', i6,
', LSWORK=', i6, / 9x,
925 $
'ISEED=(', 3( i5,
',' ), i5,
')' )