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