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