318 SUBROUTINE sdrvbd( 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,
334 INTEGER iseed( 4 ), iwork( * ), mm( * ), nn( * )
335 REAL a( lda, * ), asav( lda, * ), e( * ), s( * ),
336 $ ssav( * ), u( ldu, * ), usav( ldu, * ),
337 $ vt( ldvt, * ), vtsav( ldvt, * ), work( * )
344 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL anorm, dif, div, ovfl, ulp, ulpinv, unfl
372 INTRINSIC abs, max, min, real
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(
'SDRVBD', -info )
437 path( 1: 1 ) =
'Single precision'
441 unfl =
slamch(
'Safe minimum' )
444 ulp =
slamch(
'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
slaset(
'Full', m, n, zero, zero, a, lda )
480 ELSE IF( jtype.EQ.2 )
THEN
484 CALL
slaset(
'Full', m, n, zero, one, a, lda )
496 CALL
slatms( m, n,
'U', iseed,
'N', s, 4,
REAL( 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
slacpy(
'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
slacpy(
'F', m, n, asav, lda, a, lda )
529 CALL
sgesvd(
'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
sbdt01( 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
sort01(
'Columns', m, m, usav, ldu, work, lwork,
545 CALL
sort01(
'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
slacpy(
'F', m, n, asav, lda, a, lda )
573 CALL
sgesvd( 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
sort03(
'C', m, mnmin, m, mnmin, usav,
582 $ ldu, a, lda, work, lwork, dif,
584 ELSE IF( iju.EQ.2 )
THEN
585 CALL
sort03(
'C', m, mnmin, m, mnmin, usav,
586 $ ldu, u, ldu, work, lwork, dif,
588 ELSE IF( iju.EQ.3 )
THEN
589 CALL
sort03(
'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
sort03(
'R', n, mnmin, n, mnmin, vtsav,
602 $ ldvt, a, lda, work, lwork, dif,
604 ELSE IF( ijvt.EQ.2 )
THEN
605 CALL
sort03(
'R', n, mnmin, n, mnmin, vtsav,
606 $ ldvt, vt, ldvt, work, lwork,
608 ELSE IF( ijvt.EQ.3 )
THEN
609 CALL
sort03(
'R', n, n, n, mnmin, vtsav,
610 $ ldvt, vt, ldvt, work, lwork,
614 result( 6 ) = max( result( 6 ), dif )
619 div = max(
REAL( 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
slacpy(
'F', m, n, asav, lda, a, lda )
642 CALL
sgesdd(
'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
sbdt01( 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
sort01(
'Columns', m, m, usav, ldu, work, lwork,
658 CALL
sort01(
'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
slacpy(
'F', m, n, asav, lda, a, lda )
682 CALL
sgesdd( 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
sort03(
'C', m, mnmin, m, mnmin, usav,
692 $ ldu, a, lda, work, lwork, dif,
695 CALL
sort03(
'C', m, mnmin, m, mnmin, usav,
696 $ ldu, u, ldu, work, lwork, dif,
699 ELSE IF( ijq.EQ.2 )
THEN
700 CALL
sort03(
'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
sort03(
'R', n, mnmin, n, mnmin, vtsav,
713 $ ldvt, vt, ldvt, work, lwork,
716 CALL
sort03(
'R', n, mnmin, n, mnmin, vtsav,
717 $ ldvt, a, lda, work, lwork, dif,
720 ELSE IF( ijq.EQ.2 )
THEN
721 CALL
sort03(
'R', n, mnmin, n, mnmin, vtsav,
722 $ ldvt, vt, ldvt, work, lwork, dif,
726 result( 13 ) = max( result( 13 ), dif )
731 div = max(
REAL( 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
slacpy(
'F', m, n, asav, lda, usav, lda )
760 CALL
sgesvj(
'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
sbdt01( 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
sort01(
'Columns', m, m, usav, ldu, work,
785 $ lwork, result( 16 ) )
786 CALL
sort01(
'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
slacpy(
'F', m, n, asav, lda, vtsav, lda )
819 CALL
sgejsv(
'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
sbdt01( 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
sort01(
'Columns', m, m, usav, ldu, work,
845 $ lwork, result( 20 ) )
846 CALL
sort01(
'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 SDRVBD 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 $ / 19
x,
' S is an array, and Upartial, VTpartial, and',
893 $ / 19
x,
' 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(
' SDRVBD: ', a,
' returned INFO=', i6,
'.', / 9
x,
'M=',
921 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
923 9995
FORMAT(
' SDRVBD: ', a,
' returned INFO=', i6,
'.', / 9
x,
'M=',
924 $ i6,
', N=', i6,
', JTYPE=', i6,
', LSWORK=', i6, / 9
x,
925 $
'ISEED=(', 3( i5,
',' ), i5,
')' )
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine sdrvbd(NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH, A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S, SSAV, E, WORK, LWORK, IWORK, NOUT, INFO)
SDRVBD
REAL function slamch(CMACH)
SLAMCH
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sgejsv(JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, WORK, LWORK, IWORK, INFO)
SGEJSV
subroutine sort03(RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK, RESULT, INFO)
SORT03
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine sgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
SGESDD
subroutine sgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
SGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine sgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
SGESVJ
subroutine sbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
SBDT01