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 )
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,
',' ),