434 SUBROUTINE dchkbd( 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,
447 DOUBLE PRECISION THRESH
451 INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * )
452 DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE, TWO, HALF
462 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
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 DOUBLE PRECISION 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 DOUBLE PRECISION DUM( 1 ), DUMMA( 1 ), RESULT( 19 )
483 DOUBLE PRECISION DLAMCH, DLARND
484 EXTERNAL dlamch, dlarnd
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(
'DCHKBD', -info )
565 path( 1: 1 ) =
'Double precision'
569 unfl = dlamch(
'Safe minimum' )
570 ovfl = dlamch(
'Overflow' )
572 ulp = dlamch(
'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 dlaset(
'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 dlatms( 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 dlatms( 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 dlatms( 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 dlatmr( 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 dlatmr( 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 dlatmr( 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*dlarnd( 2, iseed ) )
725 $ be( j ) = exp( temp1*dlarnd( 2, iseed ) )
739 IF( iinfo.EQ.0 )
THEN
744 CALL dlatmr( 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 dlatmr( 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 dlacpy(
' ', m, n, a, lda, q, ldq )
777 CALL dgebrd( 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 )
'DGEBRD', iinfo, m, n,
789 CALL dlacpy(
' ', m, n, q, ldq, pt, ldpt )
801 CALL dorgbr(
'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 )
'DORGBR(Q)', iinfo, m, n,
815 CALL dorgbr(
'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 )
'DORGBR(P)', iinfo, m, n,
829 CALL dgemm(
'Transpose',
'No transpose', m, nrhs, m, one,
830 $ q, ldq, x, ldx, zero, y, ldx )
836 CALL dbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
837 $ work, result( 1 ) )
838 CALL dort01(
'Columns', m, mq, q, ldq, work, lwork,
840 CALL dort01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
847 CALL dcopy( mnmin, bd, 1, s1, 1 )
849 $
CALL dcopy( mnmin-1, be, 1, work, 1 )
850 CALL dlacpy(
' ', m, nrhs, y, ldx, z, ldx )
851 CALL dlaset(
'Full', mnmin, mnmin, zero, one, u, ldpt )
852 CALL dlaset(
'Full', mnmin, mnmin, zero, one, vt, ldpt )
854 CALL dbdsqr( 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 )
'DBDSQR(vects)', iinfo, m, n,
863 IF( iinfo.LT.0 )
THEN
874 CALL dcopy( mnmin, bd, 1, s2, 1 )
876 $
CALL dcopy( mnmin-1, be, 1, work, 1 )
878 CALL dbdsqr( 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 )
'DBDSQR(values)', iinfo, m, n,
887 IF( iinfo.LT.0 )
THEN
900 CALL dbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
901 $ work, result( 4 ) )
902 CALL dbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
904 CALL dort01(
'Columns', mnmin, mnmin, u, ldpt, work, lwork,
906 CALL dort01(
'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 dcopy( mnmin, bd, 1, s2, 1 )
958 $
CALL dcopy( mnmin-1, be, 1, work, 1 )
960 CALL dbdsqr( uplo, mnmin, n, m, nrhs, s2, work, pt, ldpt,
961 $ q, ldq, y, ldx, work( mnmin+1 ), iinfo )
968 CALL dbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
969 $ ldpt, work, result( 11 ) )
970 CALL dbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
972 CALL dort01(
'Columns', m, mq, q, ldq, work, lwork,
974 CALL dort01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
981 CALL dcopy( mnmin, bd, 1, s1, 1 )
983 $
CALL dcopy( mnmin-1, be, 1, work, 1 )
984 CALL dlaset(
'Full', mnmin, mnmin, zero, one, u, ldpt )
985 CALL dlaset(
'Full', mnmin, mnmin, zero, one, vt, ldpt )
987 CALL dbdsdc( 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 )
'DBDSDC(vects)', iinfo, m, n,
996 IF( iinfo.LT.0 )
THEN
999 result( 15 ) = ulpinv
1007 CALL dcopy( mnmin, bd, 1, s2, 1 )
1009 $
CALL dcopy( mnmin-1, be, 1, work, 1 )
1011 CALL dbdsdc( 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 )
'DBDSDC(values)', iinfo, m, n,
1020 IF( iinfo.LT.0 )
THEN
1023 result( 18 ) = ulpinv
1032 CALL dbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
1033 $ work, result( 15 ) )
1034 CALL dort01(
'Columns', mnmin, mnmin, u, ldpt, work, lwork,
1036 CALL dort01(
'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 dlahd2( 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(
' DCHKBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
1099 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
subroutine dbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
DBDSDC
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
subroutine dbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
DBDT01
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dlatmr(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)
DLATMR
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
subroutine dbdt02(M, N, B, LDB, C, LDC, U, LDU, WORK, RESID)
DBDT02
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
subroutine dlahd2(IOUNIT, PATH)
DLAHD2
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 dbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
DBDT03
subroutine dchkbd(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)
DCHKBD
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR
subroutine dort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
DORT01