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,
')' )
subroutine dort03(RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK, RESULT, INFO)
DORT03
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
subroutine dgejsv(JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, WORK, LWORK, IWORK, INFO)
DGEJSV
subroutine dgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
DGESVJ
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 dgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
DGESDD
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine ddrvbd(NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH, A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S, SSAV, E, WORK, LWORK, IWORK, NOUT, INFO)
DDRVBD
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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 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 dort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
DORT01