329 SUBROUTINE cdrvbd( 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,
346 INTEGER iseed( 4 ), iwork( * ), mm( * ), nn( * )
347 REAL e( * ), rwork( * ), s( * ), ssav( * )
348 COMPLEX a( lda, * ), asav( lda, * ), u( ldu, * ),
349 $ usav( ldu, * ), vt( ldvt, * ),
350 $ vtsav( ldvt, * ), work( * )
357 parameter( zero = 0.0e+0, one = 1.0e+0 )
359 parameter( czero = ( 0.0e+0, 0.0e+0 ),
360 $ cone = ( 1.0e+0, 0.0e+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 REAL anorm, dif, div, ovfl, ulp, ulpinv, unfl
387 INTRINSIC abs, max, min, real
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(
'CDRVBD', -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
claset(
'Full', m, n, czero, czero, a, lda )
493 DO 30 i = 1, min( m, n )
497 ELSE IF( jtype.EQ.2 )
THEN
501 CALL
claset(
'Full', m, n, czero, cone, a, lda )
502 DO 40 i = 1, min( m, n )
516 CALL
clatms( m, n,
'U', iseed,
'N', s, 4,
REAL( 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
clacpy(
'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
clacpy(
'F', m, n, asav, lda, a, lda )
550 CALL
cgesvd(
'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
cbdt01( 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
cunt01(
'Columns', mnmin, m, usav, ldu, work,
565 $ lwork, rwork, result( 2 ) )
566 CALL
cunt01(
'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
clacpy(
'F', m, n, asav, lda, a, lda )
593 CALL
cgesvd( 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
cunt03(
'C', m, mnmin, m, mnmin, usav,
602 $ ldu, a, lda, work, lwork, rwork,
604 ELSE IF( iju.EQ.2 )
THEN
605 CALL
cunt03(
'C', m, mnmin, m, mnmin, usav,
606 $ ldu, u, ldu, work, lwork, rwork,
608 ELSE IF( iju.EQ.3 )
THEN
609 CALL
cunt03(
'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
cunt03(
'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
cunt03(
'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
cunt03(
'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(
REAL( mnmin )*ulp*s( 1 ),
640 $
slamch(
'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
clacpy(
'F', m, n, asav, lda, a, lda )
664 CALL
cgesdd(
'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
cbdt01( 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
cunt01(
'Columns', mnmin, m, usav, ldu, work,
679 $ lwork, rwork, result( 9 ) )
680 CALL
cunt01(
'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
clacpy(
'F', m, n, asav, lda, a, lda )
703 CALL
cgesdd( 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
cunt03(
'C', m, mnmin, m, mnmin, usav,
713 $ ldu, a, lda, work, lwork, rwork,
716 CALL
cunt03(
'C', m, mnmin, m, mnmin, usav,
717 $ ldu, u, ldu, work, lwork, rwork,
720 ELSE IF( ijq.EQ.2 )
THEN
721 CALL
cunt03(
'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
cunt03(
'R', n, mnmin, n, mnmin, vtsav,
735 $ ldvt, vt, ldvt, work, lwork,
736 $ rwork, dif, iinfo )
738 CALL
cunt03(
'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
cunt03(
'R', n, mnmin, n, mnmin, vtsav,
744 $ ldvt, vt, ldvt, work, lwork, rwork,
748 result( 13 ) = max( result( 13 ), dif )
753 div = max(
REAL( mnmin )*ulp*s( 1 ),
754 $
slamch(
'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(
'CBD', nounit, nerrs, ntestt, 0 )
803 9999
FORMAT(
' SVD -- Complex Singular Value Decomposition Driver ',
804 $ /
' Matrix types (see CDRVBD 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(
' CDRVBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
834 $ i6,
', N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
836 9995
FORMAT(
' CDRVBD: ', a,
' returned INFO=', i6,
'.', / 9x,
'M=',
837 $ i6,
', N=', i6,
', JTYPE=', i6,
', LSWORK=', i6, / 9x,
838 $
'ISEED=(', 3( i5,
',' ), i5,
')' )
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 cunt03(RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK, RWORK, RESULT, INFO)
CUNT03
subroutine cunt01(ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, RESID)
CUNT01
subroutine cbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RWORK, RESID)
CBDT01
subroutine clatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
CLATMS
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
CGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine cdrvbd(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)
CDRVBD
subroutine cgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO)
CGESDD
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(CMACH)
SLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j