LAPACK: TESTING/EIG/cdrvbd.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
00002      $                   A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
00003      $                   SSAV, E, WORK, LWORK, RWORK, IWORK, NOUNIT,
00004      $                   INFO )
00005 *
00006 *  -- LAPACK test routine (version 3.1) --
00007 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       INTEGER            INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES,
00012      $                   NTYPES
00013       REAL               THRESH
00014 *     ..
00015 *     .. Array Arguments ..
00016       LOGICAL            DOTYPE( * )
00017       INTEGER            ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
00018       REAL               E( * ), RWORK( * ), S( * ), SSAV( * )
00019       COMPLEX            A( LDA, * ), ASAV( LDA, * ), U( LDU, * ),
00020      $                   USAV( LDU, * ), VT( LDVT, * ),
00021      $                   VTSAV( LDVT, * ), WORK( * )
00022 *     ..
00023 *
00024 *  Purpose
00025 *  =======
00026 *
00027 *  CDRVBD checks the singular value decomposition (SVD) driver CGESVD
00028 *  and CGESDD.
00029 *  CGESVD and CGESDD factors A = U diag(S) VT, where U and VT are
00030 *  unitary and diag(S) is diagonal with the entries of the array S on
00031 *  its diagonal. The entries of S are the singular values, nonnegative
00032 *  and stored in decreasing order.  U and VT can be optionally not
00033 *  computed, overwritten on A, or computed partially.
00034 *
00035 *  A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN.
00036 *  U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N.
00037 *
00038 *  When CDRVBD is called, a number of matrix "sizes" (M's and N's)
00039 *  and a number of matrix "types" are specified.  For each size (M,N)
00040 *  and each type of matrix, and for the minimal workspace as well as
00041 *  workspace adequate to permit blocking, an  M x N  matrix "A" will be
00042 *  generated and used to test the SVD routines.  For each matrix, A will
00043 *  be factored as A = U diag(S) VT and the following 12 tests computed:
00044 *
00045 *  Test for CGESVD:
00046 *
00047 *  (1)   | A - U diag(S) VT | / ( |A| max(M,N) ulp )
00048 *
00049 *  (2)   | I - U'U | / ( M ulp )
00050 *
00051 *  (3)   | I - VT VT' | / ( N ulp )
00052 *
00053 *  (4)   S contains MNMIN nonnegative values in decreasing order.
00054 *        (Return 0 if true, 1/ULP if false.)
00055 *
00056 *  (5)   | U - Upartial | / ( M ulp ) where Upartial is a partially
00057 *        computed U.
00058 *
00059 *  (6)   | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
00060 *        computed VT.
00061 *
00062 *  (7)   | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
00063 *        vector of singular values from the partial SVD
00064 *
00065 *  Test for CGESDD:
00066 *
00067 *  (1)   | A - U diag(S) VT | / ( |A| max(M,N) ulp )
00068 *
00069 *  (2)   | I - U'U | / ( M ulp )
00070 *
00071 *  (3)   | I - VT VT' | / ( N ulp )
00072 *
00073 *  (4)   S contains MNMIN nonnegative values in decreasing order.
00074 *        (Return 0 if true, 1/ULP if false.)
00075 *
00076 *  (5)   | U - Upartial | / ( M ulp ) where Upartial is a partially
00077 *        computed U.
00078 *
00079 *  (6)   | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
00080 *        computed VT.
00081 *
00082 *  (7)   | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
00083 *        vector of singular values from the partial SVD
00084 *
00085 *  The "sizes" are specified by the arrays MM(1:NSIZES) and
00086 *  NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
00087 *  specifies one size.  The "types" are specified by a logical array
00088 *  DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j"
00089 *  will be generated.
00090 *  Currently, the list of possible types is:
00091 *
00092 *  (1)  The zero matrix.
00093 *  (2)  The identity matrix.
00094 *  (3)  A matrix of the form  U D V, where U and V are unitary and
00095 *       D has evenly spaced entries 1, ..., ULP with random signs
00096 *       on the diagonal.
00097 *  (4)  Same as (3), but multiplied by the underflow-threshold / ULP.
00098 *  (5)  Same as (3), but multiplied by the overflow-threshold * ULP.
00099 *
00100 *  Arguments
00101 *  ==========
00102 *
00103 *  NSIZES  (input) INTEGER
00104 *          The number of sizes of matrices to use.  If it is zero,
00105 *          CDRVBD does nothing.  It must be at least zero.
00106 *
00107 *  MM      (input) INTEGER array, dimension (NSIZES)
00108 *          An array containing the matrix "heights" to be used.  For
00109 *          each j=1,...,NSIZES, if MM(j) is zero, then MM(j) and NN(j)
00110 *          will be ignored.  The MM(j) values must be at least zero.
00111 *
00112 *  NN      (input) INTEGER array, dimension (NSIZES)
00113 *          An array containing the matrix "widths" to be used.  For
00114 *          each j=1,...,NSIZES, if NN(j) is zero, then MM(j) and NN(j)
00115 *          will be ignored.  The NN(j) values must be at least zero.
00116 *
00117 *  NTYPES  (input) INTEGER
00118 *          The number of elements in DOTYPE.   If it is zero, CDRVBD
00119 *          does nothing.  It must be at least zero.  If it is MAXTYP+1
00120 *          and NSIZES is 1, then an additional type, MAXTYP+1 is
00121 *          defined, which is to use whatever matrices are in A and B.
00122 *          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
00123 *          DOTYPE(MAXTYP+1) is .TRUE. .
00124 *
00125 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
00126 *          If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
00127 *          of type j will be generated.  If NTYPES is smaller than the
00128 *          maximum number of types defined (PARAMETER MAXTYP), then
00129 *          types NTYPES+1 through MAXTYP will not be generated.  If
00130 *          NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
00131 *          DOTYPE(NTYPES) will be ignored.
00132 *
00133 *  ISEED   (input/output) INTEGER array, dimension (4)
00134 *          On entry ISEED specifies the seed of the random number
00135 *          generator. The array elements should be between 0 and 4095;
00136 *          if not they will be reduced mod 4096.  Also, ISEED(4) must
00137 *          be odd.  The random number generator uses a linear
00138 *          congruential sequence limited to small integers, and so
00139 *          should produce machine independent random numbers. The
00140 *          values of ISEED are changed on exit, and can be used in the
00141 *          next call to CDRVBD to continue the same random number
00142 *          sequence.
00143 *
00144 *  THRESH  (input) REAL
00145 *          A test will count as "failed" if the "error", computed as
00146 *          described above, exceeds THRESH.  Note that the error
00147 *          is scaled to be O(1), so THRESH should be a reasonably
00148 *          small multiple of 1, e.g., 10 or 100.  In particular,
00149 *          it should not depend on the precision (single vs. double)
00150 *          or the size of the matrix.  It must be at least zero.
00151 *
00152 *  NOUNIT  (input) INTEGER
00153 *          The FORTRAN unit number for printing out error messages
00154 *          (e.g., if a routine returns IINFO not equal to 0.)
00155 *
00156 *  A       (output) COMPLEX array, dimension (LDA,max(NN))
00157 *          Used to hold the matrix whose singular values are to be
00158 *          computed.  On exit, A contains the last matrix actually
00159 *          used.
00160 *
00161 *  LDA     (input) INTEGER
00162 *          The leading dimension of A.  It must be at
00163 *          least 1 and at least max( MM ).
00164 *
00165 *  U       (output) COMPLEX array, dimension (LDU,max(MM))
00166 *          Used to hold the computed matrix of right singular vectors.
00167 *          On exit, U contains the last such vectors actually computed.
00168 *
00169 *  LDU     (input) INTEGER
00170 *          The leading dimension of U.  It must be at
00171 *          least 1 and at least max( MM ).
00172 *
00173 *  VT      (output) COMPLEX array, dimension (LDVT,max(NN))
00174 *          Used to hold the computed matrix of left singular vectors.
00175 *          On exit, VT contains the last such vectors actually computed.
00176 *
00177 *  LDVT    (input) INTEGER
00178 *          The leading dimension of VT.  It must be at
00179 *          least 1 and at least max( NN ).
00180 *
00181 *  ASAV    (output) COMPLEX array, dimension (LDA,max(NN))
00182 *          Used to hold a different copy of the matrix whose singular
00183 *          values are to be computed.  On exit, A contains the last
00184 *          matrix actually used.
00185 *
00186 *  USAV    (output) COMPLEX array, dimension (LDU,max(MM))
00187 *          Used to hold a different copy of the computed matrix of
00188 *          right singular vectors. On exit, USAV contains the last such
00189 *          vectors actually computed.
00190 *
00191 *  VTSAV   (output) COMPLEX array, dimension (LDVT,max(NN))
00192 *          Used to hold a different copy of the computed matrix of
00193 *          left singular vectors. On exit, VTSAV contains the last such
00194 *          vectors actually computed.
00195 *
00196 *  S       (output) REAL array, dimension (max(min(MM,NN)))
00197 *          Contains the computed singular values.
00198 *
00199 *  SSAV    (output) REAL array, dimension (max(min(MM,NN)))
00200 *          Contains another copy of the computed singular values.
00201 *
00202 *  E       (output) REAL array, dimension (max(min(MM,NN)))
00203 *          Workspace for CGESVD.
00204 *
00205 *  WORK    (workspace) COMPLEX array, dimension (LWORK)
00206 *
00207 *  LWORK   (input) INTEGER
00208 *          The number of entries in WORK.  This must be at least
00209 *          MAX(3*MIN(M,N)+MAX(M,N)**2,5*MIN(M,N),3*MAX(M,N)) for all
00210 *          pairs  (M,N)=(MM(j),NN(j))
00211 *
00212 *  RWORK   (workspace) REAL array,
00213 *                      dimension ( 5*max(max(MM,NN)) )
00214 *
00215 *  IWORK   (workspace) INTEGER array, dimension at least 8*min(M,N)
00216 *
00217 *  RESULT  (output) REAL array, dimension (7)
00218 *          The values computed by the 7 tests described above.
00219 *          The values are currently limited to 1/ULP, to avoid
00220 *          overflow.
00221 *
00222 *  INFO    (output) INTEGER
00223 *          If 0, then everything ran OK.
00224 *           -1: NSIZES < 0
00225 *           -2: Some MM(j) < 0
00226 *           -3: Some NN(j) < 0
00227 *           -4: NTYPES < 0
00228 *           -7: THRESH < 0
00229 *          -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
00230 *          -12: LDU < 1 or LDU < MMAX.
00231 *          -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ).
00232 *          -21: LWORK too small.
00233 *          If  CLATMS, or CGESVD returns an error code, the
00234 *              absolute value of it is returned.
00235 *
00236 *  =====================================================================
00237 *
00238 *     .. Parameters ..
00239       REAL               ZERO, ONE
00240       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00241       COMPLEX            CZERO, CONE
00242       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00243      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00244       INTEGER            MAXTYP
00245       PARAMETER          ( MAXTYP = 5 )
00246 *     ..
00247 *     .. Local Scalars ..
00248       LOGICAL            BADMM, BADNN
00249       CHARACTER          JOBQ, JOBU, JOBVT
00250       INTEGER            I, IINFO, IJQ, IJU, IJVT, IWSPC, IWTMP, J,
00251      $                   JSIZE, JTYPE, LSWORK, M, MINWRK, MMAX, MNMAX,
00252      $                   MNMIN, MTYPES, N, NERRS, NFAIL, NMAX, NTEST,
00253      $                   NTESTF, NTESTT
00254       REAL               ANORM, DIF, DIV, OVFL, ULP, ULPINV, UNFL
00255 *     ..
00256 *     .. Local Arrays ..
00257       CHARACTER          CJOB( 4 )
00258       INTEGER            IOLDSD( 4 )
00259       REAL               RESULT( 14 )
00260 *     ..
00261 *     .. External Functions ..
00262       REAL               SLAMCH
00263       EXTERNAL           SLAMCH
00264 *     ..
00265 *     .. External Subroutines ..
00266       EXTERNAL           ALASVM, CBDT01, CGESDD, CGESVD, CLACPY, CLASET,
00267      $                   CLATMS, CUNT01, CUNT03, XERBLA
00268 *     ..
00269 *     .. Intrinsic Functions ..
00270       INTRINSIC          ABS, MAX, MIN, REAL
00271 *     ..
00272 *     .. Data statements ..
00273       DATA               CJOB / 'N', 'O', 'S', 'A' /
00274 *     ..
00275 *     .. Executable Statements ..
00276 *
00277 *     Check for errors
00278 *
00279       INFO = 0
00280 *
00281 *     Important constants
00282 *
00283       NERRS = 0
00284       NTESTT = 0
00285       NTESTF = 0
00286       BADMM = .FALSE.
00287       BADNN = .FALSE.
00288       MMAX = 1
00289       NMAX = 1
00290       MNMAX = 1
00291       MINWRK = 1
00292       DO 10 J = 1, NSIZES
00293          MMAX = MAX( MMAX, MM( J ) )
00294          IF( MM( J ).LT.0 )
00295      $      BADMM = .TRUE.
00296          NMAX = MAX( NMAX, NN( J ) )
00297          IF( NN( J ).LT.0 )
00298      $      BADNN = .TRUE.
00299          MNMAX = MAX( MNMAX, MIN( MM( J ), NN( J ) ) )
00300          MINWRK = MAX( MINWRK, MAX( 3*MIN( MM( J ),
00301      $            NN( J ) )+MAX( MM( J ), NN( J ) )**2, 5*MIN( MM( J ),
00302      $            NN( J ) ), 3*MAX( MM( J ), NN( J ) ) ) )
00303    10 CONTINUE
00304 *
00305 *     Check for errors
00306 *
00307       IF( NSIZES.LT.0 ) THEN
00308          INFO = -1
00309       ELSE IF( BADMM ) THEN
00310          INFO = -2
00311       ELSE IF( BADNN ) THEN
00312          INFO = -3
00313       ELSE IF( NTYPES.LT.0 ) THEN
00314          INFO = -4
00315       ELSE IF( LDA.LT.MAX( 1, MMAX ) ) THEN
00316          INFO = -10
00317       ELSE IF( LDU.LT.MAX( 1, MMAX ) ) THEN
00318          INFO = -12
00319       ELSE IF( LDVT.LT.MAX( 1, NMAX ) ) THEN
00320          INFO = -14
00321       ELSE IF( MINWRK.GT.LWORK ) THEN
00322          INFO = -21
00323       END IF
00324 *
00325       IF( INFO.NE.0 ) THEN
00326          CALL XERBLA( 'CDRVBD', -INFO )
00327          RETURN
00328       END IF
00329 *
00330 *     Quick return if nothing to do
00331 *
00332       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
00333      $   RETURN
00334 *
00335 *     More Important constants
00336 *
00337       UNFL = SLAMCH( 'S' )
00338       OVFL = ONE / UNFL
00339       ULP = SLAMCH( 'E' )
00340       ULPINV = ONE / ULP
00341 *
00342 *     Loop over sizes, types
00343 *
00344       NERRS = 0
00345 *
00346       DO 180 JSIZE = 1, NSIZES
00347          M = MM( JSIZE )
00348          N = NN( JSIZE )
00349          MNMIN = MIN( M, N )
00350 *
00351          IF( NSIZES.NE.1 ) THEN
00352             MTYPES = MIN( MAXTYP, NTYPES )
00353          ELSE
00354             MTYPES = MIN( MAXTYP+1, NTYPES )
00355          END IF
00356 *
00357          DO 170 JTYPE = 1, MTYPES
00358             IF( .NOT.DOTYPE( JTYPE ) )
00359      $         GO TO 170
00360             NTEST = 0
00361 *
00362             DO 20 J = 1, 4
00363                IOLDSD( J ) = ISEED( J )
00364    20       CONTINUE
00365 *
00366 *           Compute "A"
00367 *
00368             IF( MTYPES.GT.MAXTYP )
00369      $         GO TO 50
00370 *
00371             IF( JTYPE.EQ.1 ) THEN
00372 *
00373 *              Zero matrix
00374 *
00375                CALL CLASET( 'Full', M, N, CZERO, CZERO, A, LDA )
00376                DO 30 I = 1, MIN( M, N )
00377                   S( I ) = ZERO
00378    30          CONTINUE
00379 *
00380             ELSE IF( JTYPE.EQ.2 ) THEN
00381 *
00382 *              Identity matrix
00383 *
00384                CALL CLASET( 'Full', M, N, CZERO, CONE, A, LDA )
00385                DO 40 I = 1, MIN( M, N )
00386                   S( I ) = ONE
00387    40          CONTINUE
00388 *
00389             ELSE
00390 *
00391 *              (Scaled) random matrix
00392 *
00393                IF( JTYPE.EQ.3 )
00394      $            ANORM = ONE
00395                IF( JTYPE.EQ.4 )
00396      $            ANORM = UNFL / ULP
00397                IF( JTYPE.EQ.5 )
00398      $            ANORM = OVFL*ULP
00399                CALL CLATMS( M, N, 'U', ISEED, 'N', S, 4, REAL( MNMIN ),
00400      $                      ANORM, M-1, N-1, 'N', A, LDA, WORK, IINFO )
00401                IF( IINFO.NE.0 ) THEN
00402                   WRITE( NOUNIT, FMT = 9996 )'Generator', IINFO, M, N,
00403      $               JTYPE, IOLDSD
00404                   INFO = ABS( IINFO )
00405                   RETURN
00406                END IF
00407             END IF
00408 *
00409    50       CONTINUE
00410             CALL CLACPY( 'F', M, N, A, LDA, ASAV, LDA )
00411 *
00412 *           Do for minimal and adequate (for blocking) workspace
00413 *
00414             DO 160 IWSPC = 1, 4
00415 *
00416 *              Test for CGESVD
00417 *
00418                IWTMP = 2*MIN( M, N )+MAX( M, N )
00419                LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
00420                LSWORK = MIN( LSWORK, LWORK )
00421                LSWORK = MAX( LSWORK, 1 )
00422                IF( IWSPC.EQ.4 )
00423      $            LSWORK = LWORK
00424 *
00425                DO 60 J = 1, 14
00426                   RESULT( J ) = -ONE
00427    60          CONTINUE
00428 *
00429 *              Factorize A
00430 *
00431                IF( IWSPC.GT.1 )
00432      $            CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
00433                CALL CGESVD( 'A', 'A', M, N, A, LDA, SSAV, USAV, LDU,
00434      $                      VTSAV, LDVT, WORK, LSWORK, RWORK, IINFO )
00435                IF( IINFO.NE.0 ) THEN
00436                   WRITE( NOUNIT, FMT = 9995 )'GESVD', IINFO, M, N,
00437      $               JTYPE, LSWORK, IOLDSD
00438                   INFO = ABS( IINFO )
00439                   RETURN
00440                END IF
00441 *
00442 *              Do tests 1--4
00443 *
00444                CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
00445      $                      VTSAV, LDVT, WORK, RWORK, RESULT( 1 ) )
00446                IF( M.NE.0 .AND. N.NE.0 ) THEN
00447                   CALL CUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK,
00448      $                         LWORK, RWORK, RESULT( 2 ) )
00449                   CALL CUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK,
00450      $                         LWORK, RWORK, RESULT( 3 ) )
00451                END IF
00452                RESULT( 4 ) = 0
00453                DO 70 I = 1, MNMIN - 1
00454                   IF( SSAV( I ).LT.SSAV( I+1 ) )
00455      $               RESULT( 4 ) = ULPINV
00456                   IF( SSAV( I ).LT.ZERO )
00457      $               RESULT( 4 ) = ULPINV
00458    70          CONTINUE
00459                IF( MNMIN.GE.1 ) THEN
00460                   IF( SSAV( MNMIN ).LT.ZERO )
00461      $               RESULT( 4 ) = ULPINV
00462                END IF
00463 *
00464 *              Do partial SVDs, comparing to SSAV, USAV, and VTSAV
00465 *
00466                RESULT( 5 ) = ZERO
00467                RESULT( 6 ) = ZERO
00468                RESULT( 7 ) = ZERO
00469                DO 100 IJU = 0, 3
00470                   DO 90 IJVT = 0, 3
00471                      IF( ( IJU.EQ.3 .AND. IJVT.EQ.3 ) .OR.
00472      $                   ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 90
00473                      JOBU = CJOB( IJU+1 )
00474                      JOBVT = CJOB( IJVT+1 )
00475                      CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
00476                      CALL CGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
00477      $                            VT, LDVT, WORK, LSWORK, RWORK, IINFO )
00478 *
00479 *                    Compare U
00480 *
00481                      DIF = ZERO
00482                      IF( M.GT.0 .AND. N.GT.0 ) THEN
00483                         IF( IJU.EQ.1 ) THEN
00484                            CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
00485      $                                  LDU, A, LDA, WORK, LWORK, RWORK,
00486      $                                  DIF, IINFO )
00487                         ELSE IF( IJU.EQ.2 ) THEN
00488                            CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
00489      $                                  LDU, U, LDU, WORK, LWORK, RWORK,
00490      $                                  DIF, IINFO )
00491                         ELSE IF( IJU.EQ.3 ) THEN
00492                            CALL CUNT03( 'C', M, M, M, MNMIN, USAV, LDU,
00493      $                                  U, LDU, WORK, LWORK, RWORK, DIF,
00494      $                                  IINFO )
00495                         END IF
00496                      END IF
00497                      RESULT( 5 ) = MAX( RESULT( 5 ), DIF )
00498 *
00499 *                    Compare VT
00500 *
00501                      DIF = ZERO
00502                      IF( M.GT.0 .AND. N.GT.0 ) THEN
00503                         IF( IJVT.EQ.1 ) THEN
00504                            CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
00505      $                                  LDVT, A, LDA, WORK, LWORK,
00506      $                                  RWORK, DIF, IINFO )
00507                         ELSE IF( IJVT.EQ.2 ) THEN
00508                            CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
00509      $                                  LDVT, VT, LDVT, WORK, LWORK,
00510      $                                  RWORK, DIF, IINFO )
00511                         ELSE IF( IJVT.EQ.3 ) THEN
00512                            CALL CUNT03( 'R', N, N, N, MNMIN, VTSAV,
00513      $                                  LDVT, VT, LDVT, WORK, LWORK,
00514      $                                  RWORK, DIF, IINFO )
00515                         END IF
00516                      END IF
00517                      RESULT( 6 ) = MAX( RESULT( 6 ), DIF )
00518 *
00519 *                    Compare S
00520 *
00521                      DIF = ZERO
00522                      DIV = MAX( REAL( MNMIN )*ULP*S( 1 ),
00523      $                     SLAMCH( 'Safe minimum' ) )
00524                      DO 80 I = 1, MNMIN - 1
00525                         IF( SSAV( I ).LT.SSAV( I+1 ) )
00526      $                     DIF = ULPINV
00527                         IF( SSAV( I ).LT.ZERO )
00528      $                     DIF = ULPINV
00529                         DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
00530    80                CONTINUE
00531                      RESULT( 7 ) = MAX( RESULT( 7 ), DIF )
00532    90             CONTINUE
00533   100          CONTINUE
00534 *
00535 *              Test for CGESDD
00536 *
00537                IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
00538                LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
00539                LSWORK = MIN( LSWORK, LWORK )
00540                LSWORK = MAX( LSWORK, 1 )
00541                IF( IWSPC.EQ.4 )
00542      $            LSWORK = LWORK
00543 *
00544 *              Factorize A
00545 *
00546                CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
00547                CALL CGESDD( 'A', M, N, A, LDA, SSAV, USAV, LDU, VTSAV,
00548      $                      LDVT, WORK, LSWORK, RWORK, IWORK, IINFO )
00549                IF( IINFO.NE.0 ) THEN
00550                   WRITE( NOUNIT, FMT = 9995 )'GESDD', IINFO, M, N,
00551      $               JTYPE, LSWORK, IOLDSD
00552                   INFO = ABS( IINFO )
00553                   RETURN
00554                END IF
00555 *
00556 *              Do tests 1--4
00557 *
00558                CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
00559      $                      VTSAV, LDVT, WORK, RWORK, RESULT( 8 ) )
00560                IF( M.NE.0 .AND. N.NE.0 ) THEN
00561                   CALL CUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK,
00562      $                         LWORK, RWORK, RESULT( 9 ) )
00563                   CALL CUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK,
00564      $                         LWORK, RWORK, RESULT( 10 ) )
00565                END IF
00566                RESULT( 11 ) = 0
00567                DO 110 I = 1, MNMIN - 1
00568                   IF( SSAV( I ).LT.SSAV( I+1 ) )
00569      $               RESULT( 11 ) = ULPINV
00570                   IF( SSAV( I ).LT.ZERO )
00571      $               RESULT( 11 ) = ULPINV
00572   110          CONTINUE
00573                IF( MNMIN.GE.1 ) THEN
00574                   IF( SSAV( MNMIN ).LT.ZERO )
00575      $               RESULT( 11 ) = ULPINV
00576                END IF
00577 *
00578 *              Do partial SVDs, comparing to SSAV, USAV, and VTSAV
00579 *
00580                RESULT( 12 ) = ZERO
00581                RESULT( 13 ) = ZERO
00582                RESULT( 14 ) = ZERO
00583                DO 130 IJQ = 0, 2
00584                   JOBQ = CJOB( IJQ+1 )
00585                   CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
00586                   CALL CGESDD( JOBQ, M, N, A, LDA, S, U, LDU, VT, LDVT,
00587      $                         WORK, LSWORK, RWORK, IWORK, IINFO )
00588 *
00589 *                 Compare U
00590 *
00591                   DIF = ZERO
00592                   IF( M.GT.0 .AND. N.GT.0 ) THEN
00593                      IF( IJQ.EQ.1 ) THEN
00594                         IF( M.GE.N ) THEN
00595                            CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
00596      $                                  LDU, A, LDA, WORK, LWORK, RWORK,
00597      $                                  DIF, IINFO )
00598                         ELSE
00599                            CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
00600      $                                  LDU, U, LDU, WORK, LWORK, RWORK,
00601      $                                  DIF, IINFO )
00602                         END IF
00603                      ELSE IF( IJQ.EQ.2 ) THEN
00604                         CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV, LDU,
00605      $                               U, LDU, WORK, LWORK, RWORK, DIF,
00606      $                               IINFO )
00607                      END IF
00608                   END IF
00609                   RESULT( 12 ) = MAX( RESULT( 12 ), DIF )
00610 *
00611 *                 Compare VT
00612 *
00613                   DIF = ZERO
00614                   IF( M.GT.0 .AND. N.GT.0 ) THEN
00615                      IF( IJQ.EQ.1 ) THEN
00616                         IF( M.GE.N ) THEN
00617                            CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
00618      $                                  LDVT, VT, LDVT, WORK, LWORK,
00619      $                                  RWORK, DIF, IINFO )
00620                         ELSE
00621                            CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
00622      $                                  LDVT, A, LDA, WORK, LWORK,
00623      $                                  RWORK, DIF, IINFO )
00624                         END IF
00625                      ELSE IF( IJQ.EQ.2 ) THEN
00626                         CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
00627      $                               LDVT, VT, LDVT, WORK, LWORK, RWORK,
00628      $                               DIF, IINFO )
00629                      END IF
00630                   END IF
00631                   RESULT( 13 ) = MAX( RESULT( 13 ), DIF )
00632 *
00633 *                 Compare S
00634 *
00635                   DIF = ZERO
00636                   DIV = MAX( REAL( MNMIN )*ULP*S( 1 ),
00637      $                  SLAMCH( 'Safe minimum' ) )
00638                   DO 120 I = 1, MNMIN - 1
00639                      IF( SSAV( I ).LT.SSAV( I+1 ) )
00640      $                  DIF = ULPINV
00641                      IF( SSAV( I ).LT.ZERO )
00642      $                  DIF = ULPINV
00643                      DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
00644   120             CONTINUE
00645                   RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
00646   130          CONTINUE
00647 *
00648 *              End of Loop -- Check for RESULT(j) > THRESH
00649 *
00650                NTEST = 0
00651                NFAIL = 0
00652                DO 140 J = 1, 14
00653                   IF( RESULT( J ).GE.ZERO )
00654      $               NTEST = NTEST + 1
00655                   IF( RESULT( J ).GE.THRESH )
00656      $               NFAIL = NFAIL + 1
00657   140          CONTINUE
00658 *
00659                IF( NFAIL.GT.0 )
00660      $            NTESTF = NTESTF + 1
00661                IF( NTESTF.EQ.1 ) THEN
00662                   WRITE( NOUNIT, FMT = 9999 )
00663                   WRITE( NOUNIT, FMT = 9998 )THRESH
00664                   NTESTF = 2
00665                END IF
00666 *
00667                DO 150 J = 1, 14
00668                   IF( RESULT( J ).GE.THRESH ) THEN
00669                      WRITE( NOUNIT, FMT = 9997 )M, N, JTYPE, IWSPC,
00670      $                  IOLDSD, J, RESULT( J )
00671                   END IF
00672   150          CONTINUE
00673 *
00674                NERRS = NERRS + NFAIL
00675                NTESTT = NTESTT + NTEST
00676 *
00677   160       CONTINUE
00678 *
00679   170    CONTINUE
00680   180 CONTINUE
00681 *
00682 *     Summary
00683 *
00684       CALL ALASVM( 'CBD', NOUNIT, NERRS, NTESTT, 0 )
00685 *
00686  9999 FORMAT( ' SVD -- Complex Singular Value Decomposition Driver ',
00687      $      / ' Matrix types (see CDRVBD for details):',
00688      $      / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
00689      $      / ' 3 = Evenly spaced singular values near 1',
00690      $      / ' 4 = Evenly spaced singular values near underflow',
00691      $      / ' 5 = Evenly spaced singular values near overflow',
00692      $      / / ' Tests performed: ( A is dense, U and V are unitary,',
00693      $      / 19X, ' S is an array, and Upartial, VTpartial, and',
00694      $      / 19X, ' Spartial are partially computed U, VT and S),', / )
00695  9998 FORMAT( ' Tests performed with Test Threshold = ', F8.2,
00696      $      / ' CGESVD: ', /
00697      $      ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
00698      $      / ' 2 = | I - U**T U | / ( M ulp ) ',
00699      $      / ' 3 = | I - VT VT**T | / ( N ulp ) ',
00700      $      / ' 4 = 0 if S contains min(M,N) nonnegative values in',
00701      $      ' decreasing order, else 1/ulp',
00702      $      / ' 5 = | U - Upartial | / ( M ulp )',
00703      $      / ' 6 = | VT - VTpartial | / ( N ulp )',
00704      $      / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
00705      $      / ' CGESDD: ', /
00706      $      ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
00707      $      / ' 9 = | I - U**T U | / ( M ulp ) ',
00708      $      / '10 = | I - VT VT**T | / ( N ulp ) ',
00709      $      / '11 = 0 if S contains min(M,N) nonnegative values in',
00710      $      ' decreasing order, else 1/ulp',
00711      $      / '12 = | U - Upartial | / ( M ulp )',
00712      $      / '13 = | VT - VTpartial | / ( N ulp )',
00713      $      / '14 = | S - Spartial | / ( min(M,N) ulp |S| )', / / )
00714  9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
00715      $      ', seed=', 4( I4, ',' ), ' test(', I1, ')=', G11.4 )
00716  9996 FORMAT( ' CDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
00717      $      I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
00718      $      I5, ')' )
00719  9995 FORMAT( ' CDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
00720      $      I6, ', N=', I6, ', JTYPE=', I6, ', LSWORK=', I6, / 9X,
00721      $      'ISEED=(', 3( I5, ',' ), I5, ')' )
00722 *
00723       RETURN
00724 *
00725 *     End of CDRVBD
00726 *
00727       END