SUBROUTINE ZLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI, RSIGN, \$ UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, A, \$ LDA, WORK, INFO ) * * -- LAPACK test routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER DIST, EI, RSIGN, SIM, UPPER INTEGER INFO, KL, KU, LDA, MODE, MODES, N DOUBLE PRECISION ANORM, COND, CONDS COMPLEX*16 DMAX * .. * .. Array Arguments .. INTEGER ISEED( 4 ) DOUBLE PRECISION DS( * ) COMPLEX*16 A( LDA, * ), D( * ), WORK( * ) * .. * * Purpose * ======= * * ZLATME generates random non-symmetric square matrices with * specified eigenvalues for testing LAPACK programs. * * ZLATME operates by applying the following sequence of * operations: * * 1. Set the diagonal to D, where D may be input or * computed according to MODE, COND, DMAX, and RSIGN * as described below. * * 2. If UPPER='T', the upper triangle of A is set to random values * out of distribution DIST. * * 3. If SIM='T', A is multiplied on the left by a random matrix * X, whose singular values are specified by DS, MODES, and * CONDS, and on the right by X inverse. * * 4. If KL < N-1, the lower bandwidth is reduced to KL using * Householder transformations. If KU < N-1, the upper * bandwidth is reduced to KU. * * 5. If ANORM is not negative, the matrix is scaled to have * maximum-element-norm ANORM. * * (Note: since the matrix cannot be reduced beyond Hessenberg form, * no packing options are available.) * * Arguments * ========= * * N - INTEGER * The number of columns (or rows) of A. Not modified. * * DIST - CHARACTER*1 * On entry, DIST specifies the type of distribution to be used * to generate the random eigen-/singular values, and on the * upper triangle (see UPPER). * 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform ) * 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric ) * 'N' => NORMAL( 0, 1 ) ( 'N' for normal ) * 'D' => uniform on the complex disc |z| < 1. * Not modified. * * ISEED - INTEGER array, dimension ( 4 ) * On entry ISEED specifies the seed of the random number * generator. They should lie between 0 and 4095 inclusive, * and ISEED(4) should be odd. The random number generator * uses a linear congruential sequence limited to small * integers, and so should produce machine independent * random numbers. The values of ISEED are changed on * exit, and can be used in the next call to ZLATME * to continue the same random number sequence. * Changed on exit. * * D - COMPLEX*16 array, dimension ( N ) * This array is used to specify the eigenvalues of A. If * MODE=0, then D is assumed to contain the eigenvalues * otherwise they will be computed according to MODE, COND, * DMAX, and RSIGN and placed in D. * Modified if MODE is nonzero. * * MODE - INTEGER * On entry this describes how the eigenvalues are to * be specified: * MODE = 0 means use D as input * MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND * MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND * MODE = 3 sets D(I)=COND**(-(I-1)/(N-1)) * MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND) * MODE = 5 sets D to random numbers in the range * ( 1/COND , 1 ) such that their logarithms * are uniformly distributed. * MODE = 6 set D to random numbers from same distribution * as the rest of the matrix. * MODE < 0 has the same meaning as ABS(MODE), except that * the order of the elements of D is reversed. * Thus if MODE is between 1 and 4, D has entries ranging * from 1 to 1/COND, if between -1 and -4, D has entries * ranging from 1/COND to 1, * Not modified. * * COND - DOUBLE PRECISION * On entry, this is used as described under MODE above. * If used, it must be >= 1. Not modified. * * DMAX - COMPLEX*16 * If MODE is neither -6, 0 nor 6, the contents of D, as * computed according to MODE and COND, will be scaled by * DMAX / max(abs(D(i))). Note that DMAX need not be * positive or real: if DMAX is negative or complex (or zero), * D will be scaled by a negative or complex number (or zero). * If RSIGN='F' then the largest (absolute) eigenvalue will be * equal to DMAX. * Not modified. * * EI - CHARACTER*1 (ignored) * Not modified. * * RSIGN - CHARACTER*1 * If MODE is not 0, 6, or -6, and RSIGN='T', then the * elements of D, as computed according to MODE and COND, will * be multiplied by a random complex number from the unit * circle |z| = 1. If RSIGN='F', they will not be. RSIGN may * only have the values 'T' or 'F'. * Not modified. * * UPPER - CHARACTER*1 * If UPPER='T', then the elements of A above the diagonal * will be set to random numbers out of DIST. If UPPER='F', * they will not. UPPER may only have the values 'T' or 'F'. * Not modified. * * SIM - CHARACTER*1 * If SIM='T', then A will be operated on by a "similarity * transform", i.e., multiplied on the left by a matrix X and * on the right by X inverse. X = U S V, where U and V are * random unitary matrices and S is a (diagonal) matrix of * singular values specified by DS, MODES, and CONDS. If * SIM='F', then A will not be transformed. * Not modified. * * DS - DOUBLE PRECISION array, dimension ( N ) * This array is used to specify the singular values of X, * in the same way that D specifies the eigenvalues of A. * If MODE=0, the DS contains the singular values, which * may not be zero. * Modified if MODE is nonzero. * * MODES - INTEGER * CONDS - DOUBLE PRECISION * Similar to MODE and COND, but for specifying the diagonal * of S. MODES=-6 and +6 are not allowed (since they would * result in randomly ill-conditioned eigenvalues.) * * KL - INTEGER * This specifies the lower bandwidth of the matrix. KL=1 * specifies upper Hessenberg form. If KL is at least N-1, * then A will have full lower bandwidth. * Not modified. * * KU - INTEGER * This specifies the upper bandwidth of the matrix. KU=1 * specifies lower Hessenberg form. If KU is at least N-1, * then A will have full upper bandwidth; if KU and KL * are both at least N-1, then A will be dense. Only one of * KU and KL may be less than N-1. * Not modified. * * ANORM - DOUBLE PRECISION * If ANORM is not negative, then A will be scaled by a non- * negative real number to make the maximum-element-norm of A * to be ANORM. * Not modified. * * A - COMPLEX*16 array, dimension ( LDA, N ) * On exit A is the desired test matrix. * Modified. * * LDA - INTEGER * LDA specifies the first dimension of A as declared in the * calling program. LDA must be at least M. * Not modified. * * WORK - COMPLEX*16 array, dimension ( 3*N ) * Workspace. * Modified. * * INFO - INTEGER * Error code. On exit, INFO will be set to one of the * following values: * 0 => normal return * -1 => N negative * -2 => DIST illegal string * -5 => MODE not in range -6 to 6 * -6 => COND less than 1.0, and MODE neither -6, 0 nor 6 * -9 => RSIGN is not 'T' or 'F' * -10 => UPPER is not 'T' or 'F' * -11 => SIM is not 'T' or 'F' * -12 => MODES=0 and DS has a zero singular value. * -13 => MODES is not in the range -5 to 5. * -14 => MODES is nonzero and CONDS is less than 1. * -15 => KL is less than 1. * -16 => KU is less than 1, or KL and KU are both less than * N-1. * -19 => LDA is less than M. * 1 => Error return from ZLATM1 (computing D) * 2 => Cannot scale to DMAX (max. eigenvalue is 0) * 3 => Error return from DLATM1 (computing DS) * 4 => Error return from ZLARGE * 5 => Zero singular value from DLATM1. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) COMPLEX*16 CZERO PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) COMPLEX*16 CONE PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) * .. * .. Local Scalars .. LOGICAL BADS INTEGER I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN, \$ ISIM, IUPPER, J, JC, JCR DOUBLE PRECISION RALPHA, TEMP COMPLEX*16 ALPHA, TAU, XNORMS * .. * .. Local Arrays .. DOUBLE PRECISION TEMPA( 1 ) * .. * .. External Functions .. LOGICAL LSAME DOUBLE PRECISION ZLANGE COMPLEX*16 ZLARND EXTERNAL LSAME, ZLANGE, ZLARND * .. * .. External Subroutines .. EXTERNAL DLATM1, XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZGERC, \$ ZLACGV, ZLARFG, ZLARGE, ZLARNV, ZLASET, ZLATM1, \$ ZSCAL * .. * .. Intrinsic Functions .. INTRINSIC ABS, DCONJG, MAX, MOD * .. * .. Executable Statements .. * * 1) Decode and Test the input parameters. * Initialize flags & seed. * INFO = 0 * * Quick return if possible * IF( N.EQ.0 ) \$ RETURN * * Decode DIST * IF( LSAME( DIST, 'U' ) ) THEN IDIST = 1 ELSE IF( LSAME( DIST, 'S' ) ) THEN IDIST = 2 ELSE IF( LSAME( DIST, 'N' ) ) THEN IDIST = 3 ELSE IF( LSAME( DIST, 'D' ) ) THEN IDIST = 4 ELSE IDIST = -1 END IF * * Decode RSIGN * IF( LSAME( RSIGN, 'T' ) ) THEN IRSIGN = 1 ELSE IF( LSAME( RSIGN, 'F' ) ) THEN IRSIGN = 0 ELSE IRSIGN = -1 END IF * * Decode UPPER * IF( LSAME( UPPER, 'T' ) ) THEN IUPPER = 1 ELSE IF( LSAME( UPPER, 'F' ) ) THEN IUPPER = 0 ELSE IUPPER = -1 END IF * * Decode SIM * IF( LSAME( SIM, 'T' ) ) THEN ISIM = 1 ELSE IF( LSAME( SIM, 'F' ) ) THEN ISIM = 0 ELSE ISIM = -1 END IF * * Check DS, if MODES=0 and ISIM=1 * BADS = .FALSE. IF( MODES.EQ.0 .AND. ISIM.EQ.1 ) THEN DO 10 J = 1, N IF( DS( J ).EQ.ZERO ) \$ BADS = .TRUE. 10 CONTINUE END IF * * Set INFO if an error * IF( N.LT.0 ) THEN INFO = -1 ELSE IF( IDIST.EQ.-1 ) THEN INFO = -2 ELSE IF( ABS( MODE ).GT.6 ) THEN INFO = -5 ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE ) \$ THEN INFO = -6 ELSE IF( IRSIGN.EQ.-1 ) THEN INFO = -9 ELSE IF( IUPPER.EQ.-1 ) THEN INFO = -10 ELSE IF( ISIM.EQ.-1 ) THEN INFO = -11 ELSE IF( BADS ) THEN INFO = -12 ELSE IF( ISIM.EQ.1 .AND. ABS( MODES ).GT.5 ) THEN INFO = -13 ELSE IF( ISIM.EQ.1 .AND. MODES.NE.0 .AND. CONDS.LT.ONE ) THEN INFO = -14 ELSE IF( KL.LT.1 ) THEN INFO = -15 ELSE IF( KU.LT.1 .OR. ( KU.LT.N-1 .AND. KL.LT.N-1 ) ) THEN INFO = -16 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -19 END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'ZLATME', -INFO ) RETURN END IF * * Initialize random number generator * DO 20 I = 1, 4 ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 ) 20 CONTINUE * IF( MOD( ISEED( 4 ), 2 ).NE.1 ) \$ ISEED( 4 ) = ISEED( 4 ) + 1 * * 2) Set up diagonal of A * * Compute D according to COND and MODE * CALL ZLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, N, IINFO ) IF( IINFO.NE.0 ) THEN INFO = 1 RETURN END IF IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN * * Scale by DMAX * TEMP = ABS( D( 1 ) ) DO 30 I = 2, N TEMP = MAX( TEMP, ABS( D( I ) ) ) 30 CONTINUE * IF( TEMP.GT.ZERO ) THEN ALPHA = DMAX / TEMP ELSE INFO = 2 RETURN END IF * CALL ZSCAL( N, ALPHA, D, 1 ) * END IF * CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) CALL ZCOPY( N, D, 1, A, LDA+1 ) * * 3) If UPPER='T', set upper triangle of A to random numbers. * IF( IUPPER.NE.0 ) THEN DO 40 JC = 2, N CALL ZLARNV( IDIST, ISEED, JC-1, A( 1, JC ) ) 40 CONTINUE END IF * * 4) If SIM='T', apply similarity transformation. * * -1 * Transform is X A X , where X = U S V, thus * * it is U S V A V' (1/S) U' * IF( ISIM.NE.0 ) THEN * * Compute S (singular values of the eigenvector matrix) * according to CONDS and MODES * CALL DLATM1( MODES, CONDS, 0, 0, ISEED, DS, N, IINFO ) IF( IINFO.NE.0 ) THEN INFO = 3 RETURN END IF * * Multiply by V and V' * CALL ZLARGE( N, A, LDA, ISEED, WORK, IINFO ) IF( IINFO.NE.0 ) THEN INFO = 4 RETURN END IF * * Multiply by S and (1/S) * DO 50 J = 1, N CALL ZDSCAL( N, DS( J ), A( J, 1 ), LDA ) IF( DS( J ).NE.ZERO ) THEN CALL ZDSCAL( N, ONE / DS( J ), A( 1, J ), 1 ) ELSE INFO = 5 RETURN END IF 50 CONTINUE * * Multiply by U and U' * CALL ZLARGE( N, A, LDA, ISEED, WORK, IINFO ) IF( IINFO.NE.0 ) THEN INFO = 4 RETURN END IF END IF * * 5) Reduce the bandwidth. * IF( KL.LT.N-1 ) THEN * * Reduce bandwidth -- kill column * DO 60 JCR = KL + 1, N - 1 IC = JCR - KL IROWS = N + 1 - JCR ICOLS = N + KL - JCR * CALL ZCOPY( IROWS, A( JCR, IC ), 1, WORK, 1 ) XNORMS = WORK( 1 ) CALL ZLARFG( IROWS, XNORMS, WORK( 2 ), 1, TAU ) TAU = DCONJG( TAU ) WORK( 1 ) = CONE ALPHA = ZLARND( 5, ISEED ) * CALL ZGEMV( 'C', IROWS, ICOLS, CONE, A( JCR, IC+1 ), LDA, \$ WORK, 1, CZERO, WORK( IROWS+1 ), 1 ) CALL ZGERC( IROWS, ICOLS, -TAU, WORK, 1, WORK( IROWS+1 ), 1, \$ A( JCR, IC+1 ), LDA ) * CALL ZGEMV( 'N', N, IROWS, CONE, A( 1, JCR ), LDA, WORK, 1, \$ CZERO, WORK( IROWS+1 ), 1 ) CALL ZGERC( N, IROWS, -DCONJG( TAU ), WORK( IROWS+1 ), 1, \$ WORK, 1, A( 1, JCR ), LDA ) * A( JCR, IC ) = XNORMS CALL ZLASET( 'Full', IROWS-1, 1, CZERO, CZERO, \$ A( JCR+1, IC ), LDA ) * CALL ZSCAL( ICOLS+1, ALPHA, A( JCR, IC ), LDA ) CALL ZSCAL( N, DCONJG( ALPHA ), A( 1, JCR ), 1 ) 60 CONTINUE ELSE IF( KU.LT.N-1 ) THEN * * Reduce upper bandwidth -- kill a row at a time. * DO 70 JCR = KU + 1, N - 1 IR = JCR - KU IROWS = N + KU - JCR ICOLS = N + 1 - JCR * CALL ZCOPY( ICOLS, A( IR, JCR ), LDA, WORK, 1 ) XNORMS = WORK( 1 ) CALL ZLARFG( ICOLS, XNORMS, WORK( 2 ), 1, TAU ) TAU = DCONJG( TAU ) WORK( 1 ) = CONE CALL ZLACGV( ICOLS-1, WORK( 2 ), 1 ) ALPHA = ZLARND( 5, ISEED ) * CALL ZGEMV( 'N', IROWS, ICOLS, CONE, A( IR+1, JCR ), LDA, \$ WORK, 1, CZERO, WORK( ICOLS+1 ), 1 ) CALL ZGERC( IROWS, ICOLS, -TAU, WORK( ICOLS+1 ), 1, WORK, 1, \$ A( IR+1, JCR ), LDA ) * CALL ZGEMV( 'C', ICOLS, N, CONE, A( JCR, 1 ), LDA, WORK, 1, \$ CZERO, WORK( ICOLS+1 ), 1 ) CALL ZGERC( ICOLS, N, -DCONJG( TAU ), WORK, 1, \$ WORK( ICOLS+1 ), 1, A( JCR, 1 ), LDA ) * A( IR, JCR ) = XNORMS CALL ZLASET( 'Full', 1, ICOLS-1, CZERO, CZERO, \$ A( IR, JCR+1 ), LDA ) * CALL ZSCAL( IROWS+1, ALPHA, A( IR, JCR ), 1 ) CALL ZSCAL( N, DCONJG( ALPHA ), A( JCR, 1 ), LDA ) 70 CONTINUE END IF * * Scale the matrix to have norm ANORM * IF( ANORM.GE.ZERO ) THEN TEMP = ZLANGE( 'M', N, N, A, LDA, TEMPA ) IF( TEMP.GT.ZERO ) THEN RALPHA = ANORM / TEMP DO 80 J = 1, N CALL ZDSCAL( N, RALPHA, A( 1, J ), 1 ) 80 CONTINUE END IF END IF * RETURN * * End of ZLATME * END