413 SUBROUTINE cchkbd( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
414 $ iseed, thresh, a, lda, bd, be, s1, s2, x, ldx,
415 $ y, z, q, ldq, pt, ldpt, u, vt, work, lwork,
416 $ rwork, nout, info )
424 INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
430 INTEGER ISEED( 4 ), MVAL( * ), NVAL( * )
431 REAL BD( * ), BE( * ), RWORK( * ), S1( * ), S2( * )
432 COMPLEX A( lda, * ), PT( ldpt, * ), Q( ldq, * ),
433 $ u( ldpt, * ), vt( ldpt, * ), work( * ),
434 $ x( ldx, * ), y( ldx, * ), z( ldx, * )
440 REAL ZERO, ONE, TWO, HALF
441 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
444 parameter( czero = ( 0.0e+0, 0.0e+0 ),
445 $ cone = ( 1.0e+0, 0.0e+0 ) )
447 parameter( maxtyp = 16 )
450 LOGICAL BADMM, BADNN, BIDIAG
453 INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JSIZE, JTYPE,
454 $ log2ui, m, minwrk, mmax, mnmax, mnmin, mq,
455 $ mtypes, n, nfail, nmax, ntest
456 REAL AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
457 $ temp1, temp2, ulp, ulpinv, unfl
460 INTEGER IOLDSD( 4 ), IWORK( 1 ), KMAGN( maxtyp ),
461 $ kmode( maxtyp ), ktype( maxtyp )
462 REAL DUMMA( 1 ), RESULT( 14 )
466 EXTERNAL slamch, slarnd
474 INTRINSIC abs, exp, int, log, max, min, sqrt
482 COMMON / infoc / infot, nunit, ok, lerr
483 COMMON / srnamc / srnamt
486 DATA ktype / 1, 2, 5*4, 5*6, 3*9, 10 /
487 DATA kmagn / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
488 DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
504 mmax = max( mmax, mval( j ) )
507 nmax = max( nmax, nval( j ) )
510 mnmax = max( mnmax, min( mval( j ), nval( j ) ) )
511 minwrk = max( minwrk, 3*( mval( j )+nval( j ) ),
512 $ mval( j )*( mval( j )+max( mval( j ), nval( j ),
513 $ nrhs )+1 )+nval( j )*min( nval( j ), mval( j ) ) )
518 IF( nsizes.LT.0 )
THEN
520 ELSE IF( badmm )
THEN
522 ELSE IF( badnn )
THEN
524 ELSE IF( ntypes.LT.0 )
THEN
526 ELSE IF( nrhs.LT.0 )
THEN
528 ELSE IF( lda.LT.mmax )
THEN
530 ELSE IF( ldx.LT.mmax )
THEN
532 ELSE IF( ldq.LT.mmax )
THEN
534 ELSE IF( ldpt.LT.mnmax )
THEN
536 ELSE IF( minwrk.GT.lwork )
THEN
541 CALL xerbla(
'CCHKBD', -info )
547 path( 1: 1 ) =
'Complex precision'
551 unfl = slamch(
'Safe minimum' )
552 ovfl = slamch(
'Overflow' )
554 ulp = slamch(
'Precision' )
556 log2ui = int( log( ulpinv ) / log( two ) )
557 rtunfl = sqrt( unfl )
558 rtovfl = sqrt( ovfl )
563 DO 180 jsize = 1, nsizes
567 amninv = one / max( m, n, 1 )
569 IF( nsizes.NE.1 )
THEN
570 mtypes = min( maxtyp, ntypes )
572 mtypes = min( maxtyp+1, ntypes )
575 DO 170 jtype = 1, mtypes
576 IF( .NOT.dotype( jtype ) )
580 ioldsd( j ) = iseed( j )
605 IF( mtypes.GT.maxtyp )
608 itype = ktype( jtype )
609 imode = kmode( jtype )
613 GO TO ( 40, 50, 60 )kmagn( jtype )
620 anorm = ( rtovfl*ulp )*amninv
624 anorm = rtunfl*max( m, n )*ulpinv
629 CALL claset(
'Full', lda, n, czero, czero, a, lda )
634 IF( itype.EQ.1 )
THEN
640 ELSE IF( itype.EQ.2 )
THEN
644 DO 80 jcol = 1, mnmin
645 a( jcol, jcol ) = anorm
648 ELSE IF( itype.EQ.4 )
THEN
652 CALL clatms( mnmin, mnmin,
'S', iseed,
'N', rwork, imode,
653 $ cond, anorm, 0, 0,
'N', a, lda, work,
656 ELSE IF( itype.EQ.5 )
THEN
660 CALL clatms( mnmin, mnmin,
'S', iseed,
'S', rwork, imode,
661 $ cond, anorm, m, n,
'N', a, lda, work,
664 ELSE IF( itype.EQ.6 )
THEN
668 CALL clatms( m, n,
'S', iseed,
'N', rwork, imode, cond,
669 $ anorm, m, n,
'N', a, lda, work, iinfo )
671 ELSE IF( itype.EQ.7 )
THEN
675 CALL clatmr( mnmin, mnmin,
'S', iseed,
'N', work, 6, one,
676 $ cone,
'T',
'N', work( mnmin+1 ), 1, one,
677 $ work( 2*mnmin+1 ), 1, one,
'N', iwork, 0, 0,
678 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
680 ELSE IF( itype.EQ.8 )
THEN
684 CALL clatmr( mnmin, mnmin,
'S', iseed,
'S', work, 6, one,
685 $ cone,
'T',
'N', work( mnmin+1 ), 1, one,
686 $ work( m+mnmin+1 ), 1, one,
'N', iwork, m, n,
687 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
689 ELSE IF( itype.EQ.9 )
THEN
693 CALL clatmr( m, n,
'S', iseed,
'N', work, 6, one, cone,
694 $
'T',
'N', work( mnmin+1 ), 1, one,
695 $ work( m+mnmin+1 ), 1, one,
'N', iwork, m, n,
696 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
698 ELSE IF( itype.EQ.10 )
THEN
702 temp1 = -two*log( ulp )
704 bd( j ) = exp( temp1*slarnd( 2, iseed ) )
706 $ be( j ) = exp( temp1*slarnd( 2, iseed ) )
720 IF( iinfo.EQ.0 )
THEN
725 CALL clatmr( mnmin, nrhs,
'S', iseed,
'N', work, 6,
726 $ one, cone,
'T',
'N', work( mnmin+1 ), 1,
727 $ one, work( 2*mnmin+1 ), 1, one,
'N',
728 $ iwork, mnmin, nrhs, zero, one,
'NO', y,
729 $ ldx, iwork, iinfo )
731 CALL clatmr( m, nrhs,
'S', iseed,
'N', work, 6, one,
732 $ cone,
'T',
'N', work( m+1 ), 1, one,
733 $ work( 2*m+1 ), 1, one,
'N', iwork, m,
734 $ nrhs, zero, one,
'NO', x, ldx, iwork,
741 IF( iinfo.NE.0 )
THEN
742 WRITE( nout, fmt = 9998 )
'Generator', iinfo, m, n,
752 IF( .NOT.bidiag )
THEN
757 CALL clacpy(
' ', m, n, a, lda, q, ldq )
758 CALL cgebrd( m, n, q, ldq, bd, be, work, work( mnmin+1 ),
759 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
763 IF( iinfo.NE.0 )
THEN
764 WRITE( nout, fmt = 9998 )
'CGEBRD', iinfo, m, n,
770 CALL clacpy(
' ', m, n, q, ldq, pt, ldpt )
782 CALL cungbr(
'Q', m, mq, n, q, ldq, work,
783 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
787 IF( iinfo.NE.0 )
THEN
788 WRITE( nout, fmt = 9998 )
'CUNGBR(Q)', iinfo, m, n,
796 CALL cungbr(
'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
797 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
801 IF( iinfo.NE.0 )
THEN
802 WRITE( nout, fmt = 9998 )
'CUNGBR(P)', iinfo, m, n,
810 CALL cgemm(
'Conjugate transpose',
'No transpose', m,
811 $ nrhs, m, cone, q, ldq, x, ldx, czero, y,
818 CALL cbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
819 $ work, rwork, result( 1 ) )
820 CALL cunt01(
'Columns', m, mq, q, ldq, work, lwork,
821 $ rwork, result( 2 ) )
822 CALL cunt01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
823 $ rwork, result( 3 ) )
829 CALL scopy( mnmin, bd, 1, s1, 1 )
831 $
CALL scopy( mnmin-1, be, 1, rwork, 1 )
832 CALL clacpy(
' ', m, nrhs, y, ldx, z, ldx )
833 CALL claset(
'Full', mnmin, mnmin, czero, cone, u, ldpt )
834 CALL claset(
'Full', mnmin, mnmin, czero, cone, vt, ldpt )
836 CALL cbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, rwork, vt,
837 $ ldpt, u, ldpt, z, ldx, rwork( mnmin+1 ),
842 IF( iinfo.NE.0 )
THEN
843 WRITE( nout, fmt = 9998 )
'CBDSQR(vects)', iinfo, m, n,
846 IF( iinfo.LT.0 )
THEN
857 CALL scopy( mnmin, bd, 1, s2, 1 )
859 $
CALL scopy( mnmin-1, be, 1, rwork, 1 )
861 CALL cbdsqr( uplo, mnmin, 0, 0, 0, s2, rwork, vt, ldpt, u,
862 $ ldpt, z, ldx, rwork( mnmin+1 ), iinfo )
866 IF( iinfo.NE.0 )
THEN
867 WRITE( nout, fmt = 9998 )
'CBDSQR(values)', iinfo, m, n,
870 IF( iinfo.LT.0 )
THEN
883 CALL cbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
884 $ work, result( 4 ) )
885 CALL cbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
886 $ rwork, result( 5 ) )
887 CALL cunt01(
'Columns', mnmin, mnmin, u, ldpt, work, lwork,
888 $ rwork, result( 6 ) )
889 CALL cunt01(
'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
890 $ rwork, result( 7 ) )
896 DO 110 i = 1, mnmin - 1
897 IF( s1( i ).LT.s1( i+1 ) )
898 $ result( 8 ) = ulpinv
899 IF( s1( i ).LT.zero )
900 $ result( 8 ) = ulpinv
902 IF( mnmin.GE.1 )
THEN
903 IF( s1( mnmin ).LT.zero )
904 $ result( 8 ) = ulpinv
912 temp1 = abs( s1( j )-s2( j ) ) /
913 $ max( sqrt( unfl )*max( s1( 1 ), one ),
914 $ ulp*max( abs( s1( j ) ), abs( s2( j ) ) ) )
915 temp2 = max( temp1, temp2 )
923 temp1 = thresh*( half-ulp )
926 CALL ssvdch( mnmin, bd, be, s1, temp1, iinfo )
938 IF( .NOT.bidiag )
THEN
939 CALL scopy( mnmin, bd, 1, s2, 1 )
941 $
CALL scopy( mnmin-1, be, 1, rwork, 1 )
943 CALL cbdsqr( uplo, mnmin, n, m, nrhs, s2, rwork, pt,
944 $ ldpt, q, ldq, y, ldx, rwork( mnmin+1 ),
952 CALL cbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
953 $ ldpt, work, rwork, result( 11 ) )
954 CALL cbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
955 $ rwork, result( 12 ) )
956 CALL cunt01(
'Columns', m, mq, q, ldq, work, lwork,
957 $ rwork, result( 13 ) )
958 CALL cunt01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
959 $ rwork, result( 14 ) )
966 IF( result( j ).GE.thresh )
THEN
968 $
CALL slahd2( nout, path )
969 WRITE( nout, fmt = 9999 )m, n, jtype, ioldsd, j,
974 IF( .NOT.bidiag )
THEN
985 CALL alasum( path, nout, nfail, ntest, 0 )
991 9999
FORMAT(
' M=', i5,
', N=', i5,
', type ', i2,
', seed=',
992 $ 4( i4,
',' ),
' test(', i2,
')=', g11.4 )
993 9998
FORMAT(
' CCHKBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
994 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
subroutine cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine cunt01(ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, RESID)
CUNT01
subroutine clatmr(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)
CLATMR
subroutine cbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RWORK, RESID)
CBDT01
subroutine clatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
CLATMS
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
subroutine slahd2(IOUNIT, PATH)
SLAHD2
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cbdt02(M, N, B, LDB, C, LDC, U, LDU, WORK, RWORK, RESID)
CBDT02
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine cchkbd(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, RWORK, NOUT, INFO)
CCHKBD
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine cbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
CBDT03
subroutine ssvdch(N, S, E, SVD, TOL, INFO)
SSVDCH
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR