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,
',' ),
subroutine alasum(type, nout, nfail, nrun, nerrs)
ALASUM
subroutine cbdt01(m, n, kd, a, lda, q, ldq, d, e, pt, ldpt, work, rwork, resid)
CBDT01
subroutine cbdt02(m, n, b, ldb, c, ldc, u, ldu, work, rwork, resid)
CBDT02
subroutine cbdt03(uplo, n, kd, d, e, u, ldu, s, vt, ldvt, work, resid)
CBDT03
subroutine xerbla(srname, info)
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 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 clatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
CLATMS
subroutine cunt01(rowcol, m, n, u, ldu, work, lwork, rwork, resid)
CUNT01
subroutine cbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
CBDSQR
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
subroutine cgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
CGEBRD
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
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 cungbr(vect, m, n, k, a, lda, tau, work, lwork, info)
CUNGBR
subroutine slahd2(iounit, path)
SLAHD2
subroutine ssvdch(n, s, e, svd, tol, info)
SSVDCH