434 SUBROUTINE schkbd( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
435 $ iseed, thresh, a, lda, bd, be, s1, s2, x, ldx,
436 $ y, z, q, ldq, pt, ldpt, u, vt, work, lwork,
437 $ iwork, nout, info )
445 INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
451 INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * )
452 REAL A( lda, * ), BD( * ), BE( * ), PT( ldpt, * ),
453 $ q( ldq, * ), s1( * ), s2( * ), u( ldpt, * ),
454 $ vt( ldpt, * ), work( * ), x( ldx, * ),
455 $ y( ldx, * ), z( ldx, * )
461 REAL ZERO, ONE, TWO, HALF
462 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
465 parameter( maxtyp = 16 )
468 LOGICAL BADMM, BADNN, BIDIAG
471 INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JSIZE, JTYPE,
472 $ log2ui, m, minwrk, mmax, mnmax, mnmin, mq,
473 $ mtypes, n, nfail, nmax, ntest
474 REAL AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
475 $ temp1, temp2, ulp, ulpinv, unfl
478 INTEGER IDUM( 1 ), IOLDSD( 4 ), KMAGN( maxtyp ),
479 $ kmode( maxtyp ), ktype( maxtyp )
480 REAL DUM( 1 ), DUMMA( 1 ), RESULT( 19 )
484 EXTERNAL slamch, slarnd
492 INTRINSIC abs, exp, int, log, max, min, sqrt
500 COMMON / infoc / infot, nunit, ok, lerr
501 COMMON / srnamc / srnamt
504 DATA ktype / 1, 2, 5*4, 5*6, 3*9, 10 /
505 DATA kmagn / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
506 DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
522 mmax = max( mmax, mval( j ) )
525 nmax = max( nmax, nval( j ) )
528 mnmax = max( mnmax, min( mval( j ), nval( j ) ) )
529 minwrk = max( minwrk, 3*( mval( j )+nval( j ) ),
530 $ mval( j )*( mval( j )+max( mval( j ), nval( j ),
531 $ nrhs )+1 )+nval( j )*min( nval( j ), mval( j ) ) )
536 IF( nsizes.LT.0 )
THEN
538 ELSE IF( badmm )
THEN
540 ELSE IF( badnn )
THEN
542 ELSE IF( ntypes.LT.0 )
THEN
544 ELSE IF( nrhs.LT.0 )
THEN
546 ELSE IF( lda.LT.mmax )
THEN
548 ELSE IF( ldx.LT.mmax )
THEN
550 ELSE IF( ldq.LT.mmax )
THEN
552 ELSE IF( ldpt.LT.mnmax )
THEN
554 ELSE IF( minwrk.GT.lwork )
THEN
559 CALL xerbla(
'SCHKBD', -info )
565 path( 1: 1 ) =
'Single precision'
569 unfl = slamch(
'Safe minimum' )
570 ovfl = slamch(
'Overflow' )
572 ulp = slamch(
'Precision' )
574 log2ui = int( log( ulpinv ) / log( two ) )
575 rtunfl = sqrt( unfl )
576 rtovfl = sqrt( ovfl )
581 DO 200 jsize = 1, nsizes
585 amninv = one / max( m, n, 1 )
587 IF( nsizes.NE.1 )
THEN
588 mtypes = min( maxtyp, ntypes )
590 mtypes = min( maxtyp+1, ntypes )
593 DO 190 jtype = 1, mtypes
594 IF( .NOT.dotype( jtype ) )
598 ioldsd( j ) = iseed( j )
623 IF( mtypes.GT.maxtyp )
626 itype = ktype( jtype )
627 imode = kmode( jtype )
631 GO TO ( 40, 50, 60 )kmagn( jtype )
638 anorm = ( rtovfl*ulp )*amninv
642 anorm = rtunfl*max( m, n )*ulpinv
647 CALL slaset(
'Full', lda, n, zero, zero, a, lda )
652 IF( itype.EQ.1 )
THEN
658 ELSE IF( itype.EQ.2 )
THEN
662 DO 80 jcol = 1, mnmin
663 a( jcol, jcol ) = anorm
666 ELSE IF( itype.EQ.4 )
THEN
670 CALL slatms( mnmin, mnmin,
'S', iseed,
'N', work, imode,
671 $ cond, anorm, 0, 0,
'N', a, lda,
672 $ work( mnmin+1 ), iinfo )
674 ELSE IF( itype.EQ.5 )
THEN
678 CALL slatms( mnmin, mnmin,
'S', iseed,
'S', work, imode,
679 $ cond, anorm, m, n,
'N', a, lda,
680 $ work( mnmin+1 ), iinfo )
682 ELSE IF( itype.EQ.6 )
THEN
686 CALL slatms( m, n,
'S', iseed,
'N', work, imode, cond,
687 $ anorm, m, n,
'N', a, lda, work( mnmin+1 ),
690 ELSE IF( itype.EQ.7 )
THEN
694 CALL slatmr( mnmin, mnmin,
'S', iseed,
'N', work, 6, one,
695 $ one,
'T',
'N', work( mnmin+1 ), 1, one,
696 $ work( 2*mnmin+1 ), 1, one,
'N', iwork, 0, 0,
697 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
699 ELSE IF( itype.EQ.8 )
THEN
703 CALL slatmr( mnmin, mnmin,
'S', iseed,
'S', work, 6, one,
704 $ one,
'T',
'N', work( mnmin+1 ), 1, one,
705 $ work( m+mnmin+1 ), 1, one,
'N', iwork, m, n,
706 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
708 ELSE IF( itype.EQ.9 )
THEN
712 CALL slatmr( m, n,
'S', iseed,
'N', work, 6, one, one,
713 $
'T',
'N', work( mnmin+1 ), 1, one,
714 $ work( m+mnmin+1 ), 1, one,
'N', iwork, m, n,
715 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
717 ELSE IF( itype.EQ.10 )
THEN
721 temp1 = -two*log( ulp )
723 bd( j ) = exp( temp1*slarnd( 2, iseed ) )
725 $ be( j ) = exp( temp1*slarnd( 2, iseed ) )
739 IF( iinfo.EQ.0 )
THEN
744 CALL slatmr( mnmin, nrhs,
'S', iseed,
'N', work, 6,
745 $ one, one,
'T',
'N', work( mnmin+1 ), 1,
746 $ one, work( 2*mnmin+1 ), 1, one,
'N',
747 $ iwork, mnmin, nrhs, zero, one,
'NO', y,
748 $ ldx, iwork, iinfo )
750 CALL slatmr( m, nrhs,
'S', iseed,
'N', work, 6, one,
751 $ one,
'T',
'N', work( m+1 ), 1, one,
752 $ work( 2*m+1 ), 1, one,
'N', iwork, m,
753 $ nrhs, zero, one,
'NO', x, ldx, iwork,
760 IF( iinfo.NE.0 )
THEN
761 WRITE( nout, fmt = 9998 )
'Generator', iinfo, m, n,
771 IF( .NOT.bidiag )
THEN
776 CALL slacpy(
' ', m, n, a, lda, q, ldq )
777 CALL sgebrd( m, n, q, ldq, bd, be, work, work( mnmin+1 ),
778 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
782 IF( iinfo.NE.0 )
THEN
783 WRITE( nout, fmt = 9998 )
'SGEBRD', iinfo, m, n,
789 CALL slacpy(
' ', m, n, q, ldq, pt, ldpt )
801 CALL sorgbr(
'Q', m, mq, n, q, ldq, work,
802 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
806 IF( iinfo.NE.0 )
THEN
807 WRITE( nout, fmt = 9998 )
'SORGBR(Q)', iinfo, m, n,
815 CALL sorgbr(
'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
816 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
820 IF( iinfo.NE.0 )
THEN
821 WRITE( nout, fmt = 9998 )
'SORGBR(P)', iinfo, m, n,
829 CALL sgemm(
'Transpose',
'No transpose', m, nrhs, m, one,
830 $ q, ldq, x, ldx, zero, y, ldx )
836 CALL sbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
837 $ work, result( 1 ) )
838 CALL sort01(
'Columns', m, mq, q, ldq, work, lwork,
840 CALL sort01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
847 CALL scopy( mnmin, bd, 1, s1, 1 )
849 $
CALL scopy( mnmin-1, be, 1, work, 1 )
850 CALL slacpy(
' ', m, nrhs, y, ldx, z, ldx )
851 CALL slaset(
'Full', mnmin, mnmin, zero, one, u, ldpt )
852 CALL slaset(
'Full', mnmin, mnmin, zero, one, vt, ldpt )
854 CALL sbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, work, vt,
855 $ ldpt, u, ldpt, z, ldx, work( mnmin+1 ), iinfo )
859 IF( iinfo.NE.0 )
THEN
860 WRITE( nout, fmt = 9998 )
'SBDSQR(vects)', iinfo, m, n,
863 IF( iinfo.LT.0 )
THEN
874 CALL scopy( mnmin, bd, 1, s2, 1 )
876 $
CALL scopy( mnmin-1, be, 1, work, 1 )
878 CALL sbdsqr( uplo, mnmin, 0, 0, 0, s2, work, vt, ldpt, u,
879 $ ldpt, z, ldx, work( mnmin+1 ), iinfo )
883 IF( iinfo.NE.0 )
THEN
884 WRITE( nout, fmt = 9998 )
'SBDSQR(values)', iinfo, m, n,
887 IF( iinfo.LT.0 )
THEN
900 CALL sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
901 $ work, result( 4 ) )
902 CALL sbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
904 CALL sort01(
'Columns', mnmin, mnmin, u, ldpt, work, lwork,
906 CALL sort01(
'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
913 DO 110 i = 1, mnmin - 1
914 IF( s1( i ).LT.s1( i+1 ) )
915 $ result( 8 ) = ulpinv
916 IF( s1( i ).LT.zero )
917 $ result( 8 ) = ulpinv
919 IF( mnmin.GE.1 )
THEN
920 IF( s1( mnmin ).LT.zero )
921 $ result( 8 ) = ulpinv
929 temp1 = abs( s1( j )-s2( j ) ) /
930 $ max( sqrt( unfl )*max( s1( 1 ), one ),
931 $ ulp*max( abs( s1( j ) ), abs( s2( j ) ) ) )
932 temp2 = max( temp1, temp2 )
940 temp1 = thresh*( half-ulp )
955 IF( .NOT.bidiag )
THEN
956 CALL scopy( mnmin, bd, 1, s2, 1 )
958 $
CALL scopy( mnmin-1, be, 1, work, 1 )
960 CALL sbdsqr( uplo, mnmin, n, m, nrhs, s2, work, pt, ldpt,
961 $ q, ldq, y, ldx, work( mnmin+1 ), iinfo )
968 CALL sbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
969 $ ldpt, work, result( 11 ) )
970 CALL sbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
972 CALL sort01(
'Columns', m, mq, q, ldq, work, lwork,
974 CALL sort01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
981 CALL scopy( mnmin, bd, 1, s1, 1 )
983 $
CALL scopy( mnmin-1, be, 1, work, 1 )
984 CALL slaset(
'Full', mnmin, mnmin, zero, one, u, ldpt )
985 CALL slaset(
'Full', mnmin, mnmin, zero, one, vt, ldpt )
987 CALL sbdsdc( uplo,
'I', mnmin, s1, work, u, ldpt, vt, ldpt,
988 $ dum, idum, work( mnmin+1 ), iwork, iinfo )
992 IF( iinfo.NE.0 )
THEN
993 WRITE( nout, fmt = 9998 )
'SBDSDC(vects)', iinfo, m, n,
996 IF( iinfo.LT.0 )
THEN
999 result( 15 ) = ulpinv
1007 CALL scopy( mnmin, bd, 1, s2, 1 )
1009 $
CALL scopy( mnmin-1, be, 1, work, 1 )
1011 CALL sbdsdc( uplo,
'N', mnmin, s2, work, dum, 1, dum, 1,
1012 $ dum, idum, work( mnmin+1 ), iwork, iinfo )
1016 IF( iinfo.NE.0 )
THEN
1017 WRITE( nout, fmt = 9998 )
'SBDSDC(values)', iinfo, m, n,
1020 IF( iinfo.LT.0 )
THEN
1023 result( 18 ) = ulpinv
1032 CALL sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
1033 $ work, result( 15 ) )
1034 CALL sort01(
'Columns', mnmin, mnmin, u, ldpt, work, lwork,
1036 CALL sort01(
'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
1043 DO 150 i = 1, mnmin - 1
1044 IF( s1( i ).LT.s1( i+1 ) )
1045 $ result( 18 ) = ulpinv
1046 IF( s1( i ).LT.zero )
1047 $ result( 18 ) = ulpinv
1049 IF( mnmin.GE.1 )
THEN
1050 IF( s1( mnmin ).LT.zero )
1051 $ result( 18 ) = ulpinv
1059 temp1 = abs( s1( j )-s2( j ) ) /
1060 $ max( sqrt( unfl )*max( s1( 1 ), one ),
1061 $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1062 temp2 = max( temp1, temp2 )
1065 result( 19 ) = temp2
1071 IF( result( j ).GE.thresh )
THEN
1073 $
CALL slahd2( nout, path )
1074 WRITE( nout, fmt = 9999 )m, n, jtype, ioldsd, j,
1079 IF( .NOT.bidiag )
THEN
1090 CALL alasum( path, nout, nfail, ntest, 0 )
1096 9999
FORMAT(
' M=', i5,
', N=', i5,
', type ', i2,
', seed=',
1097 $ 4( i4,
',' ),
' test(', i2,
')=', g11.4 )
1098 9998
FORMAT(
' SCHKBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
1099 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
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 sbdt02(M, N, B, LDB, C, LDC, U, LDU, WORK, RESID)
SBDT02
subroutine schkbd(NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS, ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX, Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK, IWORK, NOUT, INFO)
SCHKBD
subroutine sbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
SBDT03
subroutine slahd2(IOUNIT, PATH)
SLAHD2
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
subroutine slatmr(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RSIGN, GRADE, DL, MODEL, CONDL, DR, MODER, CONDR, PIVTNG, IPIVOT, KL, KU, SPARSE, ANORM, PACK, A, LDA, IWORK, INFO)
SLATMR
subroutine sorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGBR
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
subroutine sbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
SBDSDC
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
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
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine sbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
SBDT01