329 SUBROUTINE zdrvbd( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
330 $ a, lda, u, ldu, vt, ldvt, asav, usav, vtsav, s,
331 $ ssav, e, work, lwork, rwork, iwork, nounit,
340 INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES,
342 DOUBLE PRECISION THRESH
346 INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
347 DOUBLE PRECISION E( * ), RWORK( * ), S( * ), SSAV( * )
348 COMPLEX*16 A( lda, * ), ASAV( lda, * ), U( ldu, * ),
349 $ usav( ldu, * ), vt( ldvt, * ),
350 $ vtsav( ldvt, * ), work( * )
356 DOUBLE PRECISION ZERO, ONE
357 parameter( zero = 0.0d+0, one = 1.0d+0 )
358 COMPLEX*16 CZERO, CONE
359 parameter( czero = ( 0.0d+0, 0.0d+0 ),
360 $ cone = ( 1.0d+0, 0.0d+0 ) )
362 parameter( maxtyp = 5 )
366 CHARACTER JOBQ, JOBU, JOBVT
367 INTEGER I, IINFO, IJQ, IJU, IJVT, IWSPC, IWTMP, J,
368 $ jsize, jtype, lswork, m, minwrk, mmax, mnmax,
369 $ mnmin, mtypes, n, nerrs, nfail, nmax, ntest,
371 DOUBLE PRECISION ANORM, DIF, DIV, OVFL, ULP, ULPINV, UNFL
376 DOUBLE PRECISION RESULT( 14 )
379 DOUBLE PRECISION DLAMCH
387 INTRINSIC abs, dble, max, min
390 DATA cjob /
'N',
'O',
'S',
'A' /
410 mmax = max( mmax, mm( j ) )
413 nmax = max( nmax, nn( j ) )
416 mnmax = max( mnmax, min( mm( j ), nn( j ) ) )
417 minwrk = max( minwrk, max( 3*min( mm( j ),
418 $ nn( j ) )+max( mm( j ), nn( j ) )**2, 5*min( mm( j ),
419 $ nn( j ) ), 3*max( mm( j ), nn( j ) ) ) )
424 IF( nsizes.LT.0 )
THEN
426 ELSE IF( badmm )
THEN
428 ELSE IF( badnn )
THEN
430 ELSE IF( ntypes.LT.0 )
THEN
432 ELSE IF( lda.LT.max( 1, mmax ) )
THEN
434 ELSE IF( ldu.LT.max( 1, mmax ) )
THEN
436 ELSE IF( ldvt.LT.max( 1, nmax ) )
THEN
438 ELSE IF( minwrk.GT.lwork )
THEN
443 CALL
xerbla(
'ZDRVBD', -info )
449 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
463 DO 180 jsize = 1, nsizes
468 IF( nsizes.NE.1 )
THEN
469 mtypes = min( maxtyp, ntypes )
471 mtypes = min( maxtyp+1, ntypes )
474 DO 170 jtype = 1, mtypes
475 IF( .NOT.dotype( jtype ) )
480 ioldsd( j ) = iseed( j )
485 IF( mtypes.GT.maxtyp )
488 IF( jtype.EQ.1 )
THEN
492 CALL
zlaset(
'Full', m, n, czero, czero, a, lda )
493 DO 30 i = 1, min( m, n )
497 ELSE IF( jtype.EQ.2 )
THEN
501 CALL
zlaset(
'Full', m, n, czero, cone, a, lda )
502 DO 40 i = 1, min( m, n )
516 CALL
zlatms( m, n,
'U', iseed,
'N', s, 4, dble( mnmin ),
517 $ anorm, m-1, n-1,
'N', a, lda, work, iinfo )
518 IF( iinfo.NE.0 )
THEN
519 WRITE( nounit, fmt = 9996 )
'Generator', iinfo, m, n,
527 CALL
zlacpy(
'F', m, n, a, lda, asav, lda )
535 iwtmp = 2*min( m, n )+max( m, n )
536 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
537 lswork = min( lswork, lwork )
538 lswork = max( lswork, 1 )
549 $ CALL
zlacpy(
'F', m, n, asav, lda, a, lda )
550 CALL
zgesvd(
'A',
'A', m, n, a, lda, ssav, usav, ldu,
551 $ vtsav, ldvt, work, lswork, rwork, iinfo )
552 IF( iinfo.NE.0 )
THEN
553 WRITE( nounit, fmt = 9995 )
'GESVD', iinfo, m, n,
554 $ jtype, lswork, ioldsd
561 CALL
zbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
562 $ vtsav, ldvt, work, rwork, result( 1 ) )
563 IF( m.NE.0 .AND. n.NE.0 )
THEN
564 CALL
zunt01(
'Columns', mnmin, m, usav, ldu, work,
565 $ lwork, rwork, result( 2 ) )
566 CALL
zunt01(
'Rows', mnmin, n, vtsav, ldvt, work,
567 $ lwork, rwork, result( 3 ) )
570 DO 70 i = 1, mnmin - 1
571 IF( ssav( i ).LT.ssav( i+1 ) )
572 $ result( 4 ) = ulpinv
573 IF( ssav( i ).LT.zero )
574 $ result( 4 ) = ulpinv
576 IF( mnmin.GE.1 )
THEN
577 IF( ssav( mnmin ).LT.zero )
578 $ result( 4 ) = ulpinv
588 IF( ( iju.EQ.3 .AND. ijvt.EQ.3 ) .OR.
589 $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )go to 90
591 jobvt = cjob( ijvt+1 )
592 CALL
zlacpy(
'F', m, n, asav, lda, a, lda )
593 CALL
zgesvd( jobu, jobvt, m, n, a, lda, s, u, ldu,
594 $ vt, ldvt, work, lswork, rwork, iinfo )
599 IF( m.GT.0 .AND. n.GT.0 )
THEN
601 CALL
zunt03(
'C', m, mnmin, m, mnmin, usav,
602 $ ldu, a, lda, work, lwork, rwork,
604 ELSE IF( iju.EQ.2 )
THEN
605 CALL
zunt03(
'C', m, mnmin, m, mnmin, usav,
606 $ ldu, u, ldu, work, lwork, rwork,
608 ELSE IF( iju.EQ.3 )
THEN
609 CALL
zunt03(
'C', m, m, m, mnmin, usav, ldu,
610 $ u, ldu, work, lwork, rwork, dif,
614 result( 5 ) = max( result( 5 ), dif )
619 IF( m.GT.0 .AND. n.GT.0 )
THEN
621 CALL
zunt03(
'R', n, mnmin, n, mnmin, vtsav,
622 $ ldvt, a, lda, work, lwork,
623 $ rwork, dif, iinfo )
624 ELSE IF( ijvt.EQ.2 )
THEN
625 CALL
zunt03(
'R', n, mnmin, n, mnmin, vtsav,
626 $ ldvt, vt, ldvt, work, lwork,
627 $ rwork, dif, iinfo )
628 ELSE IF( ijvt.EQ.3 )
THEN
629 CALL
zunt03(
'R', n, n, n, mnmin, vtsav,
630 $ ldvt, vt, ldvt, work, lwork,
631 $ rwork, dif, iinfo )
634 result( 6 ) = max( result( 6 ), dif )
639 div = max( dble( mnmin )*ulp*s( 1 ),
640 $ dlamch(
'Safe minimum' ) )
641 DO 80 i = 1, mnmin - 1
642 IF( ssav( i ).LT.ssav( i+1 ) )
644 IF( ssav( i ).LT.zero )
646 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
648 result( 7 ) = max( result( 7 ), dif )
654 iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
655 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
656 lswork = min( lswork, lwork )
657 lswork = max( lswork, 1 )
663 CALL
zlacpy(
'F', m, n, asav, lda, a, lda )
664 CALL
zgesdd(
'A', m, n, a, lda, ssav, usav, ldu, vtsav,
665 $ ldvt, work, lswork, rwork, iwork, iinfo )
666 IF( iinfo.NE.0 )
THEN
667 WRITE( nounit, fmt = 9995 )
'GESDD', iinfo, m, n,
668 $ jtype, lswork, ioldsd
675 CALL
zbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
676 $ vtsav, ldvt, work, rwork, result( 8 ) )
677 IF( m.NE.0 .AND. n.NE.0 )
THEN
678 CALL
zunt01(
'Columns', mnmin, m, usav, ldu, work,
679 $ lwork, rwork, result( 9 ) )
680 CALL
zunt01(
'Rows', mnmin, n, vtsav, ldvt, work,
681 $ lwork, rwork, result( 10 ) )
684 DO 110 i = 1, mnmin - 1
685 IF( ssav( i ).LT.ssav( i+1 ) )
686 $ result( 11 ) = ulpinv
687 IF( ssav( i ).LT.zero )
688 $ result( 11 ) = ulpinv
690 IF( mnmin.GE.1 )
THEN
691 IF( ssav( mnmin ).LT.zero )
692 $ result( 11 ) = ulpinv
702 CALL
zlacpy(
'F', m, n, asav, lda, a, lda )
703 CALL
zgesdd( jobq, m, n, a, lda, s, u, ldu, vt, ldvt,
704 $ work, lswork, rwork, iwork, iinfo )
709 IF( m.GT.0 .AND. n.GT.0 )
THEN
712 CALL
zunt03(
'C', m, mnmin, m, mnmin, usav,
713 $ ldu, a, lda, work, lwork, rwork,
716 CALL
zunt03(
'C', m, mnmin, m, mnmin, usav,
717 $ ldu, u, ldu, work, lwork, rwork,
720 ELSE IF( ijq.EQ.2 )
THEN
721 CALL
zunt03(
'C', m, mnmin, m, mnmin, usav, ldu,
722 $ u, ldu, work, lwork, rwork, dif,
726 result( 12 ) = max( result( 12 ), dif )
731 IF( m.GT.0 .AND. n.GT.0 )
THEN
734 CALL
zunt03(
'R', n, mnmin, n, mnmin, vtsav,
735 $ ldvt, vt, ldvt, work, lwork,
736 $ rwork, dif, iinfo )
738 CALL
zunt03(
'R', n, mnmin, n, mnmin, vtsav,
739 $ ldvt, a, lda, work, lwork,
740 $ rwork, dif, iinfo )
742 ELSE IF( ijq.EQ.2 )
THEN
743 CALL
zunt03(
'R', n, mnmin, n, mnmin, vtsav,
744 $ ldvt, vt, ldvt, work, lwork, rwork,
748 result( 13 ) = max( result( 13 ), dif )
753 div = max( dble( mnmin )*ulp*s( 1 ),
754 $ dlamch(
'Safe minimum' ) )
755 DO 120 i = 1, mnmin - 1
756 IF( ssav( i ).LT.ssav( i+1 ) )
758 IF( ssav( i ).LT.zero )
760 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
762 result( 14 ) = max( result( 14 ), dif )
770 IF( result( j ).GE.zero )
772 IF( result( j ).GE.thresh )
777 $ ntestf = ntestf + 1
778 IF( ntestf.EQ.1 )
THEN
779 WRITE( nounit, fmt = 9999 )
780 WRITE( nounit, fmt = 9998 )thresh
785 IF( result( j ).GE.thresh )
THEN
786 WRITE( nounit, fmt = 9997 )m, n, jtype, iwspc,
787 $ ioldsd, j, result( j )
791 nerrs = nerrs + nfail
792 ntestt = ntestt + ntest
801 CALL
alasvm(
'ZBD', nounit, nerrs, ntestt, 0 )
803 9999
FORMAT(
' SVD -- Complex Singular Value Decomposition Driver ',
804 $ /
' Matrix types (see ZDRVBD for details):',
805 $ / /
' 1 = Zero matrix', /
' 2 = Identity matrix',
806 $ /
' 3 = Evenly spaced singular values near 1',
807 $ /
' 4 = Evenly spaced singular values near underflow',
808 $ /
' 5 = Evenly spaced singular values near overflow',
809 $ / /
' Tests performed: ( A is dense, U and V are unitary,',
810 $ / 19x,
' S is an array, and Upartial, VTpartial, and',
811 $ / 19x,
' Spartial are partially computed U, VT and S),', / )
812 9998
FORMAT(
' Tests performed with Test Threshold = ', f8.2,
814 $
' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
815 $ /
' 2 = | I - U**T U | / ( M ulp ) ',
816 $ /
' 3 = | I - VT VT**T | / ( N ulp ) ',
817 $ /
' 4 = 0 if S contains min(M,N) nonnegative values in',
818 $
' decreasing order, else 1/ulp',
819 $ /
' 5 = | U - Upartial | / ( M ulp )',
820 $ /
' 6 = | VT - VTpartial | / ( N ulp )',
821 $ /
' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
823 $
' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
824 $ /
' 9 = | I - U**T U | / ( M ulp ) ',
825 $ /
'10 = | I - VT VT**T | / ( N ulp ) ',
826 $ /
'11 = 0 if S contains min(M,N) nonnegative values in',
827 $
' decreasing order, else 1/ulp',
828 $ /
'12 = | U - Upartial | / ( M ulp )',
829 $ /
'13 = | VT - VTpartial | / ( N ulp )',
830 $ /
'14 = | S - Spartial | / ( min(M,N) ulp |S| )', / / )
831 9997
FORMAT(
' M=', i5,
', N=', i5,
', type ', i1,
', IWS=', i1,
832 $
', seed=', 4( i4,
',' ),
' test(', i1,
')=', g11.4 )
833 9996
FORMAT(
' ZDRVBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
834 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
836 9995
FORMAT(
' ZDRVBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
837 $ i6,
', N=', i6,
', JTYPE=', i6,
', LSWORK=', i6, / 9x,
838 $
'ISEED=(', 3( i5,
',' ), i5,
')' )
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine zdrvbd(NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH, A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S, SSAV, E, WORK, LWORK, RWORK, IWORK, NOUNIT, INFO)
ZDRVBD
subroutine zgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO)
ZGESDD
subroutine zunt01(ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, RESID)
ZUNT01
subroutine zunt03(RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK, RWORK, RESULT, INFO)
ZUNT03
subroutine zbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RWORK, RESID)
ZBDT01
subroutine xerbla(SRNAME, INFO)
XERBLA
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 zgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
ZGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine zlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
ZLATMS