LAPACK: TESTING/EIG/dchkbd.f Source File

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