LAPACK: TESTING/MATGEN/zlatme.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI, 
00002      $  RSIGN, 
00003      $                   UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, 
00004      $  A, 
00005      $                   LDA, WORK, INFO )
00006 *
00007 *  -- LAPACK test routine (version 3.1) --
00008 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00009 *     June 2010
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          DIST, EI, RSIGN, SIM, UPPER
00013       INTEGER            INFO, KL, KU, LDA, MODE, MODES, N
00014       DOUBLE PRECISION   ANORM, COND, CONDS
00015       COMPLEX*16         DMAX
00016 *     ..
00017 *     .. Array Arguments ..
00018       INTEGER            ISEED( 4 )
00019       DOUBLE PRECISION   DS( * )
00020       COMPLEX*16         A( LDA, * ), D( * ), WORK( * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *     ZLATME generates random non-symmetric square matrices with
00027 *     specified eigenvalues for testing LAPACK programs.
00028 *
00029 *     ZLATME operates by applying the following sequence of
00030 *     operations:
00031 *
00032 *     1. Set the diagonal to D, where D may be input or
00033 *          computed according to MODE, COND, DMAX, and RSIGN
00034 *          as described below.
00035 *
00036 *     2. If UPPER='T', the upper triangle of A is set to random values
00037 *          out of distribution DIST.
00038 *
00039 *     3. If SIM='T', A is multiplied on the left by a random matrix
00040 *          X, whose singular values are specified by DS, MODES, and
00041 *          CONDS, and on the right by X inverse.
00042 *
00043 *     4. If KL < N-1, the lower bandwidth is reduced to KL using
00044 *          Householder transformations.  If KU < N-1, the upper
00045 *          bandwidth is reduced to KU.
00046 *
00047 *     5. If ANORM is not negative, the matrix is scaled to have
00048 *          maximum-element-norm ANORM.
00049 *
00050 *     (Note: since the matrix cannot be reduced beyond Hessenberg form,
00051 *      no packing options are available.)
00052 *
00053 *  Arguments
00054 *  =========
00055 *
00056 *  N        (input) INTEGER
00057 *           The number of columns (or rows) of A. Not modified.
00058 *
00059 *  DIST     (input) CHARACTER*1
00060 *           On entry, DIST specifies the type of distribution to be used
00061 *           to generate the random eigen-/singular values, and on the
00062 *           upper triangle (see UPPER).
00063 *           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
00064 *           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
00065 *           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
00066 *           'D' => uniform on the complex disc |z| < 1.
00067 *           Not modified.
00068 *
00069 *  ISEED    (input/output) INTEGER array, dimension ( 4 )
00070 *           On entry ISEED specifies the seed of the random number
00071 *           generator. They should lie between 0 and 4095 inclusive,
00072 *           and ISEED(4) should be odd. The random number generator
00073 *           uses a linear congruential sequence limited to small
00074 *           integers, and so should produce machine independent
00075 *           random numbers. The values of ISEED are changed on
00076 *           exit, and can be used in the next call to ZLATME
00077 *           to continue the same random number sequence.
00078 *           Changed on exit.
00079 *
00080 *  D        (input/output) COMPLEX*16 array, dimension ( N )
00081 *           This array is used to specify the eigenvalues of A.  If
00082 *           MODE=0, then D is assumed to contain the eigenvalues
00083 *           otherwise they will be computed according to MODE, COND,
00084 *           DMAX, and RSIGN and placed in D.
00085 *           Modified if MODE is nonzero.
00086 *
00087 *  MODE     (input) INTEGER
00088 *           On entry this describes how the eigenvalues are to
00089 *           be specified:
00090 *           MODE = 0 means use D as input
00091 *           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
00092 *           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
00093 *           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
00094 *           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
00095 *           MODE = 5 sets D to random numbers in the range
00096 *                    ( 1/COND , 1 ) such that their logarithms
00097 *                    are uniformly distributed.
00098 *           MODE = 6 set D to random numbers from same distribution
00099 *                    as the rest of the matrix.
00100 *           MODE < 0 has the same meaning as ABS(MODE), except that
00101 *              the order of the elements of D is reversed.
00102 *           Thus if MODE is between 1 and 4, D has entries ranging
00103 *              from 1 to 1/COND, if between -1 and -4, D has entries
00104 *              ranging from 1/COND to 1,
00105 *           Not modified.
00106 *
00107 *  COND     (input) DOUBLE PRECISION
00108 *           On entry, this is used as described under MODE above.
00109 *           If used, it must be >= 1. Not modified.
00110 *
00111 *  DMAX     (input) COMPLEX*16
00112 *           If MODE is neither -6, 0 nor 6, the contents of D, as
00113 *           computed according to MODE and COND, will be scaled by
00114 *           DMAX / max(abs(D(i))).  Note that DMAX need not be
00115 *           positive or real: if DMAX is negative or complex (or zero),
00116 *           D will be scaled by a negative or complex number (or zero).
00117 *           If RSIGN='F' then the largest (absolute) eigenvalue will be
00118 *           equal to DMAX.
00119 *           Not modified.
00120 *
00121 *  EI       (input) CHARACTER*1 array, dimension ( N )
00122 *           (ignored)
00123 *           Not modified.
00124 *
00125 *  RSIGN    (input) CHARACTER*1
00126 *           If MODE is not 0, 6, or -6, and RSIGN='T', then the
00127 *           elements of D, as computed according to MODE and COND, will
00128 *           be multiplied by a random complex number from the unit
00129 *           circle |z| = 1.  If RSIGN='F', they will not be.  RSIGN may
00130 *           only have the values 'T' or 'F'.
00131 *           Not modified.
00132 *
00133 *  UPPER    (input) CHARACTER*1
00134 *           If UPPER='T', then the elements of A above the diagonal
00135 *           will be set to random numbers out of DIST.  If UPPER='F',
00136 *           they will not.  UPPER may only have the values 'T' or 'F'.
00137 *           Not modified.
00138 *
00139 *  SIM      (input) CHARACTER*1
00140 *           If SIM='T', then A will be operated on by a "similarity
00141 *           transform", i.e., multiplied on the left by a matrix X and
00142 *           on the right by X inverse.  X = U S V, where U and V are
00143 *           random unitary matrices and S is a (diagonal) matrix of
00144 *           singular values specified by DS, MODES, and CONDS.  If
00145 *           SIM='F', then A will not be transformed.
00146 *           Not modified.
00147 *
00148 *  DS       (input/output) DOUBLE PRECISION array, dimension ( N )
00149 *           This array is used to specify the singular values of X,
00150 *           in the same way that D specifies the eigenvalues of A.
00151 *           If MODE=0, the DS contains the singular values, which
00152 *           may not be zero.
00153 *           Modified if MODE is nonzero.
00154 *
00155 *  MODES    (input) INTEGER
00156 *
00157 *  CONDS    (input) DOUBLE PRECISION
00158 *           Similar to MODE and COND, but for specifying the diagonal
00159 *           of S.  MODES=-6 and +6 are not allowed (since they would
00160 *           result in randomly ill-conditioned eigenvalues.)
00161 *
00162 *  KL       (input) INTEGER
00163 *           This specifies the lower bandwidth of the  matrix.  KL=1
00164 *           specifies upper Hessenberg form.  If KL is at least N-1,
00165 *           then A will have full lower bandwidth.
00166 *           Not modified.
00167 *
00168 *  KU       (input) INTEGER
00169 *           This specifies the upper bandwidth of the  matrix.  KU=1
00170 *           specifies lower Hessenberg form.  If KU is at least N-1,
00171 *           then A will have full upper bandwidth; if KU and KL
00172 *           are both at least N-1, then A will be dense.  Only one of
00173 *           KU and KL may be less than N-1.
00174 *           Not modified.
00175 *
00176 *  ANORM    (input) DOUBLE PRECISION
00177 *           If ANORM is not negative, then A will be scaled by a non-
00178 *           negative real number to make the maximum-element-norm of A
00179 *           to be ANORM.
00180 *           Not modified.
00181 *
00182 *  A        (output) COMPLEX*16 array, dimension ( LDA, N )
00183 *           On exit A is the desired test matrix.
00184 *           Modified.
00185 *
00186 *  LDA      (input) INTEGER
00187 *           LDA specifies the first dimension of A as declared in the
00188 *           calling program.  LDA must be at least M.
00189 *           Not modified.
00190 *
00191 *  WORK     (workspace) COMPLEX*16 array, dimension ( 3*N )
00192 *           Workspace.
00193 *           Modified.
00194 *
00195 *  INFO     (output) INTEGER
00196 *           Error code.  On exit, INFO will be set to one of the
00197 *           following values:
00198 *             0 => normal return
00199 *            -1 => N negative
00200 *            -2 => DIST illegal string
00201 *            -5 => MODE not in range -6 to 6
00202 *            -6 => COND less than 1.0, and MODE neither -6, 0 nor 6
00203 *            -9 => RSIGN is not 'T' or 'F'
00204 *           -10 => UPPER is not 'T' or 'F'
00205 *           -11 => SIM   is not 'T' or 'F'
00206 *           -12 => MODES=0 and DS has a zero singular value.
00207 *           -13 => MODES is not in the range -5 to 5.
00208 *           -14 => MODES is nonzero and CONDS is less than 1.
00209 *           -15 => KL is less than 1.
00210 *           -16 => KU is less than 1, or KL and KU are both less than
00211 *                  N-1.
00212 *           -19 => LDA is less than M.
00213 *            1  => Error return from ZLATM1 (computing D)
00214 *            2  => Cannot scale to DMAX (max. eigenvalue is 0)
00215 *            3  => Error return from DLATM1 (computing DS)
00216 *            4  => Error return from ZLARGE
00217 *            5  => Zero singular value from DLATM1.
00218 *
00219 *  =====================================================================
00220 *
00221 *     .. Parameters ..
00222       DOUBLE PRECISION   ZERO
00223       PARAMETER          ( ZERO = 0.0D+0 )
00224       DOUBLE PRECISION   ONE
00225       PARAMETER          ( ONE = 1.0D+0 )
00226       COMPLEX*16         CZERO
00227       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
00228       COMPLEX*16         CONE
00229       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
00230 *     ..
00231 *     .. Local Scalars ..
00232       LOGICAL            BADS
00233       INTEGER            I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN,
00234      $                   ISIM, IUPPER, J, JC, JCR
00235       DOUBLE PRECISION   RALPHA, TEMP
00236       COMPLEX*16         ALPHA, TAU, XNORMS
00237 *     ..
00238 *     .. Local Arrays ..
00239       DOUBLE PRECISION   TEMPA( 1 )
00240 *     ..
00241 *     .. External Functions ..
00242       LOGICAL            LSAME
00243       DOUBLE PRECISION   ZLANGE
00244       COMPLEX*16         ZLARND
00245       EXTERNAL           LSAME, ZLANGE, ZLARND
00246 *     ..
00247 *     .. External Subroutines ..
00248       EXTERNAL           DLATM1, XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZGERC,
00249      $                   ZLACGV, ZLARFG, ZLARGE, ZLARNV, ZLASET, ZLATM1,
00250      $                   ZSCAL
00251 *     ..
00252 *     .. Intrinsic Functions ..
00253       INTRINSIC          ABS, DCONJG, MAX, MOD
00254 *     ..
00255 *     .. Executable Statements ..
00256 *
00257 *     1)      Decode and Test the input parameters.
00258 *             Initialize flags & seed.
00259 *
00260       INFO = 0
00261 *
00262 *     Quick return if possible
00263 *
00264       IF( N.EQ.0 )
00265      $   RETURN
00266 *
00267 *     Decode DIST
00268 *
00269       IF( LSAME( DIST, 'U' ) ) THEN
00270          IDIST = 1
00271       ELSE IF( LSAME( DIST, 'S' ) ) THEN
00272          IDIST = 2
00273       ELSE IF( LSAME( DIST, 'N' ) ) THEN
00274          IDIST = 3
00275       ELSE IF( LSAME( DIST, 'D' ) ) THEN
00276          IDIST = 4
00277       ELSE
00278          IDIST = -1
00279       END IF
00280 *
00281 *     Decode RSIGN
00282 *
00283       IF( LSAME( RSIGN, 'T' ) ) THEN
00284          IRSIGN = 1
00285       ELSE IF( LSAME( RSIGN, 'F' ) ) THEN
00286          IRSIGN = 0
00287       ELSE
00288          IRSIGN = -1
00289       END IF
00290 *
00291 *     Decode UPPER
00292 *
00293       IF( LSAME( UPPER, 'T' ) ) THEN
00294          IUPPER = 1
00295       ELSE IF( LSAME( UPPER, 'F' ) ) THEN
00296          IUPPER = 0
00297       ELSE
00298          IUPPER = -1
00299       END IF
00300 *
00301 *     Decode SIM
00302 *
00303       IF( LSAME( SIM, 'T' ) ) THEN
00304          ISIM = 1
00305       ELSE IF( LSAME( SIM, 'F' ) ) THEN
00306          ISIM = 0
00307       ELSE
00308          ISIM = -1
00309       END IF
00310 *
00311 *     Check DS, if MODES=0 and ISIM=1
00312 *
00313       BADS = .FALSE.
00314       IF( MODES.EQ.0 .AND. ISIM.EQ.1 ) THEN
00315          DO 10 J = 1, N
00316             IF( DS( J ).EQ.ZERO )
00317      $         BADS = .TRUE.
00318    10    CONTINUE
00319       END IF
00320 *
00321 *     Set INFO if an error
00322 *
00323       IF( N.LT.0 ) THEN
00324          INFO = -1
00325       ELSE IF( IDIST.EQ.-1 ) THEN
00326          INFO = -2
00327       ELSE IF( ABS( MODE ).GT.6 ) THEN
00328          INFO = -5
00329       ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE )
00330      $          THEN
00331          INFO = -6
00332       ELSE IF( IRSIGN.EQ.-1 ) THEN
00333          INFO = -9
00334       ELSE IF( IUPPER.EQ.-1 ) THEN
00335          INFO = -10
00336       ELSE IF( ISIM.EQ.-1 ) THEN
00337          INFO = -11
00338       ELSE IF( BADS ) THEN
00339          INFO = -12
00340       ELSE IF( ISIM.EQ.1 .AND. ABS( MODES ).GT.5 ) THEN
00341          INFO = -13
00342       ELSE IF( ISIM.EQ.1 .AND. MODES.NE.0 .AND. CONDS.LT.ONE ) THEN
00343          INFO = -14
00344       ELSE IF( KL.LT.1 ) THEN
00345          INFO = -15
00346       ELSE IF( KU.LT.1 .OR. ( KU.LT.N-1 .AND. KL.LT.N-1 ) ) THEN
00347          INFO = -16
00348       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00349          INFO = -19
00350       END IF
00351 *
00352       IF( INFO.NE.0 ) THEN
00353          CALL XERBLA( 'ZLATME', -INFO )
00354          RETURN
00355       END IF
00356 *
00357 *     Initialize random number generator
00358 *
00359       DO 20 I = 1, 4
00360          ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
00361    20 CONTINUE
00362 *
00363       IF( MOD( ISEED( 4 ), 2 ).NE.1 )
00364      $   ISEED( 4 ) = ISEED( 4 ) + 1
00365 *
00366 *     2)      Set up diagonal of A
00367 *
00368 *             Compute D according to COND and MODE
00369 *
00370       CALL ZLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, N, IINFO )
00371       IF( IINFO.NE.0 ) THEN
00372          INFO = 1
00373          RETURN
00374       END IF
00375       IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN
00376 *
00377 *        Scale by DMAX
00378 *
00379          TEMP = ABS( D( 1 ) )
00380          DO 30 I = 2, N
00381             TEMP = MAX( TEMP, ABS( D( I ) ) )
00382    30    CONTINUE
00383 *
00384          IF( TEMP.GT.ZERO ) THEN
00385             ALPHA = DMAX / TEMP
00386          ELSE
00387             INFO = 2
00388             RETURN
00389          END IF
00390 *
00391          CALL ZSCAL( N, ALPHA, D, 1 )
00392 *
00393       END IF
00394 *
00395       CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
00396       CALL ZCOPY( N, D, 1, A, LDA+1 )
00397 *
00398 *     3)      If UPPER='T', set upper triangle of A to random numbers.
00399 *
00400       IF( IUPPER.NE.0 ) THEN
00401          DO 40 JC = 2, N
00402             CALL ZLARNV( IDIST, ISEED, JC-1, A( 1, JC ) )
00403    40    CONTINUE
00404       END IF
00405 *
00406 *     4)      If SIM='T', apply similarity transformation.
00407 *
00408 *                                -1
00409 *             Transform is  X A X  , where X = U S V, thus
00410 *
00411 *             it is  U S V A V' (1/S) U'
00412 *
00413       IF( ISIM.NE.0 ) THEN
00414 *
00415 *        Compute S (singular values of the eigenvector matrix)
00416 *        according to CONDS and MODES
00417 *
00418          CALL DLATM1( MODES, CONDS, 0, 0, ISEED, DS, N, IINFO )
00419          IF( IINFO.NE.0 ) THEN
00420             INFO = 3
00421             RETURN
00422          END IF
00423 *
00424 *        Multiply by V and V'
00425 *
00426          CALL ZLARGE( N, A, LDA, ISEED, WORK, IINFO )
00427          IF( IINFO.NE.0 ) THEN
00428             INFO = 4
00429             RETURN
00430          END IF
00431 *
00432 *        Multiply by S and (1/S)
00433 *
00434          DO 50 J = 1, N
00435             CALL ZDSCAL( N, DS( J ), A( J, 1 ), LDA )
00436             IF( DS( J ).NE.ZERO ) THEN
00437                CALL ZDSCAL( N, ONE / DS( J ), A( 1, J ), 1 )
00438             ELSE
00439                INFO = 5
00440                RETURN
00441             END IF
00442    50    CONTINUE
00443 *
00444 *        Multiply by U and U'
00445 *
00446          CALL ZLARGE( N, A, LDA, ISEED, WORK, IINFO )
00447          IF( IINFO.NE.0 ) THEN
00448             INFO = 4
00449             RETURN
00450          END IF
00451       END IF
00452 *
00453 *     5)      Reduce the bandwidth.
00454 *
00455       IF( KL.LT.N-1 ) THEN
00456 *
00457 *        Reduce bandwidth -- kill column
00458 *
00459          DO 60 JCR = KL + 1, N - 1
00460             IC = JCR - KL
00461             IROWS = N + 1 - JCR
00462             ICOLS = N + KL - JCR
00463 *
00464             CALL ZCOPY( IROWS, A( JCR, IC ), 1, WORK, 1 )
00465             XNORMS = WORK( 1 )
00466             CALL ZLARFG( IROWS, XNORMS, WORK( 2 ), 1, TAU )
00467             TAU = DCONJG( TAU )
00468             WORK( 1 ) = CONE
00469             ALPHA = ZLARND( 5, ISEED )
00470 *
00471             CALL ZGEMV( 'C', IROWS, ICOLS, CONE, A( JCR, IC+1 ), LDA,
00472      $                  WORK, 1, CZERO, WORK( IROWS+1 ), 1 )
00473             CALL ZGERC( IROWS, ICOLS, -TAU, WORK, 1, WORK( IROWS+1 ), 1,
00474      $                  A( JCR, IC+1 ), LDA )
00475 *
00476             CALL ZGEMV( 'N', N, IROWS, CONE, A( 1, JCR ), LDA, WORK, 1,
00477      $                  CZERO, WORK( IROWS+1 ), 1 )
00478             CALL ZGERC( N, IROWS, -DCONJG( TAU ), WORK( IROWS+1 ), 1,
00479      $                  WORK, 1, A( 1, JCR ), LDA )
00480 *
00481             A( JCR, IC ) = XNORMS
00482             CALL ZLASET( 'Full', IROWS-1, 1, CZERO, CZERO,
00483      $                   A( JCR+1, IC ), LDA )
00484 *
00485             CALL ZSCAL( ICOLS+1, ALPHA, A( JCR, IC ), LDA )
00486             CALL ZSCAL( N, DCONJG( ALPHA ), A( 1, JCR ), 1 )
00487    60    CONTINUE
00488       ELSE IF( KU.LT.N-1 ) THEN
00489 *
00490 *        Reduce upper bandwidth -- kill a row at a time.
00491 *
00492          DO 70 JCR = KU + 1, N - 1
00493             IR = JCR - KU
00494             IROWS = N + KU - JCR
00495             ICOLS = N + 1 - JCR
00496 *
00497             CALL ZCOPY( ICOLS, A( IR, JCR ), LDA, WORK, 1 )
00498             XNORMS = WORK( 1 )
00499             CALL ZLARFG( ICOLS, XNORMS, WORK( 2 ), 1, TAU )
00500             TAU = DCONJG( TAU )
00501             WORK( 1 ) = CONE
00502             CALL ZLACGV( ICOLS-1, WORK( 2 ), 1 )
00503             ALPHA = ZLARND( 5, ISEED )
00504 *
00505             CALL ZGEMV( 'N', IROWS, ICOLS, CONE, A( IR+1, JCR ), LDA,
00506      $                  WORK, 1, CZERO, WORK( ICOLS+1 ), 1 )
00507             CALL ZGERC( IROWS, ICOLS, -TAU, WORK( ICOLS+1 ), 1, WORK, 1,
00508      $                  A( IR+1, JCR ), LDA )
00509 *
00510             CALL ZGEMV( 'C', ICOLS, N, CONE, A( JCR, 1 ), LDA, WORK, 1,
00511      $                  CZERO, WORK( ICOLS+1 ), 1 )
00512             CALL ZGERC( ICOLS, N, -DCONJG( TAU ), WORK, 1,
00513      $                  WORK( ICOLS+1 ), 1, A( JCR, 1 ), LDA )
00514 *
00515             A( IR, JCR ) = XNORMS
00516             CALL ZLASET( 'Full', 1, ICOLS-1, CZERO, CZERO,
00517      $                   A( IR, JCR+1 ), LDA )
00518 *
00519             CALL ZSCAL( IROWS+1, ALPHA, A( IR, JCR ), 1 )
00520             CALL ZSCAL( N, DCONJG( ALPHA ), A( JCR, 1 ), LDA )
00521    70    CONTINUE
00522       END IF
00523 *
00524 *     Scale the matrix to have norm ANORM
00525 *
00526       IF( ANORM.GE.ZERO ) THEN
00527          TEMP = ZLANGE( 'M', N, N, A, LDA, TEMPA )
00528          IF( TEMP.GT.ZERO ) THEN
00529             RALPHA = ANORM / TEMP
00530             DO 80 J = 1, N
00531                CALL ZDSCAL( N, RALPHA, A( 1, J ), 1 )
00532    80       CONTINUE
00533          END IF
00534       END IF
00535 *
00536       RETURN
00537 *
00538 *     End of ZLATME
00539 *
00540       END