LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zchkbd.f
Go to the documentation of this file.
00001 *> \brief \b ZCHKBD
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 ZCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
00012 *                          ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX,
00013 *                          Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK,
00014 *                          RWORK, NOUT, INFO )
00015 * 
00016 *       .. Scalar Arguments ..
00017 *       INTEGER            INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
00018 *      $                   NSIZES, NTYPES
00019 *       DOUBLE PRECISION   THRESH
00020 *       ..
00021 *       .. Array Arguments ..
00022 *       LOGICAL            DOTYPE( * )
00023 *       INTEGER            ISEED( 4 ), MVAL( * ), NVAL( * )
00024 *       DOUBLE PRECISION   BD( * ), BE( * ), RWORK( * ), S1( * ), S2( * )
00025 *       COMPLEX*16         A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
00026 *      $                   U( LDPT, * ), VT( LDPT, * ), WORK( * ),
00027 *      $                   X( LDX, * ), Y( LDX, * ), Z( LDX, * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *> ZCHKBD checks the singular value decomposition (SVD) routines.
00037 *>
00038 *> ZGEBRD reduces a complex general m by n matrix A to real upper or
00039 *> lower bidiagonal form by an orthogonal transformation: Q' * A * P = B
00040 *> (or A = Q * B * P').  The matrix B is upper bidiagonal if m >= n
00041 *> and lower bidiagonal if m < n.
00042 *>
00043 *> ZUNGBR generates the orthogonal matrices Q and P' from ZGEBRD.
00044 *> Note that Q and P are not necessarily square.
00045 *>
00046 *> ZBDSQR computes the singular value decomposition of the bidiagonal
00047 *> matrix B as B = U S V'.  It is called three times to compute
00048 *>    1)  B = U S1 V', where S1 is the diagonal matrix of singular
00049 *>        values and the columns of the matrices U and V are the left
00050 *>        and right singular vectors, respectively, of B.
00051 *>    2)  Same as 1), but the singular values are stored in S2 and the
00052 *>        singular vectors are not computed.
00053 *>    3)  A = (UQ) S (P'V'), the SVD of the original matrix A.
00054 *> In addition, ZBDSQR has an option to apply the left orthogonal matrix
00055 *> U to a matrix X, useful in least squares applications.
00056 *>
00057 *> For each pair of matrix dimensions (M,N) and each selected matrix
00058 *> type, an M by N matrix A and an M by NRHS matrix X are generated.
00059 *> The problem dimensions are as follows
00060 *>    A:          M x N
00061 *>    Q:          M x min(M,N) (but M x M if NRHS > 0)
00062 *>    P:          min(M,N) x N
00063 *>    B:          min(M,N) x min(M,N)
00064 *>    U, V:       min(M,N) x min(M,N)
00065 *>    S1, S2      diagonal, order min(M,N)
00066 *>    X:          M x NRHS
00067 *>
00068 *> For each generated matrix, 14 tests are performed:
00069 *>
00070 *> Test ZGEBRD and ZUNGBR
00071 *>
00072 *> (1)   | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P'
00073 *>
00074 *> (2)   | I - Q' Q | / ( M ulp )
00075 *>
00076 *> (3)   | I - PT PT' | / ( N ulp )
00077 *>
00078 *> Test ZBDSQR on bidiagonal matrix B
00079 *>
00080 *> (4)   | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
00081 *>
00082 *> (5)   | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X
00083 *>                                                  and   Z = U' Y.
00084 *> (6)   | I - U' U | / ( min(M,N) ulp )
00085 *>
00086 *> (7)   | I - VT VT' | / ( min(M,N) ulp )
00087 *>
00088 *> (8)   S1 contains min(M,N) nonnegative values in decreasing order.
00089 *>       (Return 0 if true, 1/ULP if false.)
00090 *>
00091 *> (9)   0 if the true singular values of B are within THRESH of
00092 *>       those in S1.  2*THRESH if they are not.  (Tested using
00093 *>       DSVDCH)
00094 *>
00095 *> (10)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
00096 *>                                   computing U and V.
00097 *>
00098 *> Test ZBDSQR on matrix A
00099 *>
00100 *> (11)  | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp )
00101 *>
00102 *> (12)  | X - (QU) Z | / ( |X| max(M,k) ulp )
00103 *>
00104 *> (13)  | I - (QU)'(QU) | / ( M ulp )
00105 *>
00106 *> (14)  | I - (VT PT) (PT'VT') | / ( N ulp )
00107 *>
00108 *> The possible matrix types are
00109 *>
00110 *> (1)  The zero matrix.
00111 *> (2)  The identity matrix.
00112 *>
00113 *> (3)  A diagonal matrix with evenly spaced entries
00114 *>      1, ..., ULP  and random signs.
00115 *>      (ULP = (first number larger than 1) - 1 )
00116 *> (4)  A diagonal matrix with geometrically spaced entries
00117 *>      1, ..., ULP  and random signs.
00118 *> (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
00119 *>      and random signs.
00120 *>
00121 *> (6)  Same as (3), but multiplied by SQRT( overflow threshold )
00122 *> (7)  Same as (3), but multiplied by SQRT( underflow threshold )
00123 *>
00124 *> (8)  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 *>
00128 *> (9)  A matrix of the form  U D V, where U and V are orthogonal and
00129 *>      D has geometrically spaced entries 1, ..., ULP with random
00130 *>      signs on the diagonal.
00131 *>
00132 *> (10) A matrix of the form  U D V, where U and V are orthogonal and
00133 *>      D has "clustered" entries 1, ULP,..., ULP with random
00134 *>      signs on the diagonal.
00135 *>
00136 *> (11) Same as (8), but multiplied by SQRT( overflow threshold )
00137 *> (12) Same as (8), but multiplied by SQRT( underflow threshold )
00138 *>
00139 *> (13) Rectangular matrix with random entries chosen from (-1,1).
00140 *> (14) Same as (13), but multiplied by SQRT( overflow threshold )
00141 *> (15) Same as (13), but multiplied by SQRT( underflow threshold )
00142 *>
00143 *> Special case:
00144 *> (16) A bidiagonal matrix with random entries chosen from a
00145 *>      logarithmic distribution on [ulp^2,ulp^(-2)]  (I.e., each
00146 *>      entry is  e^x, where x is chosen uniformly on
00147 *>      [ 2 log(ulp), -2 log(ulp) ] .)  For *this* type:
00148 *>      (a) ZGEBRD is not called to reduce it to bidiagonal form.
00149 *>      (b) the bidiagonal is  min(M,N) x min(M,N); if M<N, the
00150 *>          matrix will be lower bidiagonal, otherwise upper.
00151 *>      (c) only tests 5--8 and 14 are performed.
00152 *>
00153 *> A subset of the full set of matrix types may be selected through
00154 *> the logical array DOTYPE.
00155 *> \endverbatim
00156 *
00157 *  Arguments:
00158 *  ==========
00159 *
00160 *> \param[in] NSIZES
00161 *> \verbatim
00162 *>          NSIZES is INTEGER
00163 *>          The number of values of M and N contained in the vectors
00164 *>          MVAL and NVAL.  The matrix sizes are used in pairs (M,N).
00165 *> \endverbatim
00166 *>
00167 *> \param[in] MVAL
00168 *> \verbatim
00169 *>          MVAL is INTEGER array, dimension (NM)
00170 *>          The values of the matrix row dimension M.
00171 *> \endverbatim
00172 *>
00173 *> \param[in] NVAL
00174 *> \verbatim
00175 *>          NVAL is INTEGER array, dimension (NM)
00176 *>          The values of the matrix column dimension N.
00177 *> \endverbatim
00178 *>
00179 *> \param[in] NTYPES
00180 *> \verbatim
00181 *>          NTYPES is INTEGER
00182 *>          The number of elements in DOTYPE.   If it is zero, ZCHKBD
00183 *>          does nothing.  It must be at least zero.  If it is MAXTYP+1
00184 *>          and NSIZES is 1, then an additional type, MAXTYP+1 is
00185 *>          defined, which is to use whatever matrices are in A and B.
00186 *>          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
00187 *>          DOTYPE(MAXTYP+1) is .TRUE. .
00188 *> \endverbatim
00189 *>
00190 *> \param[in] DOTYPE
00191 *> \verbatim
00192 *>          DOTYPE is LOGICAL array, dimension (NTYPES)
00193 *>          If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
00194 *>          of type j will be generated.  If NTYPES is smaller than the
00195 *>          maximum number of types defined (PARAMETER MAXTYP), then
00196 *>          types NTYPES+1 through MAXTYP will not be generated.  If
00197 *>          NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
00198 *>          DOTYPE(NTYPES) will be ignored.
00199 *> \endverbatim
00200 *>
00201 *> \param[in] NRHS
00202 *> \verbatim
00203 *>          NRHS is INTEGER
00204 *>          The number of columns in the "right-hand side" matrices X, Y,
00205 *>          and Z, used in testing ZBDSQR.  If NRHS = 0, then the
00206 *>          operations on the right-hand side will not be tested.
00207 *>          NRHS must be at least 0.
00208 *> \endverbatim
00209 *>
00210 *> \param[in,out] ISEED
00211 *> \verbatim
00212 *>          ISEED is INTEGER array, dimension (4)
00213 *>          On entry ISEED specifies the seed of the random number
00214 *>          generator. The array elements should be between 0 and 4095;
00215 *>          if not they will be reduced mod 4096.  Also, ISEED(4) must
00216 *>          be odd.  The values of ISEED are changed on exit, and can be
00217 *>          used in the next call to ZCHKBD to continue the same random
00218 *>          number sequence.
00219 *> \endverbatim
00220 *>
00221 *> \param[in] THRESH
00222 *> \verbatim
00223 *>          THRESH is DOUBLE PRECISION
00224 *>          The threshold value for the test ratios.  A result is
00225 *>          included in the output file if RESULT >= THRESH.  To have
00226 *>          every test ratio printed, use THRESH = 0.  Note that the
00227 *>          expected value of the test ratios is O(1), so THRESH should
00228 *>          be a reasonably small multiple of 1, e.g., 10 or 100.
00229 *> \endverbatim
00230 *>
00231 *> \param[out] A
00232 *> \verbatim
00233 *>          A is COMPLEX*16 array, dimension (LDA,NMAX)
00234 *>          where NMAX is the maximum value of N in NVAL.
00235 *> \endverbatim
00236 *>
00237 *> \param[in] LDA
00238 *> \verbatim
00239 *>          LDA is INTEGER
00240 *>          The leading dimension of the array A.  LDA >= max(1,MMAX),
00241 *>          where MMAX is the maximum value of M in MVAL.
00242 *> \endverbatim
00243 *>
00244 *> \param[out] BD
00245 *> \verbatim
00246 *>          BD is DOUBLE PRECISION array, dimension
00247 *>                      (max(min(MVAL(j),NVAL(j))))
00248 *> \endverbatim
00249 *>
00250 *> \param[out] BE
00251 *> \verbatim
00252 *>          BE is DOUBLE PRECISION array, dimension
00253 *>                      (max(min(MVAL(j),NVAL(j))))
00254 *> \endverbatim
00255 *>
00256 *> \param[out] S1
00257 *> \verbatim
00258 *>          S1 is DOUBLE PRECISION array, dimension
00259 *>                      (max(min(MVAL(j),NVAL(j))))
00260 *> \endverbatim
00261 *>
00262 *> \param[out] S2
00263 *> \verbatim
00264 *>          S2 is DOUBLE PRECISION array, dimension
00265 *>                      (max(min(MVAL(j),NVAL(j))))
00266 *> \endverbatim
00267 *>
00268 *> \param[out] X
00269 *> \verbatim
00270 *>          X is COMPLEX*16 array, dimension (LDX,NRHS)
00271 *> \endverbatim
00272 *>
00273 *> \param[in] LDX
00274 *> \verbatim
00275 *>          LDX is INTEGER
00276 *>          The leading dimension of the arrays X, Y, and Z.
00277 *>          LDX >= max(1,MMAX).
00278 *> \endverbatim
00279 *>
00280 *> \param[out] Y
00281 *> \verbatim
00282 *>          Y is COMPLEX*16 array, dimension (LDX,NRHS)
00283 *> \endverbatim
00284 *>
00285 *> \param[out] Z
00286 *> \verbatim
00287 *>          Z is COMPLEX*16 array, dimension (LDX,NRHS)
00288 *> \endverbatim
00289 *>
00290 *> \param[out] Q
00291 *> \verbatim
00292 *>          Q is COMPLEX*16 array, dimension (LDQ,MMAX)
00293 *> \endverbatim
00294 *>
00295 *> \param[in] LDQ
00296 *> \verbatim
00297 *>          LDQ is INTEGER
00298 *>          The leading dimension of the array Q.  LDQ >= max(1,MMAX).
00299 *> \endverbatim
00300 *>
00301 *> \param[out] PT
00302 *> \verbatim
00303 *>          PT is COMPLEX*16 array, dimension (LDPT,NMAX)
00304 *> \endverbatim
00305 *>
00306 *> \param[in] LDPT
00307 *> \verbatim
00308 *>          LDPT is INTEGER
00309 *>          The leading dimension of the arrays PT, U, and V.
00310 *>          LDPT >= max(1, max(min(MVAL(j),NVAL(j)))).
00311 *> \endverbatim
00312 *>
00313 *> \param[out] U
00314 *> \verbatim
00315 *>          U is COMPLEX*16 array, dimension
00316 *>                      (LDPT,max(min(MVAL(j),NVAL(j))))
00317 *> \endverbatim
00318 *>
00319 *> \param[out] VT
00320 *> \verbatim
00321 *>          VT is COMPLEX*16 array, dimension
00322 *>                      (LDPT,max(min(MVAL(j),NVAL(j))))
00323 *> \endverbatim
00324 *>
00325 *> \param[out] WORK
00326 *> \verbatim
00327 *>          WORK is COMPLEX*16 array, dimension (LWORK)
00328 *> \endverbatim
00329 *>
00330 *> \param[in] LWORK
00331 *> \verbatim
00332 *>          LWORK is INTEGER
00333 *>          The number of entries in WORK.  This must be at least
00334 *>          3(M+N) and  M(M + max(M,N,k) + 1) + N*min(M,N)  for all
00335 *>          pairs  (M,N)=(MM(j),NN(j))
00336 *> \endverbatim
00337 *>
00338 *> \param[out] RWORK
00339 *> \verbatim
00340 *>          RWORK is DOUBLE PRECISION array, dimension
00341 *>                      (5*max(min(M,N)))
00342 *> \endverbatim
00343 *>
00344 *> \param[in] NOUT
00345 *> \verbatim
00346 *>          NOUT is INTEGER
00347 *>          The FORTRAN unit number for printing out error messages
00348 *>          (e.g., if a routine returns IINFO not equal to 0.)
00349 *> \endverbatim
00350 *>
00351 *> \param[out] INFO
00352 *> \verbatim
00353 *>          INFO is INTEGER
00354 *>          If 0, then everything ran OK.
00355 *>           -1: NSIZES < 0
00356 *>           -2: Some MM(j) < 0
00357 *>           -3: Some NN(j) < 0
00358 *>           -4: NTYPES < 0
00359 *>           -6: NRHS  < 0
00360 *>           -8: THRESH < 0
00361 *>          -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
00362 *>          -17: LDB < 1 or LDB < MMAX.
00363 *>          -21: LDQ < 1 or LDQ < MMAX.
00364 *>          -23: LDP < 1 or LDP < MNMAX.
00365 *>          -27: LWORK too small.
00366 *>          If  ZLATMR, CLATMS, ZGEBRD, ZUNGBR, or ZBDSQR,
00367 *>              returns an error code, the
00368 *>              absolute value of it is returned.
00369 *>
00370 *>-----------------------------------------------------------------------
00371 *>
00372 *>     Some Local Variables and Parameters:
00373 *>     ---- ----- --------- --- ----------
00374 *>
00375 *>     ZERO, ONE       Real 0 and 1.
00376 *>     MAXTYP          The number of types defined.
00377 *>     NTEST           The number of tests performed, or which can
00378 *>                     be performed so far, for the current matrix.
00379 *>     MMAX            Largest value in NN.
00380 *>     NMAX            Largest value in NN.
00381 *>     MNMIN           min(MM(j), NN(j)) (the dimension of the bidiagonal
00382 *>                     matrix.)
00383 *>     MNMAX           The maximum value of MNMIN for j=1,...,NSIZES.
00384 *>     NFAIL           The number of tests which have exceeded THRESH
00385 *>     COND, IMODE     Values to be passed to the matrix generators.
00386 *>     ANORM           Norm of A; passed to matrix generators.
00387 *>
00388 *>     OVFL, UNFL      Overflow and underflow thresholds.
00389 *>     RTOVFL, RTUNFL  Square roots of the previous 2 values.
00390 *>     ULP, ULPINV     Finest relative precision and its inverse.
00391 *>
00392 *>             The following four arrays decode JTYPE:
00393 *>     KTYPE(j)        The general type (1-10) for type "j".
00394 *>     KMODE(j)        The MODE value to be passed to the matrix
00395 *>                     generator for type "j".
00396 *>     KMAGN(j)        The order of magnitude ( O(1),
00397 *>                     O(overflow^(1/2) ), O(underflow^(1/2) )
00398 *> \endverbatim
00399 *
00400 *  Authors:
00401 *  ========
00402 *
00403 *> \author Univ. of Tennessee 
00404 *> \author Univ. of California Berkeley 
00405 *> \author Univ. of Colorado Denver 
00406 *> \author NAG Ltd. 
00407 *
00408 *> \date November 2011
00409 *
00410 *> \ingroup complex16_eig
00411 *
00412 *  =====================================================================
00413       SUBROUTINE ZCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
00414      $                   ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX,
00415      $                   Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK,
00416      $                   RWORK, NOUT, INFO )
00417 *
00418 *  -- LAPACK test routine (version 3.4.0) --
00419 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00420 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00421 *     November 2011
00422 *
00423 *     .. Scalar Arguments ..
00424       INTEGER            INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
00425      $                   NSIZES, NTYPES
00426       DOUBLE PRECISION   THRESH
00427 *     ..
00428 *     .. Array Arguments ..
00429       LOGICAL            DOTYPE( * )
00430       INTEGER            ISEED( 4 ), MVAL( * ), NVAL( * )
00431       DOUBLE PRECISION   BD( * ), BE( * ), RWORK( * ), S1( * ), S2( * )
00432       COMPLEX*16         A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
00433      $                   U( LDPT, * ), VT( LDPT, * ), WORK( * ),
00434      $                   X( LDX, * ), Y( LDX, * ), Z( LDX, * )
00435 *     ..
00436 *
00437 * ======================================================================
00438 *
00439 *     .. Parameters ..
00440       DOUBLE PRECISION   ZERO, ONE, TWO, HALF
00441       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
00442      $                   HALF = 0.5D0 )
00443       COMPLEX*16         CZERO, CONE
00444       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00445      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00446       INTEGER            MAXTYP
00447       PARAMETER          ( MAXTYP = 16 )
00448 *     ..
00449 *     .. Local Scalars ..
00450       LOGICAL            BADMM, BADNN, BIDIAG
00451       CHARACTER          UPLO
00452       CHARACTER*3        PATH
00453       INTEGER            I, IINFO, IMODE, ITYPE, J, JCOL, JSIZE, JTYPE,
00454      $                   LOG2UI, M, MINWRK, MMAX, MNMAX, MNMIN, MQ,
00455      $                   MTYPES, N, NFAIL, NMAX, NTEST
00456       DOUBLE PRECISION   AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
00457      $                   TEMP1, TEMP2, ULP, ULPINV, UNFL
00458 *     ..
00459 *     .. Local Arrays ..
00460       INTEGER            IOLDSD( 4 ), IWORK( 1 ), KMAGN( MAXTYP ),
00461      $                   KMODE( MAXTYP ), KTYPE( MAXTYP )
00462       DOUBLE PRECISION   DUMMA( 1 ), RESULT( 14 )
00463 *     ..
00464 *     .. External Functions ..
00465       DOUBLE PRECISION   DLAMCH, DLARND
00466       EXTERNAL           DLAMCH, DLARND
00467 *     ..
00468 *     .. External Subroutines ..
00469       EXTERNAL           ALASUM, DCOPY, DLABAD, DLAHD2, DSVDCH, XERBLA,
00470      $                   ZBDSQR, ZBDT01, ZBDT02, ZBDT03, ZGEBRD, ZGEMM,
00471      $                   ZLACPY, ZLASET, ZLATMR, ZLATMS, ZUNGBR, ZUNT01
00472 *     ..
00473 *     .. Intrinsic Functions ..
00474       INTRINSIC          ABS, EXP, INT, LOG, MAX, MIN, SQRT
00475 *     ..
00476 *     .. Scalars in Common ..
00477       LOGICAL            LERR, OK
00478       CHARACTER*32       SRNAMT
00479       INTEGER            INFOT, NUNIT
00480 *     ..
00481 *     .. Common blocks ..
00482       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
00483       COMMON             / SRNAMC / SRNAMT
00484 *     ..
00485 *     .. Data statements ..
00486       DATA               KTYPE / 1, 2, 5*4, 5*6, 3*9, 10 /
00487       DATA               KMAGN / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
00488       DATA               KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
00489      $                   0, 0, 0 /
00490 *     ..
00491 *     .. Executable Statements ..
00492 *
00493 *     Check for errors
00494 *
00495       INFO = 0
00496 *
00497       BADMM = .FALSE.
00498       BADNN = .FALSE.
00499       MMAX = 1
00500       NMAX = 1
00501       MNMAX = 1
00502       MINWRK = 1
00503       DO 10 J = 1, NSIZES
00504          MMAX = MAX( MMAX, MVAL( J ) )
00505          IF( MVAL( J ).LT.0 )
00506      $      BADMM = .TRUE.
00507          NMAX = MAX( NMAX, NVAL( J ) )
00508          IF( NVAL( J ).LT.0 )
00509      $      BADNN = .TRUE.
00510          MNMAX = MAX( MNMAX, MIN( MVAL( J ), NVAL( J ) ) )
00511          MINWRK = MAX( MINWRK, 3*( MVAL( J )+NVAL( J ) ),
00512      $            MVAL( J )*( MVAL( J )+MAX( MVAL( J ), NVAL( J ),
00513      $            NRHS )+1 )+NVAL( J )*MIN( NVAL( J ), MVAL( J ) ) )
00514    10 CONTINUE
00515 *
00516 *     Check for errors
00517 *
00518       IF( NSIZES.LT.0 ) THEN
00519          INFO = -1
00520       ELSE IF( BADMM ) THEN
00521          INFO = -2
00522       ELSE IF( BADNN ) THEN
00523          INFO = -3
00524       ELSE IF( NTYPES.LT.0 ) THEN
00525          INFO = -4
00526       ELSE IF( NRHS.LT.0 ) THEN
00527          INFO = -6
00528       ELSE IF( LDA.LT.MMAX ) THEN
00529          INFO = -11
00530       ELSE IF( LDX.LT.MMAX ) THEN
00531          INFO = -17
00532       ELSE IF( LDQ.LT.MMAX ) THEN
00533          INFO = -21
00534       ELSE IF( LDPT.LT.MNMAX ) THEN
00535          INFO = -23
00536       ELSE IF( MINWRK.GT.LWORK ) THEN
00537          INFO = -27
00538       END IF
00539 *
00540       IF( INFO.NE.0 ) THEN
00541          CALL XERBLA( 'ZCHKBD', -INFO )
00542          RETURN
00543       END IF
00544 *
00545 *     Initialize constants
00546 *
00547       PATH( 1: 1 ) = 'Zomplex precision'
00548       PATH( 2: 3 ) = 'BD'
00549       NFAIL = 0
00550       NTEST = 0
00551       UNFL = DLAMCH( 'Safe minimum' )
00552       OVFL = DLAMCH( 'Overflow' )
00553       CALL DLABAD( UNFL, OVFL )
00554       ULP = DLAMCH( 'Precision' )
00555       ULPINV = ONE / ULP
00556       LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) )
00557       RTUNFL = SQRT( UNFL )
00558       RTOVFL = SQRT( OVFL )
00559       INFOT = 0
00560 *
00561 *     Loop over sizes, types
00562 *
00563       DO 180 JSIZE = 1, NSIZES
00564          M = MVAL( JSIZE )
00565          N = NVAL( JSIZE )
00566          MNMIN = MIN( M, N )
00567          AMNINV = ONE / MAX( M, N, 1 )
00568 *
00569          IF( NSIZES.NE.1 ) THEN
00570             MTYPES = MIN( MAXTYP, NTYPES )
00571          ELSE
00572             MTYPES = MIN( MAXTYP+1, NTYPES )
00573          END IF
00574 *
00575          DO 170 JTYPE = 1, MTYPES
00576             IF( .NOT.DOTYPE( JTYPE ) )
00577      $         GO TO 170
00578 *
00579             DO 20 J = 1, 4
00580                IOLDSD( J ) = ISEED( J )
00581    20       CONTINUE
00582 *
00583             DO 30 J = 1, 14
00584                RESULT( J ) = -ONE
00585    30       CONTINUE
00586 *
00587             UPLO = ' '
00588 *
00589 *           Compute "A"
00590 *
00591 *           Control parameters:
00592 *
00593 *           KMAGN  KMODE        KTYPE
00594 *       =1  O(1)   clustered 1  zero
00595 *       =2  large  clustered 2  identity
00596 *       =3  small  exponential  (none)
00597 *       =4         arithmetic   diagonal, (w/ eigenvalues)
00598 *       =5         random       symmetric, w/ eigenvalues
00599 *       =6                      nonsymmetric, w/ singular values
00600 *       =7                      random diagonal
00601 *       =8                      random symmetric
00602 *       =9                      random nonsymmetric
00603 *       =10                     random bidiagonal (log. distrib.)
00604 *
00605             IF( MTYPES.GT.MAXTYP )
00606      $         GO TO 100
00607 *
00608             ITYPE = KTYPE( JTYPE )
00609             IMODE = KMODE( JTYPE )
00610 *
00611 *           Compute norm
00612 *
00613             GO TO ( 40, 50, 60 )KMAGN( JTYPE )
00614 *
00615    40       CONTINUE
00616             ANORM = ONE
00617             GO TO 70
00618 *
00619    50       CONTINUE
00620             ANORM = ( RTOVFL*ULP )*AMNINV
00621             GO TO 70
00622 *
00623    60       CONTINUE
00624             ANORM = RTUNFL*MAX( M, N )*ULPINV
00625             GO TO 70
00626 *
00627    70       CONTINUE
00628 *
00629             CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
00630             IINFO = 0
00631             COND = ULPINV
00632 *
00633             BIDIAG = .FALSE.
00634             IF( ITYPE.EQ.1 ) THEN
00635 *
00636 *              Zero matrix
00637 *
00638                IINFO = 0
00639 *
00640             ELSE IF( ITYPE.EQ.2 ) THEN
00641 *
00642 *              Identity
00643 *
00644                DO 80 JCOL = 1, MNMIN
00645                   A( JCOL, JCOL ) = ANORM
00646    80          CONTINUE
00647 *
00648             ELSE IF( ITYPE.EQ.4 ) THEN
00649 *
00650 *              Diagonal Matrix, [Eigen]values Specified
00651 *
00652                CALL ZLATMS( MNMIN, MNMIN, 'S', ISEED, 'N', RWORK, IMODE,
00653      $                      COND, ANORM, 0, 0, 'N', A, LDA, WORK,
00654      $                      IINFO )
00655 *
00656             ELSE IF( ITYPE.EQ.5 ) THEN
00657 *
00658 *              Symmetric, eigenvalues specified
00659 *
00660                CALL ZLATMS( MNMIN, MNMIN, 'S', ISEED, 'S', RWORK, IMODE,
00661      $                      COND, ANORM, M, N, 'N', A, LDA, WORK,
00662      $                      IINFO )
00663 *
00664             ELSE IF( ITYPE.EQ.6 ) THEN
00665 *
00666 *              Nonsymmetric, singular values specified
00667 *
00668                CALL ZLATMS( M, N, 'S', ISEED, 'N', RWORK, IMODE, COND,
00669      $                      ANORM, M, N, 'N', A, LDA, WORK, IINFO )
00670 *
00671             ELSE IF( ITYPE.EQ.7 ) THEN
00672 *
00673 *              Diagonal, random entries
00674 *
00675                CALL ZLATMR( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, 6, ONE,
00676      $                      CONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
00677      $                      WORK( 2*MNMIN+1 ), 1, ONE, 'N', IWORK, 0, 0,
00678      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00679 *
00680             ELSE IF( ITYPE.EQ.8 ) THEN
00681 *
00682 *              Symmetric, random entries
00683 *
00684                CALL ZLATMR( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, 6, ONE,
00685      $                      CONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
00686      $                      WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N,
00687      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00688 *
00689             ELSE IF( ITYPE.EQ.9 ) THEN
00690 *
00691 *              Nonsymmetric, random entries
00692 *
00693                CALL ZLATMR( M, N, 'S', ISEED, 'N', WORK, 6, ONE, CONE,
00694      $                      'T', 'N', WORK( MNMIN+1 ), 1, ONE,
00695      $                      WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N,
00696      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00697 *
00698             ELSE IF( ITYPE.EQ.10 ) THEN
00699 *
00700 *              Bidiagonal, random entries
00701 *
00702                TEMP1 = -TWO*LOG( ULP )
00703                DO 90 J = 1, MNMIN
00704                   BD( J ) = EXP( TEMP1*DLARND( 2, ISEED ) )
00705                   IF( J.LT.MNMIN )
00706      $               BE( J ) = EXP( TEMP1*DLARND( 2, ISEED ) )
00707    90          CONTINUE
00708 *
00709                IINFO = 0
00710                BIDIAG = .TRUE.
00711                IF( M.GE.N ) THEN
00712                   UPLO = 'U'
00713                ELSE
00714                   UPLO = 'L'
00715                END IF
00716             ELSE
00717                IINFO = 1
00718             END IF
00719 *
00720             IF( IINFO.EQ.0 ) THEN
00721 *
00722 *              Generate Right-Hand Side
00723 *
00724                IF( BIDIAG ) THEN
00725                   CALL ZLATMR( MNMIN, NRHS, 'S', ISEED, 'N', WORK, 6,
00726      $                         ONE, CONE, 'T', 'N', WORK( MNMIN+1 ), 1,
00727      $                         ONE, WORK( 2*MNMIN+1 ), 1, ONE, 'N',
00728      $                         IWORK, MNMIN, NRHS, ZERO, ONE, 'NO', Y,
00729      $                         LDX, IWORK, IINFO )
00730                ELSE
00731                   CALL ZLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE,
00732      $                         CONE, 'T', 'N', WORK( M+1 ), 1, ONE,
00733      $                         WORK( 2*M+1 ), 1, ONE, 'N', IWORK, M,
00734      $                         NRHS, ZERO, ONE, 'NO', X, LDX, IWORK,
00735      $                         IINFO )
00736                END IF
00737             END IF
00738 *
00739 *           Error Exit
00740 *
00741             IF( IINFO.NE.0 ) THEN
00742                WRITE( NOUT, FMT = 9998 )'Generator', IINFO, M, N,
00743      $            JTYPE, IOLDSD
00744                INFO = ABS( IINFO )
00745                RETURN
00746             END IF
00747 *
00748   100       CONTINUE
00749 *
00750 *           Call ZGEBRD and ZUNGBR to compute B, Q, and P, do tests.
00751 *
00752             IF( .NOT.BIDIAG ) THEN
00753 *
00754 *              Compute transformations to reduce A to bidiagonal form:
00755 *              B := Q' * A * P.
00756 *
00757                CALL ZLACPY( ' ', M, N, A, LDA, Q, LDQ )
00758                CALL ZGEBRD( M, N, Q, LDQ, BD, BE, WORK, WORK( MNMIN+1 ),
00759      $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
00760 *
00761 *              Check error code from ZGEBRD.
00762 *
00763                IF( IINFO.NE.0 ) THEN
00764                   WRITE( NOUT, FMT = 9998 )'ZGEBRD', IINFO, M, N,
00765      $               JTYPE, IOLDSD
00766                   INFO = ABS( IINFO )
00767                   RETURN
00768                END IF
00769 *
00770                CALL ZLACPY( ' ', M, N, Q, LDQ, PT, LDPT )
00771                IF( M.GE.N ) THEN
00772                   UPLO = 'U'
00773                ELSE
00774                   UPLO = 'L'
00775                END IF
00776 *
00777 *              Generate Q
00778 *
00779                MQ = M
00780                IF( NRHS.LE.0 )
00781      $            MQ = MNMIN
00782                CALL ZUNGBR( 'Q', M, MQ, N, Q, LDQ, WORK,
00783      $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
00784 *
00785 *              Check error code from ZUNGBR.
00786 *
00787                IF( IINFO.NE.0 ) THEN
00788                   WRITE( NOUT, FMT = 9998 )'ZUNGBR(Q)', IINFO, M, N,
00789      $               JTYPE, IOLDSD
00790                   INFO = ABS( IINFO )
00791                   RETURN
00792                END IF
00793 *
00794 *              Generate P'
00795 *
00796                CALL ZUNGBR( 'P', MNMIN, N, M, PT, LDPT, WORK( MNMIN+1 ),
00797      $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
00798 *
00799 *              Check error code from ZUNGBR.
00800 *
00801                IF( IINFO.NE.0 ) THEN
00802                   WRITE( NOUT, FMT = 9998 )'ZUNGBR(P)', IINFO, M, N,
00803      $               JTYPE, IOLDSD
00804                   INFO = ABS( IINFO )
00805                   RETURN
00806                END IF
00807 *
00808 *              Apply Q' to an M by NRHS matrix X:  Y := Q' * X.
00809 *
00810                CALL ZGEMM( 'Conjugate transpose', 'No transpose', M,
00811      $                     NRHS, M, CONE, Q, LDQ, X, LDX, CZERO, Y,
00812      $                     LDX )
00813 *
00814 *              Test 1:  Check the decomposition A := Q * B * PT
00815 *                   2:  Check the orthogonality of Q
00816 *                   3:  Check the orthogonality of PT
00817 *
00818                CALL ZBDT01( M, N, 1, A, LDA, Q, LDQ, BD, BE, PT, LDPT,
00819      $                      WORK, RWORK, RESULT( 1 ) )
00820                CALL ZUNT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
00821      $                      RWORK, RESULT( 2 ) )
00822                CALL ZUNT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
00823      $                      RWORK, RESULT( 3 ) )
00824             END IF
00825 *
00826 *           Use ZBDSQR to form the SVD of the bidiagonal matrix B:
00827 *           B := U * S1 * VT, and compute Z = U' * Y.
00828 *
00829             CALL DCOPY( MNMIN, BD, 1, S1, 1 )
00830             IF( MNMIN.GT.0 )
00831      $         CALL DCOPY( MNMIN-1, BE, 1, RWORK, 1 )
00832             CALL ZLACPY( ' ', M, NRHS, Y, LDX, Z, LDX )
00833             CALL ZLASET( 'Full', MNMIN, MNMIN, CZERO, CONE, U, LDPT )
00834             CALL ZLASET( 'Full', MNMIN, MNMIN, CZERO, CONE, VT, LDPT )
00835 *
00836             CALL ZBDSQR( UPLO, MNMIN, MNMIN, MNMIN, NRHS, S1, RWORK, VT,
00837      $                   LDPT, U, LDPT, Z, LDX, RWORK( MNMIN+1 ),
00838      $                   IINFO )
00839 *
00840 *           Check error code from ZBDSQR.
00841 *
00842             IF( IINFO.NE.0 ) THEN
00843                WRITE( NOUT, FMT = 9998 )'ZBDSQR(vects)', IINFO, M, N,
00844      $            JTYPE, IOLDSD
00845                INFO = ABS( IINFO )
00846                IF( IINFO.LT.0 ) THEN
00847                   RETURN
00848                ELSE
00849                   RESULT( 4 ) = ULPINV
00850                   GO TO 150
00851                END IF
00852             END IF
00853 *
00854 *           Use ZBDSQR to compute only the singular values of the
00855 *           bidiagonal matrix B;  U, VT, and Z should not be modified.
00856 *
00857             CALL DCOPY( MNMIN, BD, 1, S2, 1 )
00858             IF( MNMIN.GT.0 )
00859      $         CALL DCOPY( MNMIN-1, BE, 1, RWORK, 1 )
00860 *
00861             CALL ZBDSQR( UPLO, MNMIN, 0, 0, 0, S2, RWORK, VT, LDPT, U,
00862      $                   LDPT, Z, LDX, RWORK( MNMIN+1 ), IINFO )
00863 *
00864 *           Check error code from ZBDSQR.
00865 *
00866             IF( IINFO.NE.0 ) THEN
00867                WRITE( NOUT, FMT = 9998 )'ZBDSQR(values)', IINFO, M, N,
00868      $            JTYPE, IOLDSD
00869                INFO = ABS( IINFO )
00870                IF( IINFO.LT.0 ) THEN
00871                   RETURN
00872                ELSE
00873                   RESULT( 9 ) = ULPINV
00874                   GO TO 150
00875                END IF
00876             END IF
00877 *
00878 *           Test 4:  Check the decomposition B := U * S1 * VT
00879 *                5:  Check the computation Z := U' * Y
00880 *                6:  Check the orthogonality of U
00881 *                7:  Check the orthogonality of VT
00882 *
00883             CALL ZBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT,
00884      $                   WORK, RESULT( 4 ) )
00885             CALL ZBDT02( MNMIN, NRHS, Y, LDX, Z, LDX, U, LDPT, WORK,
00886      $                   RWORK, RESULT( 5 ) )
00887             CALL ZUNT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK,
00888      $                   RWORK, RESULT( 6 ) )
00889             CALL ZUNT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK,
00890      $                   RWORK, RESULT( 7 ) )
00891 *
00892 *           Test 8:  Check that the singular values are sorted in
00893 *                    non-increasing order and are non-negative
00894 *
00895             RESULT( 8 ) = ZERO
00896             DO 110 I = 1, MNMIN - 1
00897                IF( S1( I ).LT.S1( I+1 ) )
00898      $            RESULT( 8 ) = ULPINV
00899                IF( S1( I ).LT.ZERO )
00900      $            RESULT( 8 ) = ULPINV
00901   110       CONTINUE
00902             IF( MNMIN.GE.1 ) THEN
00903                IF( S1( MNMIN ).LT.ZERO )
00904      $            RESULT( 8 ) = ULPINV
00905             END IF
00906 *
00907 *           Test 9:  Compare ZBDSQR with and without singular vectors
00908 *
00909             TEMP2 = ZERO
00910 *
00911             DO 120 J = 1, MNMIN
00912                TEMP1 = ABS( S1( J )-S2( J ) ) /
00913      $                 MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
00914      $                 ULP*MAX( ABS( S1( J ) ), ABS( S2( J ) ) ) )
00915                TEMP2 = MAX( TEMP1, TEMP2 )
00916   120       CONTINUE
00917 *
00918             RESULT( 9 ) = TEMP2
00919 *
00920 *           Test 10:  Sturm sequence test of singular values
00921 *                     Go up by factors of two until it succeeds
00922 *
00923             TEMP1 = THRESH*( HALF-ULP )
00924 *
00925             DO 130 J = 0, LOG2UI
00926                CALL DSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO )
00927                IF( IINFO.EQ.0 )
00928      $            GO TO 140
00929                TEMP1 = TEMP1*TWO
00930   130       CONTINUE
00931 *
00932   140       CONTINUE
00933             RESULT( 10 ) = TEMP1
00934 *
00935 *           Use ZBDSQR to form the decomposition A := (QU) S (VT PT)
00936 *           from the bidiagonal form A := Q B PT.
00937 *
00938             IF( .NOT.BIDIAG ) THEN
00939                CALL DCOPY( MNMIN, BD, 1, S2, 1 )
00940                IF( MNMIN.GT.0 )
00941      $            CALL DCOPY( MNMIN-1, BE, 1, RWORK, 1 )
00942 *
00943                CALL ZBDSQR( UPLO, MNMIN, N, M, NRHS, S2, RWORK, PT,
00944      $                      LDPT, Q, LDQ, Y, LDX, RWORK( MNMIN+1 ),
00945      $                      IINFO )
00946 *
00947 *              Test 11:  Check the decomposition A := Q*U * S2 * VT*PT
00948 *                   12:  Check the computation Z := U' * Q' * X
00949 *                   13:  Check the orthogonality of Q*U
00950 *                   14:  Check the orthogonality of VT*PT
00951 *
00952                CALL ZBDT01( M, N, 0, A, LDA, Q, LDQ, S2, DUMMA, PT,
00953      $                      LDPT, WORK, RWORK, RESULT( 11 ) )
00954                CALL ZBDT02( M, NRHS, X, LDX, Y, LDX, Q, LDQ, WORK,
00955      $                      RWORK, RESULT( 12 ) )
00956                CALL ZUNT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
00957      $                      RWORK, RESULT( 13 ) )
00958                CALL ZUNT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
00959      $                      RWORK, RESULT( 14 ) )
00960             END IF
00961 *
00962 *           End of Loop -- Check for RESULT(j) > THRESH
00963 *
00964   150       CONTINUE
00965             DO 160 J = 1, 14
00966                IF( RESULT( J ).GE.THRESH ) THEN
00967                   IF( NFAIL.EQ.0 )
00968      $               CALL DLAHD2( NOUT, PATH )
00969                   WRITE( NOUT, FMT = 9999 )M, N, JTYPE, IOLDSD, J,
00970      $               RESULT( J )
00971                   NFAIL = NFAIL + 1
00972                END IF
00973   160       CONTINUE
00974             IF( .NOT.BIDIAG ) THEN
00975                NTEST = NTEST + 14
00976             ELSE
00977                NTEST = NTEST + 5
00978             END IF
00979 *
00980   170    CONTINUE
00981   180 CONTINUE
00982 *
00983 *     Summary
00984 *
00985       CALL ALASUM( PATH, NOUT, NFAIL, NTEST, 0 )
00986 *
00987       RETURN
00988 *
00989 *     End of ZCHKBD
00990 *
00991  9999 FORMAT( ' M=', I5, ', N=', I5, ', type ', I2, ', seed=',
00992      $      4( I4, ',' ), ' test(', I2, ')=', G11.4 )
00993  9998 FORMAT( ' ZCHKBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
00994      $      I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
00995      $      I5, ')' )
00996 *
00997       END
 All Files Functions