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 )
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,
',' ),