LAPACK: TESTING/EIG/zchkbd.f Source File

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