413 SUBROUTINE zchkbd( 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,
426 DOUBLE PRECISION thresh
430 INTEGER iseed( 4 ), mval( * ), nval( * )
431 DOUBLE PRECISION bd( * ), be( * ), rwork( * ), s1( * ), s2( * )
432 COMPLEX*16 a( lda, * ), pt( ldpt, * ), q( ldq, * ),
433 $ u( ldpt, * ), vt( ldpt, * ), work( * ),
434 $ x( ldx, * ), y( ldx, * ), z( ldx, * )
440 DOUBLE PRECISION zero, one, two, half
441 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
443 COMPLEX*16 czero, cone
444 parameter( czero = ( 0.0d+0, 0.0d+0 ),
445 $ cone = ( 1.0d+0, 0.0d+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 DOUBLE PRECISION 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 DOUBLE PRECISION 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(
'ZCHKBD', -info )
547 path( 1: 1 ) =
'Zomplex precision'
551 unfl =
dlamch(
'Safe minimum' )
552 ovfl =
dlamch(
'Overflow' )
554 ulp =
dlamch(
'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
zlaset(
'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
zlatms( 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
zlatms( 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
zlatms( 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
zlatmr( 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
zlatmr( 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
zlatmr( 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*
dlarnd( 2, iseed ) )
706 $ be(
j ) = exp( temp1*
dlarnd( 2, iseed ) )
720 IF( iinfo.EQ.0 )
THEN
725 CALL
zlatmr( 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
zlatmr( 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
zlacpy(
' ', m, n, a, lda, q, ldq )
758 CALL
zgebrd( 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 )
'ZGEBRD', iinfo, m, n,
770 CALL
zlacpy(
' ', m, n, q, ldq, pt, ldpt )
782 CALL
zungbr(
'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 )
'ZUNGBR(Q)', iinfo, m, n,
796 CALL
zungbr(
'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 )
'ZUNGBR(P)', iinfo, m, n,
810 CALL
zgemm(
'Conjugate transpose',
'No transpose', m,
811 $ nrhs, m, cone, q, ldq, x, ldx, czero, y,
818 CALL
zbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
819 $ work, rwork, result( 1 ) )
820 CALL
zunt01(
'Columns', m, mq, q, ldq, work, lwork,
821 $ rwork, result( 2 ) )
822 CALL
zunt01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
823 $ rwork, result( 3 ) )
829 CALL
dcopy( mnmin, bd, 1, s1, 1 )
831 $ CALL
dcopy( mnmin-1, be, 1, rwork, 1 )
832 CALL
zlacpy(
' ', m, nrhs, y, ldx, z, ldx )
833 CALL
zlaset(
'Full', mnmin, mnmin, czero, cone, u, ldpt )
834 CALL
zlaset(
'Full', mnmin, mnmin, czero, cone, vt, ldpt )
836 CALL
zbdsqr( 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 )
'ZBDSQR(vects)', iinfo, m, n,
846 IF( iinfo.LT.0 )
THEN
857 CALL
dcopy( mnmin, bd, 1, s2, 1 )
859 $ CALL
dcopy( mnmin-1, be, 1, rwork, 1 )
861 CALL
zbdsqr( 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 )
'ZBDSQR(values)', iinfo, m, n,
870 IF( iinfo.LT.0 )
THEN
883 CALL
zbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
884 $ work, result( 4 ) )
885 CALL
zbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
886 $ rwork, result( 5 ) )
887 CALL
zunt01(
'Columns', mnmin, mnmin, u, ldpt, work, lwork,
888 $ rwork, result( 6 ) )
889 CALL
zunt01(
'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
dsvdch( mnmin, bd, be, s1, temp1, iinfo )
938 IF( .NOT.bidiag )
THEN
939 CALL
dcopy( mnmin, bd, 1, s2, 1 )
941 $ CALL
dcopy( mnmin-1, be, 1, rwork, 1 )
943 CALL
zbdsqr( uplo, mnmin, n, m, nrhs, s2, rwork, pt,
944 $ ldpt, q, ldq, y, ldx, rwork( mnmin+1 ),
952 CALL
zbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
953 $ ldpt, work, rwork, result( 11 ) )
954 CALL
zbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
955 $ rwork, result( 12 ) )
956 CALL
zunt01(
'Columns', m, mq, q, ldq, work, lwork,
957 $ rwork, result( 13 ) )
958 CALL
zunt01(
'Rows', mnmin, n, pt, ldpt, work, lwork,
959 $ rwork, result( 14 ) )
966 IF( result(
j ).GE.thresh )
THEN
968 $ CALL
dlahd2( 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(
' ZCHKBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
994 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
double precision function dlarnd(IDIST, ISEED)
DLARND
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zunt01(ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, RESID)
ZUNT01
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RWORK, RESID)
ZBDT01
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zbdt02(M, N, B, LDB, C, LDC, U, LDU, WORK, RWORK, RESID)
ZBDT02
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
subroutine zlatmr(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)
ZLATMR
subroutine dlahd2(IOUNIT, PATH)
DLAHD2
subroutine zbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
ZBDT03
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine zlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
ZLATMS
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
subroutine zchkbd(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)
ZCHKBD
subroutine dsvdch(N, S, E, SVD, TOL, INFO)
DSVDCH