path: root/TESTING/MATGEN/zlatme.f

      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        (input) INTEGER
*           The number of columns (or rows) of A. Not modified.
*
*  DIST     (input) 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    (input/output) 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        (input/output) 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     (input) 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     (input) DOUBLE PRECISION
*           On entry, this is used as described under MODE above.
*           If used, it must be >= 1. Not modified.
*
*  DMAX     (input) 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       (input) CHARACTER*1 array, dimension ( N )
*           (ignored)
*           Not modified.
*
*  RSIGN    (input) 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    (input) 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      (input) 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       (input/output) 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    (input) INTEGER
*
*  CONDS    (input) 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       (input) 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       (input) 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    (input) 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        (output) COMPLEX*16 array, dimension ( LDA, N )
*           On exit A is the desired test matrix.
*           Modified.
*
*  LDA      (input) INTEGER
*           LDA specifies the first dimension of A as declared in the
*           calling program.  LDA must be at least M.
*           Not modified.
*
*  WORK     (workspace) COMPLEX*16 array, dimension ( 3*N )
*           Workspace.
*           Modified.
*
*  INFO     (output) 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