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
475 INTRINSIC abs, exp, int, log, max, min, sqrt
483 COMMON / infoc / infot, nunit, ok, lerr
484 COMMON / srnamc / srnamt
487 DATA ktype / 1, 2, 5*4, 5*6, 3*9, 10 /
488 DATA kmagn / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
489 DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
505 mmax = max( mmax, mval( j ) )
508 nmax = max( nmax, nval( j ) )
511 mnmax = max( mnmax, min( mval( j ), nval( j ) ) )
512 minwrk = max( minwrk, 3*( mval( j )+nval( j ) ),
513 $ mval( j )*( mval( j )+max( mval( j ), nval( j ),
514 $ nrhs )+1 )+nval( j )*min( nval( j ), mval( j ) ) )
519 IF( nsizes.LT.0 )
THEN
521 ELSE IF( badmm )
THEN
523 ELSE IF( badnn )
THEN
525 ELSE IF( ntypes.LT.0 )
THEN
527 ELSE IF( nrhs.LT.0 )
THEN
529 ELSE IF( lda.LT.mmax )
THEN
531 ELSE IF( ldx.LT.mmax )
THEN
533 ELSE IF( ldq.LT.mmax )
THEN
535 ELSE IF( ldpt.LT.mnmax )
THEN
537 ELSE IF( minwrk.GT.lwork )
THEN
542 CALL xerbla(
'CCHKBD', -info )
548 path( 1: 1 ) =
'Complex precision'
552 unfl = slamch(
'Safe minimum' )
553 ovfl = slamch(
'Overflow' )
555 ulp = slamch(
'Precision' )
557 log2ui = int( log( ulpinv ) / log( two ) )
558 rtunfl = sqrt( unfl )
559 rtovfl = sqrt( ovfl )
564 DO 180 jsize = 1, nsizes
568 amninv = one / max( m, n, 1 )
570 IF( nsizes.NE.1 )
THEN
571 mtypes = min( maxtyp, ntypes )
573 mtypes = min( maxtyp+1, ntypes )
576 DO 170 jtype = 1, mtypes
577 IF( .NOT.dotype( jtype ) )
581 ioldsd( j ) = iseed( j )
606 IF( mtypes.GT.maxtyp )
609 itype = ktype( jtype )
610 imode = kmode( jtype )
614 GO TO ( 40, 50, 60 )kmagn( jtype )
621 anorm = ( rtovfl*ulp )*amninv
625 anorm = rtunfl*max( m, n )*ulpinv
630 CALL claset(
'Full', lda, n, czero, czero, a, lda )
635 IF( itype.EQ.1 )
THEN
641 ELSE IF( itype.EQ.2 )
THEN
645 DO 80 jcol = 1, mnmin
646 a( jcol, jcol ) = anorm
649 ELSE IF( itype.EQ.4 )
THEN
653 CALL clatms( mnmin, mnmin,
'S', iseed,
'N', rwork, imode,
654 $ cond, anorm, 0, 0,
'N', a, lda, work,
657 ELSE IF( itype.EQ.5 )
THEN
661 CALL clatms( mnmin, mnmin,
'S', iseed,
'S', rwork, imode,
662 $ cond, anorm, m, n,
'N', a, lda, work,
665 ELSE IF( itype.EQ.6 )
THEN
669 CALL clatms( m, n,
'S', iseed,
'N', rwork, imode, cond,
670 $ anorm, m, n,
'N', a, lda, work, iinfo )
672 ELSE IF( itype.EQ.7 )
THEN
676 CALL clatmr( mnmin, mnmin,
'S', iseed,
'N', work, 6, one,
677 $ cone,
'T',
'N', work( mnmin+1 ), 1, one,
678 $ work( 2*mnmin+1 ), 1, one,
'N', iwork, 0, 0,
679 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
681 ELSE IF( itype.EQ.8 )
THEN
685 CALL clatmr( mnmin, mnmin,
'S', iseed,
'S', work, 6, one,
686 $ cone,
'T',
'N', work( mnmin+1 ), 1, one,
687 $ work( m+mnmin+1 ), 1, one,
'N', iwork, m, n,
688 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
690 ELSE IF( itype.EQ.9 )
THEN
694 CALL clatmr( m, n,
'S', iseed,
'N', work, 6, one, cone,
695 $
'T',
'N', work( mnmin+1 ), 1, one,
696 $ work( m+mnmin+1 ), 1, one,
'N', iwork, m, n,
697 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
699 ELSE IF( itype.EQ.10 )
THEN
703 temp1 = -two*log( ulp )
705 bd( j ) = exp( temp1*slarnd( 2, iseed ) )
707 $ be( j ) = exp( temp1*slarnd( 2, iseed ) )
721 IF( iinfo.EQ.0 )
THEN
726 CALL clatmr( mnmin, nrhs,
'S', iseed,
'N', work, 6,
727 $ one, cone,
'T',
'N', work( mnmin+1 ), 1,
728 $ one, work( 2*mnmin+1 ), 1, one,
'N',
729 $ iwork, mnmin, nrhs, zero, one,
'NO', y,
730 $ ldx, iwork, iinfo )
732 CALL clatmr( m, nrhs,
'S', iseed,
'N', work, 6, one,
733 $ cone,
'T',
'N', work( m+1 ), 1, one,
734 $ work( 2*m+1 ), 1, one,
'N', iwork, m,
735 $ nrhs, zero, one,
'NO', x, ldx, iwork,
742 IF( iinfo.NE.0 )
THEN
743 WRITE( nout, fmt = 9998 )
'Generator', iinfo, m, n,
753 IF( .NOT.bidiag )
THEN
758 CALL clacpy(
' ', m, n, a, lda, q, ldq )
759 CALL cgebrd( m, n, q, ldq, bd, be, work, work( mnmin+1 ),
760 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
764 IF( iinfo.NE.0 )
THEN
765 WRITE( nout, fmt = 9998 )
'CGEBRD', iinfo, m, n,
771 CALL clacpy(
' ', m, n, q, ldq, pt, ldpt )
783 CALL cungbr(
'Q', m, mq, n, q, ldq, work,
784 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
788 IF( iinfo.NE.0 )
THEN
789 WRITE( nout, fmt = 9998 )
'CUNGBR(Q)', iinfo, m, n,
797 CALL cungbr(
'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
798 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
802 IF( iinfo.NE.0 )
THEN
803 WRITE( nout, fmt = 9998 )
'CUNGBR(P)', iinfo, m, n,
811 CALL cgemm(
'Conjugate transpose',
'No transpose', m,
812 $ nrhs, m, cone, q, ldq, x, ldx, czero, y,
819 CALL cbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
820 $ work, rwork, result( 1 ) )
821 CALL cunt01(
'Columns', m, mq, q, ldq, work, lwork,
822 $ rwork, result( 2 ) )
823 CALL cunt01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
824 $ rwork, result( 3 ) )
830 CALL scopy( mnmin, bd, 1, s1, 1 )
832 $
CALL scopy( mnmin-1, be, 1, rwork, 1 )
833 CALL clacpy(
' ', m, nrhs, y, ldx, z, ldx )
834 CALL claset(
'Full', mnmin, mnmin, czero, cone, u, ldpt )
835 CALL claset(
'Full', mnmin, mnmin, czero, cone, vt, ldpt )
837 CALL cbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, rwork, vt,
838 $ ldpt, u, ldpt, z, ldx, rwork( mnmin+1 ),
843 IF( iinfo.NE.0 )
THEN
844 WRITE( nout, fmt = 9998 )
'CBDSQR(vects)', iinfo, m, n,
847 IF( iinfo.LT.0 )
THEN
858 CALL scopy( mnmin, bd, 1, s2, 1 )
860 $
CALL scopy( mnmin-1, be, 1, rwork, 1 )
862 CALL cbdsqr( uplo, mnmin, 0, 0, 0, s2, rwork, vt, ldpt, u,
863 $ ldpt, z, ldx, rwork( mnmin+1 ), iinfo )
867 IF( iinfo.NE.0 )
THEN
868 WRITE( nout, fmt = 9998 )
'CBDSQR(values)', iinfo, m, n,
871 IF( iinfo.LT.0 )
THEN
884 CALL cbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
885 $ work, result( 4 ) )
886 CALL cbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
887 $ rwork, result( 5 ) )
888 CALL cunt01(
'Columns', mnmin, mnmin, u, ldpt, work, lwork,
889 $ rwork, result( 6 ) )
890 CALL cunt01(
'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
891 $ rwork, result( 7 ) )
897 DO 110 i = 1, mnmin - 1
898 IF( s1( i ).LT.s1( i+1 ) )
899 $ result( 8 ) = ulpinv
900 IF( s1( i ).LT.zero )
901 $ result( 8 ) = ulpinv
903 IF( mnmin.GE.1 )
THEN
904 IF( s1( mnmin ).LT.zero )
905 $ result( 8 ) = ulpinv
913 temp1 = abs( s1( j )-s2( j ) ) /
914 $ max( sqrt( unfl )*max( s1( 1 ), one ),
915 $ ulp*max( abs( s1( j ) ), abs( s2( j ) ) ) )
916 temp2 = max( temp1, temp2 )
924 temp1 = thresh*( half-ulp )
927 CALL ssvdch( mnmin, bd, be, s1, temp1, iinfo )
939 IF( .NOT.bidiag )
THEN
940 CALL scopy( mnmin, bd, 1, s2, 1 )
942 $
CALL scopy( mnmin-1, be, 1, rwork, 1 )
944 CALL cbdsqr( uplo, mnmin, n, m, nrhs, s2, rwork, pt,
945 $ ldpt, q, ldq, y, ldx, rwork( mnmin+1 ),
953 CALL cbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
954 $ ldpt, work, rwork, result( 11 ) )
955 CALL cbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
956 $ rwork, result( 12 ) )
957 CALL cunt01(
'Columns', m, mq, q, ldq, work, lwork,
958 $ rwork, result( 13 ) )
959 CALL cunt01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
960 $ rwork, result( 14 ) )
967 IF( result( j ).GE.thresh )
THEN
969 $
CALL slahd2( nout, path )
970 WRITE( nout, fmt = 9999 )m, n, jtype, ioldsd, j,
975 IF( .NOT.bidiag )
THEN
986 CALL alasum( path, nout, nfail, ntest, 0 )
992 9999
FORMAT(
' M=', i5,
', N=', i5,
', type ', i2,
', seed=',
993 $ 4( i4,
',' ),
' test(', i2,
')=', g11.4 )
994 9998
FORMAT(
' CCHKBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
995 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
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 cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slahd2(IOUNIT, PATH)
SLAHD2
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RWORK, RESID)
CBDT01
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 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 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 clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine ssvdch(N, S, E, SVD, TOL, INFO)
SSVDCH
subroutine cbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
CBDT03
subroutine cbdt02(M, N, B, LDB, C, LDC, U, LDU, WORK, RWORK, RESID)
CBDT02
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
subroutine cunt01(ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, RESID)
CUNT01