LAPACK: TESTING/EIG/ddrvbd.f Source File

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