slatme()

subroutine slatme	(	integer	n,
		character	dist,
		integer, dimension( 4 )	iseed,
		real, dimension( * )	d,
		integer	mode,
		real	cond,
		real	dmax,
		character, dimension( * )	ei,
		character	rsign,
		character	upper,
		character	sim,
		real, dimension( * )	ds,
		integer	modes,
		real	conds,
		integer	kl,
		integer	ku,
		real	anorm,
		real, dimension( lda, * )	a,
		integer	lda,
		real, dimension( * )	work,
		integer	info )

SLATME

Purpose:

!>
!>    SLATME generates random non-symmetric square matrices with
!>    specified eigenvalues for testing LAPACK programs.
!>
!>    SLATME 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 complex conjugate pairs are desired (MODE=0 and EI(1)='R',
!>         or MODE=5), certain pairs of adjacent elements of D are
!>         interpreted as the real and complex parts of a complex
!>         conjugate pair; A thus becomes block diagonal, with 1x1
!>         and 2x2 blocks.
!>
!>    3. If UPPER='T', the upper triangle of A is set to random values
!>         out of distribution DIST.
!>
!>    4. 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.
!>
!>    5. 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.
!>
!>    6. 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.)
!> 

Parameters

[in]	N	!> N is INTEGER !> The number of columns (or rows) of A. Not modified. !>
[in]	DIST	!> DIST is CHARACTER*1 !> On entry, DIST specifies the type of distribution to be used !> to generate the random eigen-/singular values, and for 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 ) !> Not modified. !>
[in,out]	ISEED	!> ISEED is 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 SLATME !> to continue the same random number sequence. !> Changed on exit. !>
[in,out]	D	!> D is REAL array, dimension ( N ) !> This array is used to specify the eigenvalues of A. If !> MODE=0, then D is assumed to contain the eigenvalues (but !> see the description of EI), otherwise they will be !> computed according to MODE, COND, DMAX, and RSIGN and !> placed in D. !> Modified if MODE is nonzero. !>
[in]	MODE	!> MODE is INTEGER !> On entry this describes how the eigenvalues are to !> be specified: !> MODE = 0 means use D (with EI) 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. Each odd-even pair !> of elements will be either used as two real !> eigenvalues or as the real and imaginary part !> of a complex conjugate pair of eigenvalues; !> the choice of which is done is random, with !> 50-50 probability, for each pair. !> 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. !>
[in]	COND	!> COND is REAL !> On entry, this is used as described under MODE above. !> If used, it must be >= 1. Not modified. !>
[in]	DMAX	!> DMAX is REAL !> 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: if DMAX is negative (or zero), D will be !> scaled by a negative number (or zero). !> Not modified. !>
[in]	EI	!> EI is CHARACTER*1 array, dimension ( N ) !> If MODE is 0, and EI(1) is not ' ' (space character), !> this array specifies which elements of D (on input) are !> real eigenvalues and which are the real and imaginary parts !> of a complex conjugate pair of eigenvalues. The elements !> of EI may then only have the values 'R' and 'I'. If !> EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is !> CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex !> conjugate thereof. If EI(j)=EI(j+1)='R', then the j-th !> eigenvalue is D(j) (i.e., real). EI(1) may not be 'I', !> nor may two adjacent elements of EI both have the value 'I'. !> If MODE is not 0, then EI is ignored. If MODE is 0 and !> EI(1)=' ', then the eigenvalues will all be real. !> Not modified. !>
[in]	RSIGN	!> RSIGN is 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 sign (+1 or -1). If RSIGN='F', !> they will not be. RSIGN may only have the values 'T' or !> 'F'. !> Not modified. !>
[in]	UPPER	!> UPPER is CHARACTER*1 !> If UPPER='T', then the elements of A above the diagonal !> (and above the 2x2 diagonal blocks, if A has complex !> eigenvalues) 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. !>
[in]	SIM	!> SIM is CHARACTER*1 !> If SIM='T', then A will be operated on by a , 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. !>
[in,out]	DS	!> DS is REAL 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. !>
[in]	MODES	!> MODES is INTEGER !>
[in]	CONDS	!> CONDS is REAL !> Same as 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.) !>
[in]	KL	!> KL is 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. KL must be at !> least 1. !> Not modified. !>
[in]	KU	!> KU is 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. KU must be at least 1. !> Not modified. !>
[in]	ANORM	!> ANORM is REAL !> 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. !>
[out]	A	!> A is REAL array, dimension ( LDA, N ) !> On exit A is the desired test matrix. !> Modified. !>
[in]	LDA	!> LDA is INTEGER !> LDA specifies the first dimension of A as declared in the !> calling program. LDA must be at least N. !> Not modified. !>
[out]	WORK	!> WORK is REAL array, dimension ( 3*N ) !> Workspace. !> Modified. !>
[out]	INFO	!> INFO is 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 !> -8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or !> two adjacent elements of EI are 'I'. !> -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 N. !> 1 => Error return from SLATM1 (computing D) !> 2 => Cannot scale to DMAX (max. eigenvalue is 0) !> 3 => Error return from SLATM1 (computing DS) !> 4 => Error return from SLARGE !> 5 => Zero singular value from SLATM1. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 327 of file slatme.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          DIST, RSIGN, SIM, UPPER
      INTEGER            INFO, KL, KU, LDA, MODE, MODES, N
      REAL               ANORM, COND, CONDS, DMAX
*     ..
*     .. Array Arguments ..
      CHARACTER          EI( * )
      INTEGER            ISEED( 4 )
      REAL               A( LDA, * ), D( * ), DS( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      parameter( zero = 0.0e0 )
      REAL               ONE
      parameter( one = 1.0e0 )
      REAL               HALF
      parameter( half = 1.0e0 / 2.0e0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            BADEI, BADS, USEEI
      INTEGER            I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN,
     $                   ISIM, IUPPER, J, JC, JCR, JR
      REAL               ALPHA, TAU, TEMP, XNORMS
*     ..
*     .. Local Arrays ..
      REAL               TEMPA( 1 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLANGE, SLARAN
      EXTERNAL           lsame, slange, slaran
*     ..
*     .. External Subroutines ..
      EXTERNAL           scopy, sgemv, sger, slarfg, slarge, slarnv,
     $                   slatm1, slaset, sscal, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, 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
         idist = -1
      END IF
*
*     Check EI
*
      useei = .true.
      badei = .false.
      IF( lsame( ei( 1 ), ' ' ) .OR. mode.NE.0 ) THEN
         useei = .false.
      ELSE
         IF( lsame( ei( 1 ), 'R' ) ) THEN
            DO 10 j = 2, n
               IF( lsame( ei( j ), 'I' ) ) THEN
                  IF( lsame( ei( j-1 ), 'I' ) )
     $               badei = .true.
               ELSE
                  IF( .NOT.lsame( ei( j ), 'R' ) )
     $               badei = .true.
               END IF
   10       CONTINUE
         ELSE
            badei = .true.
         END IF
      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 20 j = 1, n
            IF( ds( j ).EQ.zero )
     $         bads = .true.
   20    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( badei ) THEN
         info = -8
      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( 'SLATME', -info )
         RETURN
      END IF
*
*     Initialize random number generator
*
      DO 30 i = 1, 4
         iseed( i ) = mod( abs( iseed( i ) ), 4096 )
   30 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 slatm1( 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 40 i = 2, n
            temp = max( temp, abs( d( i ) ) )
   40    CONTINUE
*
         IF( temp.GT.zero ) THEN
            alpha = dmax / temp
         ELSE IF( dmax.NE.zero ) THEN
            info = 2
            RETURN
         ELSE
            alpha = zero
         END IF
*
         CALL sscal( n, alpha, d, 1 )
*
      END IF
*
      CALL slaset( 'Full', n, n, zero, zero, a, lda )
      CALL scopy( n, d, 1, a, lda+1 )
*
*     Set up complex conjugate pairs
*
      IF( mode.EQ.0 ) THEN
         IF( useei ) THEN
            DO 50 j = 2, n
               IF( lsame( ei( j ), 'I' ) ) THEN
                  a( j-1, j ) = a( j, j )
                  a( j, j-1 ) = -a( j, j )
                  a( j, j ) = a( j-1, j-1 )
               END IF
   50       CONTINUE
         END IF
*
      ELSE IF( abs( mode ).EQ.5 ) THEN
*
         DO 60 j = 2, n, 2
            IF( slaran( iseed ).GT.half ) THEN
               a( j-1, j ) = a( j, j )
               a( j, j-1 ) = -a( j, j )
               a( j, j ) = a( j-1, j-1 )
            END IF
   60    CONTINUE
      END IF
*
*     3)      If UPPER='T', set upper triangle of A to random numbers.
*             (but don't modify the corners of 2x2 blocks.)
*
      IF( iupper.NE.0 ) THEN
         DO 70 jc = 2, n
            IF( a( jc-1, jc ).NE.zero ) THEN
               jr = jc - 2
            ELSE
               jr = jc - 1
            END IF
            CALL slarnv( idist, iseed, jr, a( 1, jc ) )
   70    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 slatm1( modes, conds, 0, 0, iseed, ds, n, iinfo )
         IF( iinfo.NE.0 ) THEN
            info = 3
            RETURN
         END IF
*
*        Multiply by V and V'
*
         CALL slarge( n, a, lda, iseed, work, iinfo )
         IF( iinfo.NE.0 ) THEN
            info = 4
            RETURN
         END IF
*
*        Multiply by S and (1/S)
*
         DO 80 j = 1, n
            CALL sscal( n, ds( j ), a( j, 1 ), lda )
            IF( ds( j ).NE.zero ) THEN
               CALL sscal( n, one / ds( j ), a( 1, j ), 1 )
            ELSE
               info = 5
               RETURN
            END IF
   80    CONTINUE
*
*        Multiply by U and U'
*
         CALL slarge( 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 90 jcr = kl + 1, n - 1
            ic = jcr - kl
            irows = n + 1 - jcr
            icols = n + kl - jcr
*
            CALL scopy( irows, a( jcr, ic ), 1, work, 1 )
            xnorms = work( 1 )
            CALL slarfg( irows, xnorms, work( 2 ), 1, tau )
            work( 1 ) = one
*
            CALL sgemv( 'T', irows, icols, one, a( jcr, ic+1 ), lda,
     $                  work, 1, zero, work( irows+1 ), 1 )
            CALL sger( irows, icols, -tau, work, 1, work( irows+1 ), 1,
     $                 a( jcr, ic+1 ), lda )
*
            CALL sgemv( 'N', n, irows, one, a( 1, jcr ), lda, work, 1,
     $                  zero, work( irows+1 ), 1 )
            CALL sger( n, irows, -tau, work( irows+1 ), 1, work, 1,
     $                 a( 1, jcr ), lda )
*
            a( jcr, ic ) = xnorms
            CALL slaset( 'Full', irows-1, 1, zero, zero, a( jcr+1, ic ),
     $                   lda )
   90    CONTINUE
      ELSE IF( ku.LT.n-1 ) THEN
*
*        Reduce upper bandwidth -- kill a row at a time.
*
         DO 100 jcr = ku + 1, n - 1
            ir = jcr - ku
            irows = n + ku - jcr
            icols = n + 1 - jcr
*
            CALL scopy( icols, a( ir, jcr ), lda, work, 1 )
            xnorms = work( 1 )
            CALL slarfg( icols, xnorms, work( 2 ), 1, tau )
            work( 1 ) = one
*
            CALL sgemv( 'N', irows, icols, one, a( ir+1, jcr ), lda,
     $                  work, 1, zero, work( icols+1 ), 1 )
            CALL sger( irows, icols, -tau, work( icols+1 ), 1, work, 1,
     $                 a( ir+1, jcr ), lda )
*
            CALL sgemv( 'C', icols, n, one, a( jcr, 1 ), lda, work, 1,
     $                  zero, work( icols+1 ), 1 )
            CALL sger( icols, n, -tau, work, 1, work( icols+1 ), 1,
     $                 a( jcr, 1 ), lda )
*
            a( ir, jcr ) = xnorms
            CALL slaset( 'Full', 1, icols-1, zero, zero, a( ir, jcr+1 ),
     $                   lda )
  100    CONTINUE
      END IF
*
*     Scale the matrix to have norm ANORM
*
      IF( anorm.GE.zero ) THEN
         temp = slange( 'M', n, n, a, lda, tempa )
         IF( temp.GT.zero ) THEN
            alpha = anorm / temp
            DO 110 j = 1, n
               CALL sscal( n, alpha, a( 1, j ), 1 )
  110       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     End of SLATME
*

Here is the call graph for this function:

Here is the caller graph for this function: