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