TESTING/EIG/dchkbd.f

   1 *> \brief \b DCHKBD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE DCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
  12 *                          ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX,
  13 *                          Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK,
  14 *                          IWORK, NOUT, INFO )
  15 *
  16 *       .. Scalar Arguments ..
  17 *       INTEGER            INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
  18 *      $                   NSIZES, NTYPES
  19 *       DOUBLE PRECISION   THRESH
  20 *       ..
  21 *       .. Array Arguments ..
  22 *       LOGICAL            DOTYPE( * )
  23 *       INTEGER            ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * )
  24 *       DOUBLE PRECISION   A( LDA, * ), BD( * ), BE( * ), PT( LDPT, * ),
  25 *      $                   Q( LDQ, * ), S1( * ), S2( * ), U( LDPT, * ),
  26 *      $                   VT( LDPT, * ), WORK( * ), X( LDX, * ),
  27 *      $                   Y( LDX, * ), Z( LDX, * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> DCHKBD checks the singular value decomposition (SVD) routines.
  37 *>
  38 *> DGEBRD reduces a real general m by n matrix A to upper or lower
  39 *> bidiagonal form B by an orthogonal transformation:  Q' * A * P = B
  40 *> (or A = Q * B * P').  The matrix B is upper bidiagonal if m >= n
  41 *> and lower bidiagonal if m < n.
  42 *>
  43 *> DORGBR generates the orthogonal matrices Q and P' from DGEBRD.
  44 *> Note that Q and P are not necessarily square.
  45 *>
  46 *> DBDSQR computes the singular value decomposition of the bidiagonal
  47 *> matrix B as B = U S V'.  It is called three times to compute
  48 *>    1)  B = U S1 V', where S1 is the diagonal matrix of singular
  49 *>        values and the columns of the matrices U and V are the left
  50 *>        and right singular vectors, respectively, of B.
  51 *>    2)  Same as 1), but the singular values are stored in S2 and the
  52 *>        singular vectors are not computed.
  53 *>    3)  A = (UQ) S (P'V'), the SVD of the original matrix A.
  54 *> In addition, DBDSQR has an option to apply the left orthogonal matrix
  55 *> U to a matrix X, useful in least squares applications.
  56 *>
  57 *> DBDSDC computes the singular value decomposition of the bidiagonal
  58 *> matrix B as B = U S V' using divide-and-conquer. It is called twice
  59 *> to compute
  60 *>    1) B = U S1 V', where S1 is the diagonal matrix of singular
  61 *>        values and the columns of the matrices U and V are the left
  62 *>        and right singular vectors, respectively, of B.
  63 *>    2) Same as 1), but the singular values are stored in S2 and the
  64 *>        singular vectors are not computed.
  65 *>
  66 *>  DBDSVDX computes the singular value decomposition of the bidiagonal
  67 *>  matrix B as B = U S V' using bisection and inverse iteration. It is
  68 *>  called six times to compute
  69 *>     1) B = U S1 V', RANGE='A', where S1 is the diagonal matrix of singular
  70 *>         values and the columns of the matrices U and V are the left
  71 *>         and right singular vectors, respectively, of B.
  72 *>     2) Same as 1), but the singular values are stored in S2 and the
  73 *>         singular vectors are not computed.
  74 *>     3) B = U S1 V', RANGE='I', with where S1 is the diagonal matrix of singular
  75 *>         values and the columns of the matrices U and V are the left
  76 *>         and right singular vectors, respectively, of B
  77 *>     4) Same as 3), but the singular values are stored in S2 and the
  78 *>         singular vectors are not computed.
  79 *>     5) B = U S1 V', RANGE='V', with where S1 is the diagonal matrix of singular
  80 *>         values and the columns of the matrices U and V are the left
  81 *>         and right singular vectors, respectively, of B
  82 *>     6) Same as 5), but the singular values are stored in S2 and the
  83 *>         singular vectors are not computed.
  84 *>
  85 *> For each pair of matrix dimensions (M,N) and each selected matrix
  86 *> type, an M by N matrix A and an M by NRHS matrix X are generated.
  87 *> The problem dimensions are as follows
  88 *>    A:          M x N
  89 *>    Q:          M x min(M,N) (but M x M if NRHS > 0)
  90 *>    P:          min(M,N) x N
  91 *>    B:          min(M,N) x min(M,N)
  92 *>    U, V:       min(M,N) x min(M,N)
  93 *>    S1, S2      diagonal, order min(M,N)
  94 *>    X:          M x NRHS
  95 *>
  96 *> For each generated matrix, 14 tests are performed:
  97 *>
  98 *> Test DGEBRD and DORGBR
  99 *>
 100 *> (1)   | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P'
 101 *>
 102 *> (2)   | I - Q' Q | / ( M ulp )
 103 *>
 104 *> (3)   | I - PT PT' | / ( N ulp )
 105 *>
 106 *> Test DBDSQR on bidiagonal matrix B
 107 *>
 108 *> (4)   | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
 109 *>
 110 *> (5)   | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X
 111 *>                                                  and   Z = U' Y.
 112 *> (6)   | I - U' U | / ( min(M,N) ulp )
 113 *>
 114 *> (7)   | I - VT VT' | / ( min(M,N) ulp )
 115 *>
 116 *> (8)   S1 contains min(M,N) nonnegative values in decreasing order.
 117 *>       (Return 0 if true, 1/ULP if false.)
 118 *>
 119 *> (9)   | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
 120 *>                                   computing U and V.
 121 *>
 122 *> (10)  0 if the true singular values of B are within THRESH of
 123 *>       those in S1.  2*THRESH if they are not.  (Tested using
 124 *>       DSVDCH)
 125 *>
 126 *> Test DBDSQR on matrix A
 127 *>
 128 *> (11)  | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp )
 129 *>
 130 *> (12)  | X - (QU) Z | / ( |X| max(M,k) ulp )
 131 *>
 132 *> (13)  | I - (QU)'(QU) | / ( M ulp )
 133 *>
 134 *> (14)  | I - (VT PT) (PT'VT') | / ( N ulp )
 135 *>
 136 *> Test DBDSDC on bidiagonal matrix B
 137 *>
 138 *> (15)  | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
 139 *>
 140 *> (16)  | I - U' U | / ( min(M,N) ulp )
 141 *>
 142 *> (17)  | I - VT VT' | / ( min(M,N) ulp )
 143 *>
 144 *> (18)  S1 contains min(M,N) nonnegative values in decreasing order.
 145 *>       (Return 0 if true, 1/ULP if false.)
 146 *>
 147 *> (19)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
 148 *>                                   computing U and V.
 149 *>  Test DBDSVDX on bidiagonal matrix B
 150 *>
 151 *>  (20)  | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
 152 *>
 153 *>  (21)  | I - U' U | / ( min(M,N) ulp )
 154 *>
 155 *>  (22)  | I - VT VT' | / ( min(M,N) ulp )
 156 *>
 157 *>  (23)  S1 contains min(M,N) nonnegative values in decreasing order.
 158 *>        (Return 0 if true, 1/ULP if false.)
 159 *>
 160 *>  (24)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
 161 *>                                    computing U and V.
 162 *>
 163 *>  (25)  | S1 - U' B VT' | / ( |S| n ulp )    DBDSVDX('V', 'I')
 164 *>
 165 *>  (26)  | I - U' U | / ( min(M,N) ulp )
 166 *>
 167 *>  (27)  | I - VT VT' | / ( min(M,N) ulp )
 168 *>
 169 *>  (28)  S1 contains min(M,N) nonnegative values in decreasing order.
 170 *>        (Return 0 if true, 1/ULP if false.)
 171 *>
 172 *>  (29)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
 173 *>                                    computing U and V.
 174 *>
 175 *>  (30)  | S1 - U' B VT' | / ( |S1| n ulp )   DBDSVDX('V', 'V')
 176 *>
 177 *>  (31)  | I - U' U | / ( min(M,N) ulp )
 178 *>
 179 *>  (32)  | I - VT VT' | / ( min(M,N) ulp )
 180 *>
 181 *>  (33)  S1 contains min(M,N) nonnegative values in decreasing order.
 182 *>        (Return 0 if true, 1/ULP if false.)
 183 *>
 184 *>  (34)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
 185 *>                                    computing U and V.
 186 *>
 187 *> The possible matrix types are
 188 *>
 189 *> (1)  The zero matrix.
 190 *> (2)  The identity matrix.
 191 *>
 192 *> (3)  A diagonal matrix with evenly spaced entries
 193 *>      1, ..., ULP  and random signs.
 194 *>      (ULP = (first number larger than 1) - 1 )
 195 *> (4)  A diagonal matrix with geometrically spaced entries
 196 *>      1, ..., ULP  and random signs.
 197 *> (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
 198 *>      and random signs.
 199 *>
 200 *> (6)  Same as (3), but multiplied by SQRT( overflow threshold )
 201 *> (7)  Same as (3), but multiplied by SQRT( underflow threshold )
 202 *>
 203 *> (8)  A matrix of the form  U D V, where U and V are orthogonal and
 204 *>      D has evenly spaced entries 1, ..., ULP with random signs
 205 *>      on the diagonal.
 206 *>
 207 *> (9)  A matrix of the form  U D V, where U and V are orthogonal and
 208 *>      D has geometrically spaced entries 1, ..., ULP with random
 209 *>      signs on the diagonal.
 210 *>
 211 *> (10) A matrix of the form  U D V, where U and V are orthogonal and
 212 *>      D has "clustered" entries 1, ULP,..., ULP with random
 213 *>      signs on the diagonal.
 214 *>
 215 *> (11) Same as (8), but multiplied by SQRT( overflow threshold )
 216 *> (12) Same as (8), but multiplied by SQRT( underflow threshold )
 217 *>
 218 *> (13) Rectangular matrix with random entries chosen from (-1,1).
 219 *> (14) Same as (13), but multiplied by SQRT( overflow threshold )
 220 *> (15) Same as (13), but multiplied by SQRT( underflow threshold )
 221 *>
 222 *> Special case:
 223 *> (16) A bidiagonal matrix with random entries chosen from a
 224 *>      logarithmic distribution on [ulp^2,ulp^(-2)]  (I.e., each
 225 *>      entry is  e^x, where x is chosen uniformly on
 226 *>      [ 2 log(ulp), -2 log(ulp) ] .)  For *this* type:
 227 *>      (a) DGEBRD is not called to reduce it to bidiagonal form.
 228 *>      (b) the bidiagonal is  min(M,N) x min(M,N); if M<N, the
 229 *>          matrix will be lower bidiagonal, otherwise upper.
 230 *>      (c) only tests 5--8 and 14 are performed.
 231 *>
 232 *> A subset of the full set of matrix types may be selected through
 233 *> the logical array DOTYPE.
 234 *> \endverbatim
 235 *
 236 *  Arguments:
 237 *  ==========
 238 *
 239 *> \param[in] NSIZES
 240 *> \verbatim
 241 *>          NSIZES is INTEGER
 242 *>          The number of values of M and N contained in the vectors
 243 *>          MVAL and NVAL.  The matrix sizes are used in pairs (M,N).
 244 *> \endverbatim
 245 *>
 246 *> \param[in] MVAL
 247 *> \verbatim
 248 *>          MVAL is INTEGER array, dimension (NM)
 249 *>          The values of the matrix row dimension M.
 250 *> \endverbatim
 251 *>
 252 *> \param[in] NVAL
 253 *> \verbatim
 254 *>          NVAL is INTEGER array, dimension (NM)
 255 *>          The values of the matrix column dimension N.
 256 *> \endverbatim
 257 *>
 258 *> \param[in] NTYPES
 259 *> \verbatim
 260 *>          NTYPES is INTEGER
 261 *>          The number of elements in DOTYPE.   If it is zero, DCHKBD
 262 *>          does nothing.  It must be at least zero.  If it is MAXTYP+1
 263 *>          and NSIZES is 1, then an additional type, MAXTYP+1 is
 264 *>          defined, which is to use whatever matrices are in A and B.
 265 *>          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
 266 *>          DOTYPE(MAXTYP+1) is .TRUE. .
 267 *> \endverbatim
 268 *>
 269 *> \param[in] DOTYPE
 270 *> \verbatim
 271 *>          DOTYPE is LOGICAL array, dimension (NTYPES)
 272 *>          If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
 273 *>          of type j will be generated.  If NTYPES is smaller than the
 274 *>          maximum number of types defined (PARAMETER MAXTYP), then
 275 *>          types NTYPES+1 through MAXTYP will not be generated.  If
 276 *>          NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
 277 *>          DOTYPE(NTYPES) will be ignored.
 278 *> \endverbatim
 279 *>
 280 *> \param[in] NRHS
 281 *> \verbatim
 282 *>          NRHS is INTEGER
 283 *>          The number of columns in the "right-hand side" matrices X, Y,
 284 *>          and Z, used in testing DBDSQR.  If NRHS = 0, then the
 285 *>          operations on the right-hand side will not be tested.
 286 *>          NRHS must be at least 0.
 287 *> \endverbatim
 288 *>
 289 *> \param[in,out] ISEED
 290 *> \verbatim
 291 *>          ISEED is INTEGER array, dimension (4)
 292 *>          On entry ISEED specifies the seed of the random number
 293 *>          generator. The array elements should be between 0 and 4095;
 294 *>          if not they will be reduced mod 4096.  Also, ISEED(4) must
 295 *>          be odd.  The values of ISEED are changed on exit, and can be
 296 *>          used in the next call to DCHKBD to continue the same random
 297 *>          number sequence.
 298 *> \endverbatim
 299 *>
 300 *> \param[in] THRESH
 301 *> \verbatim
 302 *>          THRESH is DOUBLE PRECISION
 303 *>          The threshold value for the test ratios.  A result is
 304 *>          included in the output file if RESULT >= THRESH.  To have
 305 *>          every test ratio printed, use THRESH = 0.  Note that the
 306 *>          expected value of the test ratios is O(1), so THRESH should
 307 *>          be a reasonably small multiple of 1, e.g., 10 or 100.
 308 *> \endverbatim
 309 *>
 310 *> \param[out] A
 311 *> \verbatim
 312 *>          A is DOUBLE PRECISION array, dimension (LDA,NMAX)
 313 *>          where NMAX is the maximum value of N in NVAL.
 314 *> \endverbatim
 315 *>
 316 *> \param[in] LDA
 317 *> \verbatim
 318 *>          LDA is INTEGER
 319 *>          The leading dimension of the array A.  LDA >= max(1,MMAX),
 320 *>          where MMAX is the maximum value of M in MVAL.
 321 *> \endverbatim
 322 *>
 323 *> \param[out] BD
 324 *> \verbatim
 325 *>          BD is DOUBLE PRECISION array, dimension
 326 *>                      (max(min(MVAL(j),NVAL(j))))
 327 *> \endverbatim
 328 *>
 329 *> \param[out] BE
 330 *> \verbatim
 331 *>          BE is DOUBLE PRECISION array, dimension
 332 *>                      (max(min(MVAL(j),NVAL(j))))
 333 *> \endverbatim
 334 *>
 335 *> \param[out] S1
 336 *> \verbatim
 337 *>          S1 is DOUBLE PRECISION array, dimension
 338 *>                      (max(min(MVAL(j),NVAL(j))))
 339 *> \endverbatim
 340 *>
 341 *> \param[out] S2
 342 *> \verbatim
 343 *>          S2 is DOUBLE PRECISION array, dimension
 344 *>                      (max(min(MVAL(j),NVAL(j))))
 345 *> \endverbatim
 346 *>
 347 *> \param[out] X
 348 *> \verbatim
 349 *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
 350 *> \endverbatim
 351 *>
 352 *> \param[in] LDX
 353 *> \verbatim
 354 *>          LDX is INTEGER
 355 *>          The leading dimension of the arrays X, Y, and Z.
 356 *>          LDX >= max(1,MMAX)
 357 *> \endverbatim
 358 *>
 359 *> \param[out] Y
 360 *> \verbatim
 361 *>          Y is DOUBLE PRECISION array, dimension (LDX,NRHS)
 362 *> \endverbatim
 363 *>
 364 *> \param[out] Z
 365 *> \verbatim
 366 *>          Z is DOUBLE PRECISION array, dimension (LDX,NRHS)
 367 *> \endverbatim
 368 *>
 369 *> \param[out] Q
 370 *> \verbatim
 371 *>          Q is DOUBLE PRECISION array, dimension (LDQ,MMAX)
 372 *> \endverbatim
 373 *>
 374 *> \param[in] LDQ
 375 *> \verbatim
 376 *>          LDQ is INTEGER
 377 *>          The leading dimension of the array Q.  LDQ >= max(1,MMAX).
 378 *> \endverbatim
 379 *>
 380 *> \param[out] PT
 381 *> \verbatim
 382 *>          PT is DOUBLE PRECISION array, dimension (LDPT,NMAX)
 383 *> \endverbatim
 384 *>
 385 *> \param[in] LDPT
 386 *> \verbatim
 387 *>          LDPT is INTEGER
 388 *>          The leading dimension of the arrays PT, U, and V.
 389 *>          LDPT >= max(1, max(min(MVAL(j),NVAL(j)))).
 390 *> \endverbatim
 391 *>
 392 *> \param[out] U
 393 *> \verbatim
 394 *>          U is DOUBLE PRECISION array, dimension
 395 *>                      (LDPT,max(min(MVAL(j),NVAL(j))))
 396 *> \endverbatim
 397 *>
 398 *> \param[out] VT
 399 *> \verbatim
 400 *>          VT is DOUBLE PRECISION array, dimension
 401 *>                      (LDPT,max(min(MVAL(j),NVAL(j))))
 402 *> \endverbatim
 403 *>
 404 *> \param[out] WORK
 405 *> \verbatim
 406 *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
 407 *> \endverbatim
 408 *>
 409 *> \param[in] LWORK
 410 *> \verbatim
 411 *>          LWORK is INTEGER
 412 *>          The number of entries in WORK.  This must be at least
 413 *>          3(M+N) and  M(M + max(M,N,k) + 1) + N*min(M,N)  for all
 414 *>          pairs  (M,N)=(MM(j),NN(j))
 415 *> \endverbatim
 416 *>
 417 *> \param[out] IWORK
 418 *> \verbatim
 419 *>          IWORK is INTEGER array, dimension at least 8*min(M,N)
 420 *> \endverbatim
 421 *>
 422 *> \param[in] NOUT
 423 *> \verbatim
 424 *>          NOUT is INTEGER
 425 *>          The FORTRAN unit number for printing out error messages
 426 *>          (e.g., if a routine returns IINFO not equal to 0.)
 427 *> \endverbatim
 428 *>
 429 *> \param[out] INFO
 430 *> \verbatim
 431 *>          INFO is INTEGER
 432 *>          If 0, then everything ran OK.
 433 *>           -1: NSIZES < 0
 434 *>           -2: Some MM(j) < 0
 435 *>           -3: Some NN(j) < 0
 436 *>           -4: NTYPES < 0
 437 *>           -6: NRHS  < 0
 438 *>           -8: THRESH < 0
 439 *>          -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
 440 *>          -17: LDB < 1 or LDB < MMAX.
 441 *>          -21: LDQ < 1 or LDQ < MMAX.
 442 *>          -23: LDPT< 1 or LDPT< MNMAX.
 443 *>          -27: LWORK too small.
 444 *>          If  DLATMR, SLATMS, DGEBRD, DORGBR, or DBDSQR,
 445 *>              returns an error code, the
 446 *>              absolute value of it is returned.
 447 *>
 448 *>-----------------------------------------------------------------------
 449 *>
 450 *>     Some Local Variables and Parameters:
 451 *>     ---- ----- --------- --- ----------
 452 *>
 453 *>     ZERO, ONE       Real 0 and 1.
 454 *>     MAXTYP          The number of types defined.
 455 *>     NTEST           The number of tests performed, or which can
 456 *>                     be performed so far, for the current matrix.
 457 *>     MMAX            Largest value in NN.
 458 *>     NMAX            Largest value in NN.
 459 *>     MNMIN           min(MM(j), NN(j)) (the dimension of the bidiagonal
 460 *>                     matrix.)
 461 *>     MNMAX           The maximum value of MNMIN for j=1,...,NSIZES.
 462 *>     NFAIL           The number of tests which have exceeded THRESH
 463 *>     COND, IMODE     Values to be passed to the matrix generators.
 464 *>     ANORM           Norm of A; passed to matrix generators.
 465 *>
 466 *>     OVFL, UNFL      Overflow and underflow thresholds.
 467 *>     RTOVFL, RTUNFL  Square roots of the previous 2 values.
 468 *>     ULP, ULPINV     Finest relative precision and its inverse.
 469 *>
 470 *>             The following four arrays decode JTYPE:
 471 *>     KTYPE(j)        The general type (1-10) for type "j".
 472 *>     KMODE(j)        The MODE value to be passed to the matrix
 473 *>                     generator for type "j".
 474 *>     KMAGN(j)        The order of magnitude ( O(1),
 475 *>                     O(overflow^(1/2) ), O(underflow^(1/2) )
 476 *> \endverbatim
 477 *
 478 *  Authors:
 479 *  ========
 480 *
 481 *> \author Univ. of Tennessee
 482 *> \author Univ. of California Berkeley
 483 *> \author Univ. of Colorado Denver
 484 *> \author NAG Ltd.
 485 *
 486 *> \date June 2016
 487 *
 488 *> \ingroup double_eig
 489 *
 490 *  =====================================================================
 491       SUBROUTINE DCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
 492      $                   ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX,
 493      $                   Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK,
 494      $                   IWORK, NOUT, INFO )
 495 *
 496 *  -- LAPACK test routine (version 3.6.1) --
 497 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 498 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 499 *     June 2016
 500 *
 501 *     .. Scalar Arguments ..
 502       INTEGER            INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
 503      $                   NSIZES, NTYPES
 504       DOUBLE PRECISION   THRESH
 505 *     ..
 506 *     .. Array Arguments ..
 507       LOGICAL            DOTYPE( * )
 508       INTEGER            ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * )
 509       DOUBLE PRECISION   A( LDA, * ), BD( * ), BE( * ), PT( LDPT, * ),
 510      $                   Q( LDQ, * ), S1( * ), S2( * ), U( LDPT, * ),
 511      $                   VT( LDPT, * ), WORK( * ), X( LDX, * ),
 512      $                   Y( LDX, * ), Z( LDX, * )
 513 *     ..
 514 *
 515 * ======================================================================
 516 *
 517 *     .. Parameters ..
 518       DOUBLE PRECISION   ZERO, ONE, TWO, HALF
 519       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
 520      $                   HALF = 0.5D0 )
 521       INTEGER            MAXTYP
 522       PARAMETER          ( MAXTYP = 16 )
 523 *     ..
 524 *     .. Local Scalars ..
 525       LOGICAL            BADMM, BADNN, BIDIAG
 526       CHARACTER          UPLO
 527       CHARACTER*3        PATH
 528       INTEGER            I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, IWBD,
 529      $                   IWBE, IWBS, IWBZ, IWWORK, J, JCOL, JSIZE,
 530      $                   JTYPE, LOG2UI, M, MINWRK, MMAX, MNMAX, MNMIN,
 531      $                   MNMIN2, MQ, MTYPES, N, NFAIL, NMAX,
 532      $                   NS1, NS2, NTEST
 533       DOUBLE PRECISION   ABSTOL, AMNINV, ANORM, COND, OVFL, RTOVFL,
 534      $                   RTUNFL, TEMP1, TEMP2, TEMP3, ULP, ULPINV,
 535      $                   UNFL, VL, VU
 536 *     ..
 537 *     .. Local Arrays ..
 538       INTEGER            IDUM( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
 539      $                   KMAGN( MAXTYP ), KMODE( MAXTYP ),
 540      $                   KTYPE( MAXTYP )
 541       DOUBLE PRECISION   DUM( 1 ), DUMMA( 1 ), RESULT( 40 )
 542 *     ..
 543 *     .. External Functions ..
 544       DOUBLE PRECISION   DLAMCH, DLARND, DSXT1
 545       EXTERNAL           DLAMCH, DLARND, DSXT1
 546 *     ..
 547 *     .. External Subroutines ..
 548       EXTERNAL           ALASUM, DBDSDC, DBDSQR, DBDSVDX, DBDT01,
 549      $                   DBDT02, DBDT03, DBDT04, DCOPY, DGEBRD,
 550      $                   DGEMM, DLABAD, DLACPY, DLAHD2, DLASET,
 551      $                   DLATMR, DLATMS, DORGBR, DORT01, XERBLA
 552 *     ..
 553 *     .. Intrinsic Functions ..
 554       INTRINSIC          ABS, EXP, INT, LOG, MAX, MIN, SQRT
 555 *     ..
 556 *     .. Scalars in Common ..
 557       LOGICAL            LERR, OK
 558       CHARACTER*32       SRNAMT
 559       INTEGER            INFOT, NUNIT
 560 *     ..
 561 *     .. Common blocks ..
 562       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
 563       COMMON             / SRNAMC / SRNAMT
 564 *     ..
 565 *     .. Data statements ..
 566       DATA            KTYPE / 1, 2, 5*4, 5*6, 3*9, 10 /
 567       DATA            KMAGN / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
 568       DATA            KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
 569      $                0, 0, 0 /
 570 *     ..
 571 *     .. Executable Statements ..
 572 *
 573 *     Check for errors
 574 *
 575       INFO = 0
 576 *
 577       BADMM = .FALSE.
 578       BADNN = .FALSE.
 579       MMAX = 1
 580       NMAX = 1
 581       MNMAX = 1
 582       MINWRK = 1
 583       DO 10 J = 1, NSIZES
 584          MMAX = MAX( MMAX, MVAL( J ) )
 585          IF( MVAL( J ).LT.0 )
 586      $      BADMM = .TRUE.
 587          NMAX = MAX( NMAX, NVAL( J ) )
 588          IF( NVAL( J ).LT.0 )
 589      $      BADNN = .TRUE.
 590          MNMAX = MAX( MNMAX, MIN( MVAL( J ), NVAL( J ) ) )
 591          MINWRK = MAX( MINWRK, 3*( MVAL( J )+NVAL( J ) ),
 592      $            MVAL( J )*( MVAL( J )+MAX( MVAL( J ), NVAL( J ),
 593      $            NRHS )+1 )+NVAL( J )*MIN( NVAL( J ), MVAL( J ) ) )
 594    10 CONTINUE
 595 *
 596 *     Check for errors
 597 *
 598       IF( NSIZES.LT.0 ) THEN
 599          INFO = -1
 600       ELSE IF( BADMM ) THEN
 601          INFO = -2
 602       ELSE IF( BADNN ) THEN
 603          INFO = -3
 604       ELSE IF( NTYPES.LT.0 ) THEN
 605          INFO = -4
 606       ELSE IF( NRHS.LT.0 ) THEN
 607          INFO = -6
 608       ELSE IF( LDA.LT.MMAX ) THEN
 609          INFO = -11
 610       ELSE IF( LDX.LT.MMAX ) THEN
 611          INFO = -17
 612       ELSE IF( LDQ.LT.MMAX ) THEN
 613          INFO = -21
 614       ELSE IF( LDPT.LT.MNMAX ) THEN
 615          INFO = -23
 616       ELSE IF( MINWRK.GT.LWORK ) THEN
 617          INFO = -27
 618       END IF
 619 *
 620       IF( INFO.NE.0 ) THEN
 621          CALL XERBLA( 'DCHKBD', -INFO )
 622          RETURN
 623       END IF
 624 *
 625 *     Initialize constants
 626 *
 627       PATH( 1: 1 ) = 'Double precision'
 628       PATH( 2: 3 ) = 'BD'
 629       NFAIL = 0
 630       NTEST = 0
 631       UNFL = DLAMCH( 'Safe minimum' )
 632       OVFL = DLAMCH( 'Overflow' )
 633       CALL DLABAD( UNFL, OVFL )
 634       ULP = DLAMCH( 'Precision' )
 635       ULPINV = ONE / ULP
 636       LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) )
 637       RTUNFL = SQRT( UNFL )
 638       RTOVFL = SQRT( OVFL )
 639       INFOT = 0
 640       ABSTOL = 2*UNFL
 641 *
 642 *     Loop over sizes, types
 643 *
 644       DO 300 JSIZE = 1, NSIZES
 645          M = MVAL( JSIZE )
 646          N = NVAL( JSIZE )
 647          MNMIN = MIN( M, N )
 648          AMNINV = ONE / MAX( M, N, 1 )
 649 *
 650          IF( NSIZES.NE.1 ) THEN
 651             MTYPES = MIN( MAXTYP, NTYPES )
 652          ELSE
 653             MTYPES = MIN( MAXTYP+1, NTYPES )
 654          END IF
 655 *
 656          DO 290 JTYPE = 1, MTYPES
 657             IF( .NOT.DOTYPE( JTYPE ) )
 658      $         GO TO 290
 659 *
 660             DO 20 J = 1, 4
 661                IOLDSD( J ) = ISEED( J )
 662    20       CONTINUE
 663 *
 664             DO 30 J = 1, 34
 665                RESULT( J ) = -ONE
 666    30       CONTINUE
 667 *
 668             UPLO = ' '
 669 *
 670 *           Compute "A"
 671 *
 672 *           Control parameters:
 673 *
 674 *           KMAGN  KMODE        KTYPE
 675 *       =1  O(1)   clustered 1  zero
 676 *       =2  large  clustered 2  identity
 677 *       =3  small  exponential  (none)
 678 *       =4         arithmetic   diagonal, (w/ eigenvalues)
 679 *       =5         random       symmetric, w/ eigenvalues
 680 *       =6                      nonsymmetric, w/ singular values
 681 *       =7                      random diagonal
 682 *       =8                      random symmetric
 683 *       =9                      random nonsymmetric
 684 *       =10                     random bidiagonal (log. distrib.)
 685 *
 686             IF( MTYPES.GT.MAXTYP )
 687      $         GO TO 100
 688 *
 689             ITYPE = KTYPE( JTYPE )
 690             IMODE = KMODE( JTYPE )
 691 *
 692 *           Compute norm
 693 *
 694             GO TO ( 40, 50, 60 )KMAGN( JTYPE )
 695 *
 696    40       CONTINUE
 697             ANORM = ONE
 698             GO TO 70
 699 *
 700    50       CONTINUE
 701             ANORM = ( RTOVFL*ULP )*AMNINV
 702             GO TO 70
 703 *
 704    60       CONTINUE
 705             ANORM = RTUNFL*MAX( M, N )*ULPINV
 706             GO TO 70
 707 *
 708    70       CONTINUE
 709 *
 710             CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
 711             IINFO = 0
 712             COND = ULPINV
 713 *
 714             BIDIAG = .FALSE.
 715             IF( ITYPE.EQ.1 ) THEN
 716 *
 717 *              Zero matrix
 718 *
 719                IINFO = 0
 720 *
 721             ELSE IF( ITYPE.EQ.2 ) THEN
 722 *
 723 *              Identity
 724 *
 725                DO 80 JCOL = 1, MNMIN
 726                   A( JCOL, JCOL ) = ANORM
 727    80          CONTINUE
 728 *
 729             ELSE IF( ITYPE.EQ.4 ) THEN
 730 *
 731 *              Diagonal Matrix, [Eigen]values Specified
 732 *
 733                CALL DLATMS( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, IMODE,
 734      $                      COND, ANORM, 0, 0, 'N', A, LDA,
 735      $                      WORK( MNMIN+1 ), IINFO )
 736 *
 737             ELSE IF( ITYPE.EQ.5 ) THEN
 738 *
 739 *              Symmetric, eigenvalues specified
 740 *
 741                CALL DLATMS( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, IMODE,
 742      $                      COND, ANORM, M, N, 'N', A, LDA,
 743      $                      WORK( MNMIN+1 ), IINFO )
 744 *
 745             ELSE IF( ITYPE.EQ.6 ) THEN
 746 *
 747 *              Nonsymmetric, singular values specified
 748 *
 749                CALL DLATMS( M, N, 'S', ISEED, 'N', WORK, IMODE, COND,
 750      $                      ANORM, M, N, 'N', A, LDA, WORK( MNMIN+1 ),
 751      $                      IINFO )
 752 *
 753             ELSE IF( ITYPE.EQ.7 ) THEN
 754 *
 755 *              Diagonal, random entries
 756 *
 757                CALL DLATMR( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, 6, ONE,
 758      $                      ONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
 759      $                      WORK( 2*MNMIN+1 ), 1, ONE, 'N', IWORK, 0, 0,
 760      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
 761 *
 762             ELSE IF( ITYPE.EQ.8 ) THEN
 763 *
 764 *              Symmetric, random entries
 765 *
 766                CALL DLATMR( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, 6, ONE,
 767      $                      ONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
 768      $                      WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N,
 769      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
 770 *
 771             ELSE IF( ITYPE.EQ.9 ) THEN
 772 *
 773 *              Nonsymmetric, random entries
 774 *
 775                CALL DLATMR( M, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
 776      $                      'T', 'N', WORK( MNMIN+1 ), 1, ONE,
 777      $                      WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N,
 778      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
 779 *
 780             ELSE IF( ITYPE.EQ.10 ) THEN
 781 *
 782 *              Bidiagonal, random entries
 783 *
 784                TEMP1 = -TWO*LOG( ULP )
 785                DO 90 J = 1, MNMIN
 786                   BD( J ) = EXP( TEMP1*DLARND( 2, ISEED ) )
 787                   IF( J.LT.MNMIN )
 788      $               BE( J ) = EXP( TEMP1*DLARND( 2, ISEED ) )
 789    90          CONTINUE
 790 *
 791                IINFO = 0
 792                BIDIAG = .TRUE.
 793                IF( M.GE.N ) THEN
 794                   UPLO = 'U'
 795                ELSE
 796                   UPLO = 'L'
 797                END IF
 798             ELSE
 799                IINFO = 1
 800             END IF
 801 *
 802             IF( IINFO.EQ.0 ) THEN
 803 *
 804 *              Generate Right-Hand Side
 805 *
 806                IF( BIDIAG ) THEN
 807                   CALL DLATMR( MNMIN, NRHS, 'S', ISEED, 'N', WORK, 6,
 808      $                         ONE, ONE, 'T', 'N', WORK( MNMIN+1 ), 1,
 809      $                         ONE, WORK( 2*MNMIN+1 ), 1, ONE, 'N',
 810      $                         IWORK, MNMIN, NRHS, ZERO, ONE, 'NO', Y,
 811      $                         LDX, IWORK, IINFO )
 812                ELSE
 813                   CALL DLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE,
 814      $                         ONE, 'T', 'N', WORK( M+1 ), 1, ONE,
 815      $                         WORK( 2*M+1 ), 1, ONE, 'N', IWORK, M,
 816      $                         NRHS, ZERO, ONE, 'NO', X, LDX, IWORK,
 817      $                         IINFO )
 818                END IF
 819             END IF
 820 *
 821 *           Error Exit
 822 *
 823             IF( IINFO.NE.0 ) THEN
 824                WRITE( NOUT, FMT = 9998 )'Generator', IINFO, M, N,
 825      $            JTYPE, IOLDSD
 826                INFO = ABS( IINFO )
 827                RETURN
 828             END IF
 829 *
 830   100       CONTINUE
 831 *
 832 *           Call DGEBRD and DORGBR to compute B, Q, and P, do tests.
 833 *
 834             IF( .NOT.BIDIAG ) THEN
 835 *
 836 *              Compute transformations to reduce A to bidiagonal form:
 837 *              B := Q' * A * P.
 838 *
 839                CALL DLACPY( ' ', M, N, A, LDA, Q, LDQ )
 840                CALL DGEBRD( M, N, Q, LDQ, BD, BE, WORK, WORK( MNMIN+1 ),
 841      $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
 842 *
 843 *              Check error code from DGEBRD.
 844 *
 845                IF( IINFO.NE.0 ) THEN
 846                   WRITE( NOUT, FMT = 9998 )'DGEBRD', IINFO, M, N,
 847      $               JTYPE, IOLDSD
 848                   INFO = ABS( IINFO )
 849                   RETURN
 850                END IF
 851 *
 852                CALL DLACPY( ' ', M, N, Q, LDQ, PT, LDPT )
 853                IF( M.GE.N ) THEN
 854                   UPLO = 'U'
 855                ELSE
 856                   UPLO = 'L'
 857                END IF
 858 *
 859 *              Generate Q
 860 *
 861                MQ = M
 862                IF( NRHS.LE.0 )
 863      $            MQ = MNMIN
 864                CALL DORGBR( 'Q', M, MQ, N, Q, LDQ, WORK,
 865      $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
 866 *
 867 *              Check error code from DORGBR.
 868 *
 869                IF( IINFO.NE.0 ) THEN
 870                   WRITE( NOUT, FMT = 9998 )'DORGBR(Q)', IINFO, M, N,
 871      $               JTYPE, IOLDSD
 872                   INFO = ABS( IINFO )
 873                   RETURN
 874                END IF
 875 *
 876 *              Generate P'
 877 *
 878                CALL DORGBR( 'P', MNMIN, N, M, PT, LDPT, WORK( MNMIN+1 ),
 879      $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
 880 *
 881 *              Check error code from DORGBR.
 882 *
 883                IF( IINFO.NE.0 ) THEN
 884                   WRITE( NOUT, FMT = 9998 )'DORGBR(P)', IINFO, M, N,
 885      $               JTYPE, IOLDSD
 886                   INFO = ABS( IINFO )
 887                   RETURN
 888                END IF
 889 *
 890 *              Apply Q' to an M by NRHS matrix X:  Y := Q' * X.
 891 *
 892                CALL DGEMM( 'Transpose', 'No transpose', M, NRHS, M, ONE,
 893      $                     Q, LDQ, X, LDX, ZERO, Y, LDX )
 894 *
 895 *              Test 1:  Check the decomposition A := Q * B * PT
 896 *                   2:  Check the orthogonality of Q
 897 *                   3:  Check the orthogonality of PT
 898 *
 899                CALL DBDT01( M, N, 1, A, LDA, Q, LDQ, BD, BE, PT, LDPT,
 900      $                      WORK, RESULT( 1 ) )
 901                CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
 902      $                      RESULT( 2 ) )
 903                CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
 904      $                      RESULT( 3 ) )
 905             END IF
 906 *
 907 *           Use DBDSQR to form the SVD of the bidiagonal matrix B:
 908 *           B := U * S1 * VT, and compute Z = U' * Y.
 909 *
 910             CALL DCOPY( MNMIN, BD, 1, S1, 1 )
 911             IF( MNMIN.GT.0 )
 912      $         CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
 913             CALL DLACPY( ' ', M, NRHS, Y, LDX, Z, LDX )
 914             CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT )
 915             CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT )
 916 *
 917             CALL DBDSQR( UPLO, MNMIN, MNMIN, MNMIN, NRHS, S1, WORK, VT,
 918      $                   LDPT, U, LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO )
 919 *
 920 *           Check error code from DBDSQR.
 921 *
 922             IF( IINFO.NE.0 ) THEN
 923                WRITE( NOUT, FMT = 9998 )'DBDSQR(vects)', IINFO, M, N,
 924      $            JTYPE, IOLDSD
 925                INFO = ABS( IINFO )
 926                IF( IINFO.LT.0 ) THEN
 927                   RETURN
 928                ELSE
 929                   RESULT( 4 ) = ULPINV
 930                   GO TO 270
 931                END IF
 932             END IF
 933 *
 934 *           Use DBDSQR to compute only the singular values of the
 935 *           bidiagonal matrix B;  U, VT, and Z should not be modified.
 936 *
 937             CALL DCOPY( MNMIN, BD, 1, S2, 1 )
 938             IF( MNMIN.GT.0 )
 939      $         CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
 940 *
 941             CALL DBDSQR( UPLO, MNMIN, 0, 0, 0, S2, WORK, VT, LDPT, U,
 942      $                   LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO )
 943 *
 944 *           Check error code from DBDSQR.
 945 *
 946             IF( IINFO.NE.0 ) THEN
 947                WRITE( NOUT, FMT = 9998 )'DBDSQR(values)', IINFO, M, N,
 948      $            JTYPE, IOLDSD
 949                INFO = ABS( IINFO )
 950                IF( IINFO.LT.0 ) THEN
 951                   RETURN
 952                ELSE
 953                   RESULT( 9 ) = ULPINV
 954                   GO TO 270
 955                END IF
 956             END IF
 957 *
 958 *           Test 4:  Check the decomposition B := U * S1 * VT
 959 *                5:  Check the computation Z := U' * Y
 960 *                6:  Check the orthogonality of U
 961 *                7:  Check the orthogonality of VT
 962 *
 963             CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT,
 964      $                   WORK, RESULT( 4 ) )
 965             CALL DBDT02( MNMIN, NRHS, Y, LDX, Z, LDX, U, LDPT, WORK,
 966      $                   RESULT( 5 ) )
 967             CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK,
 968      $                   RESULT( 6 ) )
 969             CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK,
 970      $                   RESULT( 7 ) )
 971 *
 972 *           Test 8:  Check that the singular values are sorted in
 973 *                    non-increasing order and are non-negative
 974 *
 975             RESULT( 8 ) = ZERO
 976             DO 110 I = 1, MNMIN - 1
 977                IF( S1( I ).LT.S1( I+1 ) )
 978      $            RESULT( 8 ) = ULPINV
 979                IF( S1( I ).LT.ZERO )
 980      $            RESULT( 8 ) = ULPINV
 981   110       CONTINUE
 982             IF( MNMIN.GE.1 ) THEN
 983                IF( S1( MNMIN ).LT.ZERO )
 984      $            RESULT( 8 ) = ULPINV
 985             END IF
 986 *
 987 *           Test 9:  Compare DBDSQR with and without singular vectors
 988 *
 989             TEMP2 = ZERO
 990 *
 991             DO 120 J = 1, MNMIN
 992                TEMP1 = ABS( S1( J )-S2( J ) ) /
 993      $                 MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
 994      $                 ULP*MAX( ABS( S1( J ) ), ABS( S2( J ) ) ) )
 995                TEMP2 = MAX( TEMP1, TEMP2 )
 996   120       CONTINUE
 997 *
 998             RESULT( 9 ) = TEMP2
 999 *
1000 *           Test 10:  Sturm sequence test of singular values
1001 *                     Go up by factors of two until it succeeds
1002 *
1003             TEMP1 = THRESH*( HALF-ULP )
1004 *
1005             DO 130 J = 0, LOG2UI
1006 *               CALL DSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO )
1007                IF( IINFO.EQ.0 )
1008      $            GO TO 140
1009                TEMP1 = TEMP1*TWO
1010   130       CONTINUE
1011 *
1012   140       CONTINUE
1013             RESULT( 10 ) = TEMP1
1014 *
1015 *           Use DBDSQR to form the decomposition A := (QU) S (VT PT)
1016 *           from the bidiagonal form A := Q B PT.
1017 *
1018             IF( .NOT.BIDIAG ) THEN
1019                CALL DCOPY( MNMIN, BD, 1, S2, 1 )
1020                IF( MNMIN.GT.0 )
1021      $            CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
1022 *
1023                CALL DBDSQR( UPLO, MNMIN, N, M, NRHS, S2, WORK, PT, LDPT,
1024      $                      Q, LDQ, Y, LDX, WORK( MNMIN+1 ), IINFO )
1025 *
1026 *              Test 11:  Check the decomposition A := Q*U * S2 * VT*PT
1027 *                   12:  Check the computation Z := U' * Q' * X
1028 *                   13:  Check the orthogonality of Q*U
1029 *                   14:  Check the orthogonality of VT*PT
1030 *
1031                CALL DBDT01( M, N, 0, A, LDA, Q, LDQ, S2, DUMMA, PT,
1032      $                      LDPT, WORK, RESULT( 11 ) )
1033                CALL DBDT02( M, NRHS, X, LDX, Y, LDX, Q, LDQ, WORK,
1034      $                      RESULT( 12 ) )
1035                CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
1036      $                      RESULT( 13 ) )
1037                CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
1038      $                      RESULT( 14 ) )
1039             END IF
1040 *
1041 *           Use DBDSDC to form the SVD of the bidiagonal matrix B:
1042 *           B := U * S1 * VT
1043 *
1044             CALL DCOPY( MNMIN, BD, 1, S1, 1 )
1045             IF( MNMIN.GT.0 )
1046      $         CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
1047             CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT )
1048             CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT )
1049 *
1050             CALL DBDSDC( UPLO, 'I', MNMIN, S1, WORK, U, LDPT, VT, LDPT,
1051      $                   DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO )
1052 *
1053 *           Check error code from DBDSDC.
1054 *
1055             IF( IINFO.NE.0 ) THEN
1056                WRITE( NOUT, FMT = 9998 )'DBDSDC(vects)', IINFO, M, N,
1057      $            JTYPE, IOLDSD
1058                INFO = ABS( IINFO )
1059                IF( IINFO.LT.0 ) THEN
1060                   RETURN
1061                ELSE
1062                   RESULT( 15 ) = ULPINV
1063                   GO TO 270
1064                END IF
1065             END IF
1066 *
1067 *           Use DBDSDC to compute only the singular values of the
1068 *           bidiagonal matrix B;  U and VT should not be modified.
1069 *
1070             CALL DCOPY( MNMIN, BD, 1, S2, 1 )
1071             IF( MNMIN.GT.0 )
1072      $         CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
1073 *
1074             CALL DBDSDC( UPLO, 'N', MNMIN, S2, WORK, DUM, 1, DUM, 1,
1075      $                   DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO )
1076 *
1077 *           Check error code from DBDSDC.
1078 *
1079             IF( IINFO.NE.0 ) THEN
1080                WRITE( NOUT, FMT = 9998 )'DBDSDC(values)', IINFO, M, N,
1081      $            JTYPE, IOLDSD
1082                INFO = ABS( IINFO )
1083                IF( IINFO.LT.0 ) THEN
1084                   RETURN
1085                ELSE
1086                   RESULT( 18 ) = ULPINV
1087                   GO TO 270
1088                END IF
1089             END IF
1090 *
1091 *           Test 15:  Check the decomposition B := U * S1 * VT
1092 *                16:  Check the orthogonality of U
1093 *                17:  Check the orthogonality of VT
1094 *
1095             CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT,
1096      $                   WORK, RESULT( 15 ) )
1097             CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK,
1098      $                   RESULT( 16 ) )
1099             CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK,
1100      $                   RESULT( 17 ) )
1101 *
1102 *           Test 18:  Check that the singular values are sorted in
1103 *                     non-increasing order and are non-negative
1104 *
1105             RESULT( 18 ) = ZERO
1106             DO 150 I = 1, MNMIN - 1
1107                IF( S1( I ).LT.S1( I+1 ) )
1108      $            RESULT( 18 ) = ULPINV
1109                IF( S1( I ).LT.ZERO )
1110      $            RESULT( 18 ) = ULPINV
1111   150       CONTINUE
1112             IF( MNMIN.GE.1 ) THEN
1113                IF( S1( MNMIN ).LT.ZERO )
1114      $            RESULT( 18 ) = ULPINV
1115             END IF
1116 *
1117 *           Test 19:  Compare DBDSQR with and without singular vectors
1118 *
1119             TEMP2 = ZERO
1120 *
1121             DO 160 J = 1, MNMIN
1122                TEMP1 = ABS( S1( J )-S2( J ) ) /
1123      $                 MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
1124      $                 ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) )
1125                TEMP2 = MAX( TEMP1, TEMP2 )
1126   160       CONTINUE
1127 *
1128             RESULT( 19 ) = TEMP2
1129 *
1130 *
1131 *           Use DBDSVDX to compute the SVD of the bidiagonal matrix B:
1132 *           B := U * S1 * VT
1133 *
1134             IF( JTYPE.EQ.10 .OR. JTYPE.EQ.16 ) THEN
1135 *              =================================
1136 *              Matrix types temporarily disabled
1137 *              =================================
1138                RESULT( 20:34 ) = ZERO
1139                GO TO 270
1140             END IF
1141 *
1142             IWBS = 1
1143             IWBD = IWBS + MNMIN
1144             IWBE = IWBD + MNMIN
1145             IWBZ = IWBE + MNMIN
1146             IWWORK = IWBZ + 2*MNMIN*(MNMIN+1)
1147             MNMIN2 = MAX( 1,MNMIN*2 )
1148 *
1149             CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 )
1150             IF( MNMIN.GT.0 )
1151      $         CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 )
1152 *
1153             CALL DBDSVDX( UPLO, 'V', 'A', MNMIN, WORK( IWBD ),
1154      $                    WORK( IWBE ), ZERO, ZERO, 0, 0, NS1, S1,
1155      $                    WORK( IWBZ ), MNMIN2, WORK( IWWORK ),
1156      $                    IWORK, IINFO)
1157 *
1158 *           Check error code from DBDSVDX.
1159 *
1160             IF( IINFO.NE.0 ) THEN
1161                WRITE( NOUT, FMT = 9998 )'DBDSVDX(vects,A)', IINFO, M, N,
1162      $            JTYPE, IOLDSD
1163                INFO = ABS( IINFO )
1164                IF( IINFO.LT.0 ) THEN
1165                   RETURN
1166                ELSE
1167                   RESULT( 20 ) = ULPINV
1168                   GO TO 270
1169                END IF
1170             END IF
1171 *
1172             J = IWBZ
1173             DO 170 I = 1, NS1
1174                CALL DCOPY( MNMIN, WORK( J ), 1, U( 1,I ), 1 )
1175                J = J + MNMIN
1176                CALL DCOPY( MNMIN, WORK( J ), 1, VT( I,1 ), LDPT )
1177                J = J + MNMIN
1178   170       CONTINUE
1179 *
1180 *           Use DBDSVDX to compute only the singular values of the
1181 *           bidiagonal matrix B;  U and VT should not be modified.
1182 *
1183             IF( JTYPE.EQ.9 ) THEN
1184 *              =================================
1185 *              Matrix types temporarily disabled
1186 *              =================================
1187                RESULT( 24 ) = ZERO
1188                GO TO 270
1189             END IF
1190 *
1191             CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 )
1192             IF( MNMIN.GT.0 )
1193      $         CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 )
1194 *
1195             CALL DBDSVDX( UPLO, 'N', 'A', MNMIN, WORK( IWBD ),
1196      $                    WORK( IWBE ), ZERO, ZERO, 0, 0, NS2, S2,
1197      $                    WORK( IWBZ ), MNMIN2, WORK( IWWORK ),
1198      $                    IWORK, IINFO )
1199 *
1200 *           Check error code from DBDSVDX.
1201 *
1202             IF( IINFO.NE.0 ) THEN
1203                WRITE( NOUT, FMT = 9998 )'DBDSVDX(values,A)', IINFO,
1204      $            M, N, JTYPE, IOLDSD
1205                INFO = ABS( IINFO )
1206                IF( IINFO.LT.0 ) THEN
1207                   RETURN
1208                ELSE
1209                   RESULT( 24 ) = ULPINV
1210                   GO TO 270
1211                END IF
1212             END IF
1213 *
1214 *           Save S1 for tests 30-34.
1215 *
1216             CALL DCOPY( MNMIN, S1, 1, WORK( IWBS ), 1 )
1217 *
1218 *           Test 20:  Check the decomposition B := U * S1 * VT
1219 *                21:  Check the orthogonality of U
1220 *                22:  Check the orthogonality of VT
1221 *                23:  Check that the singular values are sorted in
1222 *                     non-increasing order and are non-negative
1223 *                24:  Compare DBDSVDX with and without singular vectors
1224 *
1225             CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT,
1226      $                   LDPT, WORK( IWBS+MNMIN ), RESULT( 20 ) )
1227             CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT,
1228      $                   WORK( IWBS+MNMIN ), LWORK-MNMIN,
1229      $                   RESULT( 21 ) )
1230             CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT,
1231      $                   WORK( IWBS+MNMIN ), LWORK-MNMIN,
1232      $                   RESULT( 22) )
1233 *
1234             RESULT( 23 ) = ZERO
1235             DO 180 I = 1, MNMIN - 1
1236                IF( S1( I ).LT.S1( I+1 ) )
1237      $            RESULT( 23 ) = ULPINV
1238                IF( S1( I ).LT.ZERO )
1239      $            RESULT( 23 ) = ULPINV
1240   180       CONTINUE
1241             IF( MNMIN.GE.1 ) THEN
1242                IF( S1( MNMIN ).LT.ZERO )
1243      $            RESULT( 23 ) = ULPINV
1244             END IF
1245 *
1246             TEMP2 = ZERO
1247             DO 190 J = 1, MNMIN
1248                TEMP1 = ABS( S1( J )-S2( J ) ) /
1249      $                 MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
1250      $                 ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) )
1251                TEMP2 = MAX( TEMP1, TEMP2 )
1252   190       CONTINUE
1253             RESULT( 24 ) = TEMP2
1254             ANORM = S1( 1 )
1255 *
1256 *           Use DBDSVDX with RANGE='I': choose random values for IL and
1257 *           IU, and ask for the IL-th through IU-th singular values
1258 *           and corresponding vectors.
1259 *
1260             DO 200 I = 1, 4
1261                ISEED2( I ) = ISEED( I )
1262   200       CONTINUE
1263             IF( MNMIN.LE.1 ) THEN
1264                IL = 1
1265                IU = MNMIN
1266             ELSE
1267                IL = 1 + INT( ( MNMIN-1 )*DLARND( 1, ISEED2 ) )
1268                IU = 1 + INT( ( MNMIN-1 )*DLARND( 1, ISEED2 ) )
1269                IF( IU.LT.IL ) THEN
1270                   ITEMP = IU
1271                   IU = IL
1272                   IL = ITEMP
1273                END IF
1274             END IF
1275 *
1276             CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 )
1277             IF( MNMIN.GT.0 )
1278      $         CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 )
1279 *
1280             CALL DBDSVDX( UPLO, 'V', 'I', MNMIN, WORK( IWBD ),
1281      $                    WORK( IWBE ), ZERO, ZERO, IL, IU, NS1, S1,
1282      $                    WORK( IWBZ ), MNMIN2, WORK( IWWORK ),
1283      $                    IWORK, IINFO)
1284 *
1285 *           Check error code from DBDSVDX.
1286 *
1287             IF( IINFO.NE.0 ) THEN
1288                WRITE( NOUT, FMT = 9998 )'DBDSVDX(vects,I)', IINFO,
1289      $            M, N, JTYPE, IOLDSD
1290                INFO = ABS( IINFO )
1291                IF( IINFO.LT.0 ) THEN
1292                   RETURN
1293                ELSE
1294                   RESULT( 25 ) = ULPINV
1295                   GO TO 270
1296                END IF
1297             END IF
1298 *
1299             J = IWBZ
1300             DO 210 I = 1, NS1
1301                CALL DCOPY( MNMIN, WORK( J ), 1, U( 1,I ), 1 )
1302                J = J + MNMIN
1303                CALL DCOPY( MNMIN, WORK( J ), 1, VT( I,1 ), LDPT )
1304                J = J + MNMIN
1305   210       CONTINUE
1306 *
1307 *           Use DBDSVDX to compute only the singular values of the
1308 *           bidiagonal matrix B;  U and VT should not be modified.
1309 *
1310             CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 )
1311             IF( MNMIN.GT.0 )
1312      $         CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 )
1313 *
1314             CALL DBDSVDX( UPLO, 'N', 'I', MNMIN, WORK( IWBD ),
1315      $                    WORK( IWBE ), ZERO, ZERO, IL, IU, NS2, S2,
1316      $                    WORK( IWBZ ), MNMIN2, WORK( IWWORK ),
1317      $                    IWORK, IINFO )
1318 *
1319 *           Check error code from DBDSVDX.
1320 *
1321             IF( IINFO.NE.0 ) THEN
1322                WRITE( NOUT, FMT = 9998 )'DBDSVDX(values,I)', IINFO,
1323      $            M, N, JTYPE, IOLDSD
1324                INFO = ABS( IINFO )
1325                IF( IINFO.LT.0 ) THEN
1326                   RETURN
1327                ELSE
1328                   RESULT( 29 ) = ULPINV
1329                   GO TO 270
1330                END IF
1331             END IF
1332 *
1333 *           Test 25:  Check S1 - U' * B * VT'
1334 *                26:  Check the orthogonality of U
1335 *                27:  Check the orthogonality of VT
1336 *                28:  Check that the singular values are sorted in
1337 *                     non-increasing order and are non-negative
1338 *                29:  Compare DBDSVDX with and without singular vectors
1339 *
1340             CALL DBDT04( UPLO, MNMIN, BD, BE, S1, NS1, U,
1341      $                   LDPT, VT, LDPT, WORK( IWBS+MNMIN ),
1342      $                   RESULT( 25 ) )
1343             CALL DORT01( 'Columns', MNMIN, NS1, U, LDPT,
1344      $                   WORK( IWBS+MNMIN ), LWORK-MNMIN,
1345      $                   RESULT( 26 ) )
1346             CALL DORT01( 'Rows', NS1, MNMIN, VT, LDPT,
1347      $                   WORK( IWBS+MNMIN ), LWORK-MNMIN,
1348      $                   RESULT( 27 ) )
1349 *
1350             RESULT( 28 ) = ZERO
1351             DO 220 I = 1, NS1 - 1
1352                IF( S1( I ).LT.S1( I+1 ) )
1353      $            RESULT( 28 ) = ULPINV
1354                IF( S1( I ).LT.ZERO )
1355      $            RESULT( 28 ) = ULPINV
1356   220       CONTINUE
1357             IF( NS1.GE.1 ) THEN
1358                IF( S1( NS1 ).LT.ZERO )
1359      $            RESULT( 28 ) = ULPINV
1360             END IF
1361 *
1362             TEMP2 = ZERO
1363             DO 230 J = 1, NS1
1364                TEMP1 = ABS( S1( J )-S2( J ) ) /
1365      $                 MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
1366      $                 ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) )
1367                TEMP2 = MAX( TEMP1, TEMP2 )
1368   230       CONTINUE
1369             RESULT( 29 ) = TEMP2
1370 *
1371 *           Use DBDSVDX with RANGE='V': determine the values VL and VU
1372 *           of the IL-th and IU-th singular values and ask for all
1373 *           singular values in this range.
1374 *
1375             CALL DCOPY( MNMIN, WORK( IWBS ), 1, S1, 1 )
1376 *
1377             IF( MNMIN.GT.0 ) THEN
1378                IF( IL.NE.1 ) THEN
1379                   VU = S1( IL ) + MAX( HALF*ABS( S1( IL )-S1( IL-1 ) ),
1380      $                 ULP*ANORM, TWO*RTUNFL )
1381                ELSE
1382                   VU = S1( 1 ) + MAX( HALF*ABS( S1( MNMIN )-S1( 1 ) ),
1383      $                 ULP*ANORM, TWO*RTUNFL )
1384                END IF
1385                IF( IU.NE.NS1 ) THEN
1386                   VL = S1( IU ) - MAX( ULP*ANORM, TWO*RTUNFL,
1387      $                 HALF*ABS( S1( IU+1 )-S1( IU ) ) )
1388                ELSE
1389                   VL = S1( NS1 ) - MAX( ULP*ANORM, TWO*RTUNFL,
1390      $                 HALF*ABS( S1( MNMIN )-S1( 1 ) ) )
1391                END IF
1392                VL = MAX( VL,ZERO )
1393                VU = MAX( VU,ZERO )
1394                IF( VL.GE.VU ) VU = MAX( VU*2, VU+VL+HALF )
1395             ELSE
1396                VL = ZERO
1397                VU = ONE
1398             END IF
1399 *
1400             CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 )
1401             IF( MNMIN.GT.0 )
1402      $         CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 )
1403 *
1404             CALL DBDSVDX( UPLO, 'V', 'V', MNMIN, WORK( IWBD ),
1405      $                    WORK( IWBE ), VL, VU, 0, 0, NS1, S1,
1406      $                    WORK( IWBZ ), MNMIN2, WORK( IWWORK ),
1407      $                    IWORK, IINFO )
1408 *
1409 *           Check error code from DBDSVDX.
1410 *
1411             IF( IINFO.NE.0 ) THEN
1412                WRITE( NOUT, FMT = 9998 )'DBDSVDX(vects,V)', IINFO,
1413      $            M, N, JTYPE, IOLDSD
1414                INFO = ABS( IINFO )
1415                IF( IINFO.LT.0 ) THEN
1416                   RETURN
1417                ELSE
1418                   RESULT( 30 ) = ULPINV
1419                   GO TO 270
1420                END IF
1421             END IF
1422 *
1423             J = IWBZ
1424             DO 240 I = 1, NS1
1425                CALL DCOPY( MNMIN, WORK( J ), 1, U( 1,I ), 1 )
1426                J = J + MNMIN
1427                CALL DCOPY( MNMIN, WORK( J ), 1, VT( I,1 ), LDPT )
1428                J = J + MNMIN
1429   240       CONTINUE
1430 *
1431 *           Use DBDSVDX to compute only the singular values of the
1432 *           bidiagonal matrix B;  U and VT should not be modified.
1433 *
1434             CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 )
1435             IF( MNMIN.GT.0 )
1436      $         CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 )
1437 *
1438             CALL DBDSVDX( UPLO, 'N', 'V', MNMIN, WORK( IWBD ),
1439      $                    WORK( IWBE ), VL, VU, 0, 0, NS2, S2,
1440      $                    WORK( IWBZ ), MNMIN2, WORK( IWWORK ),
1441      $                    IWORK, IINFO )
1442 *
1443 *           Check error code from DBDSVDX.
1444 *
1445             IF( IINFO.NE.0 ) THEN
1446                WRITE( NOUT, FMT = 9998 )'DBDSVDX(values,V)', IINFO,
1447      $            M, N, JTYPE, IOLDSD
1448                INFO = ABS( IINFO )
1449                IF( IINFO.LT.0 ) THEN
1450                   RETURN
1451                ELSE
1452                   RESULT( 34 ) = ULPINV
1453                   GO TO 270
1454                END IF
1455             END IF
1456 *
1457 *           Test 30:  Check S1 - U' * B * VT'
1458 *                31:  Check the orthogonality of U
1459 *                32:  Check the orthogonality of VT
1460 *                33:  Check that the singular values are sorted in
1461 *                     non-increasing order and are non-negative
1462 *                34:  Compare DBDSVDX with and without singular vectors
1463 *
1464             CALL DBDT04( UPLO, MNMIN, BD, BE, S1, NS1, U,
1465      $                   LDPT, VT, LDPT, WORK( IWBS+MNMIN ),
1466      $                   RESULT( 30 ) )
1467             CALL DORT01( 'Columns', MNMIN, NS1, U, LDPT,
1468      $                   WORK( IWBS+MNMIN ), LWORK-MNMIN,
1469      $                   RESULT( 31 ) )
1470             CALL DORT01( 'Rows', NS1, MNMIN, VT, LDPT,
1471      $                   WORK( IWBS+MNMIN ), LWORK-MNMIN,
1472      $                   RESULT( 32 ) )
1473 *
1474             RESULT( 33 ) = ZERO
1475             DO 250 I = 1, NS1 - 1
1476                IF( S1( I ).LT.S1( I+1 ) )
1477      $            RESULT( 28 ) = ULPINV
1478                IF( S1( I ).LT.ZERO )
1479      $            RESULT( 28 ) = ULPINV
1480   250       CONTINUE
1481             IF( NS1.GE.1 ) THEN
1482                IF( S1( NS1 ).LT.ZERO )
1483      $            RESULT( 28 ) = ULPINV
1484             END IF
1485 *
1486             TEMP2 = ZERO
1487             DO 260 J = 1, NS1
1488                TEMP1 = ABS( S1( J )-S2( J ) ) /
1489      $                 MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
1490      $                 ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) )
1491                TEMP2 = MAX( TEMP1, TEMP2 )
1492   260       CONTINUE
1493             RESULT( 34 ) = TEMP2
1494 *
1495 *           End of Loop -- Check for RESULT(j) > THRESH
1496 *
1497   270       CONTINUE
1498 *
1499             DO 280 J = 1, 34
1500                IF( RESULT( J ).GE.THRESH ) THEN
1501                   IF( NFAIL.EQ.0 )
1502      $               CALL DLAHD2( NOUT, PATH )
1503                   WRITE( NOUT, FMT = 9999 )M, N, JTYPE, IOLDSD, J,
1504      $               RESULT( J )
1505                   NFAIL = NFAIL + 1
1506                END IF
1507   280       CONTINUE
1508             IF( .NOT.BIDIAG ) THEN
1509                NTEST = NTEST + 34
1510             ELSE
1511                NTEST = NTEST + 30
1512             END IF
1513 *
1514   290    CONTINUE
1515   300 CONTINUE
1516 *
1517 *     Summary
1518 *
1519       CALL ALASUM( PATH, NOUT, NFAIL, NTEST, 0 )
1520 *
1521       RETURN
1522 *
1523 *     End of DCHKBD
1524 *
1525  9999 FORMAT( ' M=', I5, ', N=', I5, ', type ', I2, ', seed=',
1526      $      4( I4, ',' ), ' test(', I2, ')=', G11.4 )
1527  9998 FORMAT( ' DCHKBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
1528      $      I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
1529      $      I5, ')' )
1530 *
1531       END