411 SUBROUTINE cchkbd( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
412 $ ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX,
413 $ Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK,
414 $ RWORK, NOUT, INFO )
421 INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
427 INTEGER ISEED( 4 ), MVAL( * ), NVAL( * )
428 REAL BD( * ), BE( * ), RWORK( * ), S1( * ), S2( * )
429 COMPLEX A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
430 $ u( ldpt, * ), vt( ldpt, * ), work( * ),
431 $ x( ldx, * ), y( ldx, * ), z( ldx, * )
437 REAL ZERO, ONE, TWO, HALF
438 PARAMETER ( ZERO = 0.0e0, one = 1.0e0, two = 2.0e0,
441 parameter( czero = ( 0.0e+0, 0.0e+0 ),
442 $ cone = ( 1.0e+0, 0.0e+0 ) )
444 parameter( maxtyp = 16 )
447 LOGICAL BADMM, BADNN, BIDIAG
450 INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JSIZE, JTYPE,
451 $ log2ui, m, minwrk, mmax, mnmax, mnmin, mq,
452 $ mtypes, n, nfail, nmax, ntest
453 REAL AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
454 $ TEMP1, TEMP2, ULP, ULPINV, UNFL
457 INTEGER IOLDSD( 4 ), IWORK( 1 ), KMAGN( MAXTYP ),
458 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
459 REAL DUMMA( 1 ), RESULT( 14 )
463 EXTERNAL SLAMCH, SLARND
472 INTRINSIC abs, exp, int, log, max, min, sqrt
480 COMMON / infoc / infot, nunit, ok, lerr
481 COMMON / srnamc / srnamt
484 DATA ktype / 1, 2, 5*4, 5*6, 3*9, 10 /
485 DATA kmagn / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
486 DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
502 mmax = max( mmax, mval( j ) )
505 nmax = max( nmax, nval( j ) )
508 mnmax = max( mnmax, min( mval( j ), nval( j ) ) )
509 minwrk = max( minwrk, 3*( mval( j )+nval( j ) ),
510 $ mval( j )*( mval( j )+max( mval( j ), nval( j ),
511 $ nrhs )+1 )+nval( j )*min( nval( j ), mval( j ) ) )
516 IF( nsizes.LT.0 )
THEN
518 ELSE IF( badmm )
THEN
520 ELSE IF( badnn )
THEN
522 ELSE IF( ntypes.LT.0 )
THEN
524 ELSE IF( nrhs.LT.0 )
THEN
526 ELSE IF( lda.LT.mmax )
THEN
528 ELSE IF( ldx.LT.mmax )
THEN
530 ELSE IF( ldq.LT.mmax )
THEN
532 ELSE IF( ldpt.LT.mnmax )
THEN
534 ELSE IF( minwrk.GT.lwork )
THEN
539 CALL xerbla(
'CCHKBD', -info )
545 path( 1: 1 ) =
'Complex precision'
549 unfl = slamch(
'Safe minimum' )
550 ovfl = slamch(
'Overflow' )
551 ulp = slamch(
'Precision' )
553 log2ui = int( log( ulpinv ) / log( two ) )
554 rtunfl = sqrt( unfl )
555 rtovfl = sqrt( ovfl )
560 DO 180 jsize = 1, nsizes
564 amninv = one / max( m, n, 1 )
566 IF( nsizes.NE.1 )
THEN
567 mtypes = min( maxtyp, ntypes )
569 mtypes = min( maxtyp+1, ntypes )
572 DO 170 jtype = 1, mtypes
573 IF( .NOT.dotype( jtype ) )
577 ioldsd( j ) = iseed( j )
602 IF( mtypes.GT.maxtyp )
605 itype = ktype( jtype )
606 imode = kmode( jtype )
610 GO TO ( 40, 50, 60 )kmagn( jtype )
617 anorm = ( rtovfl*ulp )*amninv
621 anorm = rtunfl*max( m, n )*ulpinv
626 CALL claset(
'Full', lda, n, czero, czero, a, lda )
631 IF( itype.EQ.1 )
THEN
637 ELSE IF( itype.EQ.2 )
THEN
641 DO 80 jcol = 1, mnmin
642 a( jcol, jcol ) = anorm
645 ELSE IF( itype.EQ.4 )
THEN
649 CALL clatms( mnmin, mnmin,
'S', iseed,
'N', rwork, imode,
650 $ cond, anorm, 0, 0,
'N', a, lda, work,
653 ELSE IF( itype.EQ.5 )
THEN
657 CALL clatms( mnmin, mnmin,
'S', iseed,
'S', rwork, imode,
658 $ cond, anorm, m, n,
'N', a, lda, work,
661 ELSE IF( itype.EQ.6 )
THEN
665 CALL clatms( m, n,
'S', iseed,
'N', rwork, imode, cond,
666 $ anorm, m, n,
'N', a, lda, work, iinfo )
668 ELSE IF( itype.EQ.7 )
THEN
672 CALL clatmr( mnmin, mnmin,
'S', iseed,
'N', work, 6, one,
673 $ cone,
'T',
'N', work( mnmin+1 ), 1, one,
674 $ work( 2*mnmin+1 ), 1, one,
'N', iwork, 0, 0,
675 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
677 ELSE IF( itype.EQ.8 )
THEN
681 CALL clatmr( mnmin, mnmin,
'S', iseed,
'S', work, 6, one,
682 $ cone,
'T',
'N', work( mnmin+1 ), 1, one,
683 $ work( m+mnmin+1 ), 1, one,
'N', iwork, m, n,
684 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
686 ELSE IF( itype.EQ.9 )
THEN
690 CALL clatmr( m, n,
'S', iseed,
'N', work, 6, one, cone,
691 $
'T',
'N', work( mnmin+1 ), 1, one,
692 $ work( m+mnmin+1 ), 1, one,
'N', iwork, m, n,
693 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
695 ELSE IF( itype.EQ.10 )
THEN
699 temp1 = -two*log( ulp )
701 bd( j ) = exp( temp1*slarnd( 2, iseed ) )
703 $ be( j ) = exp( temp1*slarnd( 2, iseed ) )
717 IF( iinfo.EQ.0 )
THEN
722 CALL clatmr( mnmin, nrhs,
'S', iseed,
'N', work, 6,
723 $ one, cone,
'T',
'N', work( mnmin+1 ), 1,
724 $ one, work( 2*mnmin+1 ), 1, one,
'N',
725 $ iwork, mnmin, nrhs, zero, one,
'NO', y,
726 $ ldx, iwork, iinfo )
728 CALL clatmr( m, nrhs,
'S', iseed,
'N', work, 6, one,
729 $ cone,
'T',
'N', work( m+1 ), 1, one,
730 $ work( 2*m+1 ), 1, one,
'N', iwork, m,
731 $ nrhs, zero, one,
'NO', x, ldx, iwork,
738 IF( iinfo.NE.0 )
THEN
739 WRITE( nout, fmt = 9998 )
'Generator', iinfo, m, n,
749 IF( .NOT.bidiag )
THEN
754 CALL clacpy(
' ', m, n, a, lda, q, ldq )
755 CALL cgebrd( m, n, q, ldq, bd, be, work, work( mnmin+1 ),
756 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
760 IF( iinfo.NE.0 )
THEN
761 WRITE( nout, fmt = 9998 )
'CGEBRD', iinfo, m, n,
767 CALL clacpy(
' ', m, n, q, ldq, pt, ldpt )
779 CALL cungbr(
'Q', m, mq, n, q, ldq, work,
780 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
784 IF( iinfo.NE.0 )
THEN
785 WRITE( nout, fmt = 9998 )
'CUNGBR(Q)', iinfo, m, n,
793 CALL cungbr(
'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
794 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
798 IF( iinfo.NE.0 )
THEN
799 WRITE( nout, fmt = 9998 )
'CUNGBR(P)', iinfo, m, n,
807 CALL cgemm(
'Conjugate transpose',
'No transpose', m,
808 $ nrhs, m, cone, q, ldq, x, ldx, czero, y,
815 CALL cbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
816 $ work, rwork, result( 1 ) )
817 CALL cunt01(
'Columns', m, mq, q, ldq, work, lwork,
818 $ rwork, result( 2 ) )
819 CALL cunt01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
820 $ rwork, result( 3 ) )
826 CALL scopy( mnmin, bd, 1, s1, 1 )
828 $
CALL scopy( mnmin-1, be, 1, rwork, 1 )
829 CALL clacpy(
' ', m, nrhs, y, ldx, z, ldx )
830 CALL claset(
'Full', mnmin, mnmin, czero, cone, u, ldpt )
831 CALL claset(
'Full', mnmin, mnmin, czero, cone, vt, ldpt )
833 CALL cbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, rwork, vt,
834 $ ldpt, u, ldpt, z, ldx, rwork( mnmin+1 ),
839 IF( iinfo.NE.0 )
THEN
840 WRITE( nout, fmt = 9998 )
'CBDSQR(vects)', iinfo, m, n,
843 IF( iinfo.LT.0 )
THEN
854 CALL scopy( mnmin, bd, 1, s2, 1 )
856 $
CALL scopy( mnmin-1, be, 1, rwork, 1 )
858 CALL cbdsqr( uplo, mnmin, 0, 0, 0, s2, rwork, vt, ldpt, u,
859 $ ldpt, z, ldx, rwork( mnmin+1 ), iinfo )
863 IF( iinfo.NE.0 )
THEN
864 WRITE( nout, fmt = 9998 )
'CBDSQR(values)', iinfo, m, n,
867 IF( iinfo.LT.0 )
THEN
880 CALL cbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
881 $ work, result( 4 ) )
882 CALL cbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
883 $ rwork, result( 5 ) )
884 CALL cunt01(
'Columns', mnmin, mnmin, u, ldpt, work, lwork,
885 $ rwork, result( 6 ) )
886 CALL cunt01(
'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
887 $ rwork, result( 7 ) )
893 DO 110 i = 1, mnmin - 1
894 IF( s1( i ).LT.s1( i+1 ) )
895 $ result( 8 ) = ulpinv
896 IF( s1( i ).LT.zero )
897 $ result( 8 ) = ulpinv
899 IF( mnmin.GE.1 )
THEN
900 IF( s1( mnmin ).LT.zero )
901 $ result( 8 ) = ulpinv
909 temp1 = abs( s1( j )-s2( j ) ) /
910 $ max( sqrt( unfl )*max( s1( 1 ), one ),
911 $ ulp*max( abs( s1( j ) ), abs( s2( j ) ) ) )
912 temp2 = max( temp1, temp2 )
920 temp1 = thresh*( half-ulp )
923 CALL ssvdch( mnmin, bd, be, s1, temp1, iinfo )
935 IF( .NOT.bidiag )
THEN
936 CALL scopy( mnmin, bd, 1, s2, 1 )
938 $
CALL scopy( mnmin-1, be, 1, rwork, 1 )
940 CALL cbdsqr( uplo, mnmin, n, m, nrhs, s2, rwork, pt,
941 $ ldpt, q, ldq, y, ldx, rwork( mnmin+1 ),
949 CALL cbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
950 $ ldpt, work, rwork, result( 11 ) )
951 CALL cbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
952 $ rwork, result( 12 ) )
953 CALL cunt01(
'Columns', m, mq, q, ldq, work, lwork,
954 $ rwork, result( 13 ) )
955 CALL cunt01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
956 $ rwork, result( 14 ) )
963 IF( result( j ).GE.thresh )
THEN
965 $
CALL slahd2( nout, path )
966 WRITE( nout, fmt = 9999 )m, n, jtype, ioldsd, j,
971 IF( .NOT.bidiag )
THEN
982 CALL alasum( path, nout, nfail, ntest, 0 )
988 9999
FORMAT(
' M=', i5,
', N=', i5,
', type ', i2,
', seed=',
989 $ 4( i4,
',' ),
' test(', i2,
')=', g11.4 )
990 9998
FORMAT(
' CCHKBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
991 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),