dlatm4.f

Go to the documentation of this file.

*> \brief \b DLATM4
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLATM4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND,
*                          TRIANG, IDIST, ISEED, A, LDA )
*
*       .. Scalar Arguments ..
*       INTEGER            IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2
*       DOUBLE PRECISION   AMAGN, RCOND, TRIANG
*       ..
*       .. Array Arguments ..
*       INTEGER            ISEED( 4 )
*       DOUBLE PRECISION   A( LDA, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLATM4 generates basic square matrices, which may later be
*> multiplied by others in order to produce test matrices.  It is
*> intended mainly to be used to test the generalized eigenvalue
*> routines.
*>
*> It first generates the diagonal and (possibly) subdiagonal,
*> according to the value of ITYPE, NZ1, NZ2, ISIGN, AMAGN, and RCOND.
*> It then fills in the upper triangle with random numbers, if TRIANG is
*> non-zero.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] ITYPE
*> \verbatim
*>          ITYPE is INTEGER
*>          The "type" of matrix on the diagonal and sub-diagonal.
*>          If ITYPE < 0, then type abs(ITYPE) is generated and then
*>             swapped end for end (A(I,J) := A'(N-J,N-I).)  See also
*>             the description of AMAGN and ISIGN.
*>
*>          Special types:
*>          = 0:  the zero matrix.
*>          = 1:  the identity.
*>          = 2:  a transposed Jordan block.
*>          = 3:  If N is odd, then a k+1 x k+1 transposed Jordan block
*>                followed by a k x k identity block, where k=(N-1)/2.
*>                If N is even, then k=(N-2)/2, and a zero diagonal entry
*>                is tacked onto the end.
*>
*>          Diagonal types.  The diagonal consists of NZ1 zeros, then
*>             k=N-NZ1-NZ2 nonzeros.  The subdiagonal is zero.  ITYPE
*>             specifies the nonzero diagonal entries as follows:
*>          = 4:  1, ..., k
*>          = 5:  1, RCOND, ..., RCOND
*>          = 6:  1, ..., 1, RCOND
*>          = 7:  1, a, a^2, ..., a^(k-1)=RCOND
*>          = 8:  1, 1-d, 1-2*d, ..., 1-(k-1)*d=RCOND
*>          = 9:  random numbers chosen from (RCOND,1)
*>          = 10: random numbers with distribution IDIST (see DLARND.)
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix.
*> \endverbatim
*>
*> \param[in] NZ1
*> \verbatim
*>          NZ1 is INTEGER
*>          If abs(ITYPE) > 3, then the first NZ1 diagonal entries will
*>          be zero.
*> \endverbatim
*>
*> \param[in] NZ2
*> \verbatim
*>          NZ2 is INTEGER
*>          If abs(ITYPE) > 3, then the last NZ2 diagonal entries will
*>          be zero.
*> \endverbatim
*>
*> \param[in] ISIGN
*> \verbatim
*>          ISIGN is INTEGER
*>          = 0: The sign of the diagonal and subdiagonal entries will
*>               be left unchanged.
*>          = 1: The diagonal and subdiagonal entries will have their
*>               sign changed at random.
*>          = 2: If ITYPE is 2 or 3, then the same as ISIGN=1.
*>               Otherwise, with probability 0.5, odd-even pairs of
*>               diagonal entries A(2*j-1,2*j-1), A(2*j,2*j) will be
*>               converted to a 2x2 block by pre- and post-multiplying
*>               by distinct random orthogonal rotations.  The remaining
*>               diagonal entries will have their sign changed at random.
*> \endverbatim
*>
*> \param[in] AMAGN
*> \verbatim
*>          AMAGN is DOUBLE PRECISION
*>          The diagonal and subdiagonal entries will be multiplied by
*>          AMAGN.
*> \endverbatim
*>
*> \param[in] RCOND
*> \verbatim
*>          RCOND is DOUBLE PRECISION
*>          If abs(ITYPE) > 4, then the smallest diagonal entry will be
*>          entry will be RCOND.  RCOND must be between 0 and 1.
*> \endverbatim
*>
*> \param[in] TRIANG
*> \verbatim
*>          TRIANG is DOUBLE PRECISION
*>          The entries above the diagonal will be random numbers with
*>          magnitude bounded by TRIANG (i.e., random numbers multiplied
*>          by TRIANG.)
*> \endverbatim
*>
*> \param[in] IDIST
*> \verbatim
*>          IDIST is INTEGER
*>          Specifies the type of distribution to be used to generate a
*>          random matrix.
*>          = 1:  UNIFORM( 0, 1 )
*>          = 2:  UNIFORM( -1, 1 )
*>          = 3:  NORMAL ( 0, 1 )
*> \endverbatim
*>
*> \param[in,out] ISEED
*> \verbatim
*>          ISEED is INTEGER array, dimension (4)
*>          On entry ISEED specifies the seed of the random number
*>          generator.  The values of ISEED are changed on exit, and can
*>          be used in the next call to DLATM4 to continue the same
*>          random number sequence.
*>          Note: ISEED(4) should be odd, for the random number generator
*>          used at present.
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, N)
*>          Array to be computed.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          Leading dimension of A.  Must be at least 1 and at least N.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE dlatm4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND,
     $                   TRIANG, IDIST, ISEED, A, LDA )
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2
      DOUBLE PRECISION   AMAGN, RCOND, TRIANG
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( 4 )
      DOUBLE PRECISION   A( LDA, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO
      parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
      DOUBLE PRECISION   HALF
      parameter( half = one / two )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IOFF, ISDB, ISDE, JC, JD, JR, K, KBEG, KEND,
     $                   klen
      DOUBLE PRECISION   ALPHA, CL, CR, SAFMIN, SL, SR, SV1, SV2, TEMP
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLARAN, DLARND
      EXTERNAL           dlamch, dlaran, dlarnd
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlaset
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, exp, log, max, min, mod, sqrt
*     ..
*     .. Executable Statements ..
*
      IF( n.LE.0 )
     $   RETURN
      CALL dlaset( 'Full', n, n, zero, zero, a, lda )
*
*     Insure a correct ISEED
*
      IF( mod( iseed( 4 ), 2 ).NE.1 )
     $   iseed( 4 ) = iseed( 4 ) + 1
*
*     Compute diagonal and subdiagonal according to ITYPE, NZ1, NZ2,
*     and RCOND
*
      IF( itype.NE.0 ) THEN
         IF( abs( itype ).GE.4 ) THEN
            kbeg = max( 1, min( n, nz1+1 ) )
            kend = max( kbeg, min( n, n-nz2 ) )
            klen = kend + 1 - kbeg
         ELSE
            kbeg = 1
            kend = n
            klen = n
         END IF
         isdb = 1
         isde = 0
         GO TO ( 10, 30, 50, 80, 100, 120, 140, 160,
     $           180, 200 )abs( itype )
*
*        abs(ITYPE) = 1: Identity
*
   10    CONTINUE
         DO 20 jd = 1, n
            a( jd, jd ) = one
   20    CONTINUE
         GO TO 220
*
*        abs(ITYPE) = 2: Transposed Jordan block
*
   30    CONTINUE
         DO 40 jd = 1, n - 1
            a( jd+1, jd ) = one
   40    CONTINUE
         isdb = 1
         isde = n - 1
         GO TO 220
*
*        abs(ITYPE) = 3: Transposed Jordan block, followed by the
*                        identity.
*
   50    CONTINUE
         k = ( n-1 ) / 2
         DO 60 jd = 1, k
            a( jd+1, jd ) = one
   60    CONTINUE
         isdb = 1
         isde = k
         DO 70 jd = k + 2, 2*k + 1
            a( jd, jd ) = one
   70    CONTINUE
         GO TO 220
*
*        abs(ITYPE) = 4: 1,...,k
*
   80    CONTINUE
         DO 90 jd = kbeg, kend
            a( jd, jd ) = dble( jd-nz1 )
   90    CONTINUE
         GO TO 220
*
*        abs(ITYPE) = 5: One large D value:
*
  100    CONTINUE
         DO 110 jd = kbeg + 1, kend
            a( jd, jd ) = rcond
  110    CONTINUE
         a( kbeg, kbeg ) = one
         GO TO 220
*
*        abs(ITYPE) = 6: One small D value:
*
  120    CONTINUE
         DO 130 jd = kbeg, kend - 1
            a( jd, jd ) = one
  130    CONTINUE
         a( kend, kend ) = rcond
         GO TO 220
*
*        abs(ITYPE) = 7: Exponentially distributed D values:
*
  140    CONTINUE
         a( kbeg, kbeg ) = one
         IF( klen.GT.1 ) THEN
            alpha = rcond**( one / dble( klen-1 ) )
            DO 150 i = 2, klen
               a( nz1+i, nz1+i ) = alpha**dble( i-1 )
  150       CONTINUE
         END IF
         GO TO 220
*
*        abs(ITYPE) = 8: Arithmetically distributed D values:
*
  160    CONTINUE
         a( kbeg, kbeg ) = one
         IF( klen.GT.1 ) THEN
            alpha = ( one-rcond ) / dble( klen-1 )
            DO 170 i = 2, klen
               a( nz1+i, nz1+i ) = dble( klen-i )*alpha + rcond
  170       CONTINUE
         END IF
         GO TO 220
*
*        abs(ITYPE) = 9: Randomly distributed D values on ( RCOND, 1):
*
  180    CONTINUE
         alpha = log( rcond )
         DO 190 jd = kbeg, kend
            a( jd, jd ) = exp( alpha*dlaran( iseed ) )
  190    CONTINUE
         GO TO 220
*
*        abs(ITYPE) = 10: Randomly distributed D values from DIST
*
  200    CONTINUE
         DO 210 jd = kbeg, kend
            a( jd, jd ) = dlarnd( idist, iseed )
  210    CONTINUE
*
  220    CONTINUE
*
*        Scale by AMAGN
*
         DO 230 jd = kbeg, kend
            a( jd, jd ) = amagn*dble( a( jd, jd ) )
  230    CONTINUE
         DO 240 jd = isdb, isde
            a( jd+1, jd ) = amagn*dble( a( jd+1, jd ) )
  240    CONTINUE
*
*        If ISIGN = 1 or 2, assign random signs to diagonal and
*        subdiagonal
*
         IF( isign.GT.0 ) THEN
            DO 250 jd = kbeg, kend
               IF( dble( a( jd, jd ) ).NE.zero ) THEN
                  IF( dlaran( iseed ).GT.half )
     $               a( jd, jd ) = -a( jd, jd )
               END IF
  250       CONTINUE
            DO 260 jd = isdb, isde
               IF( dble( a( jd+1, jd ) ).NE.zero ) THEN
                  IF( dlaran( iseed ).GT.half )
     $               a( jd+1, jd ) = -a( jd+1, jd )
               END IF
  260       CONTINUE
         END IF
*
*        Reverse if ITYPE < 0
*
         IF( itype.LT.0 ) THEN
            DO 270 jd = kbeg, ( kbeg+kend-1 ) / 2
               temp = a( jd, jd )
               a( jd, jd ) = a( kbeg+kend-jd, kbeg+kend-jd )
               a( kbeg+kend-jd, kbeg+kend-jd ) = temp
  270       CONTINUE
            DO 280 jd = 1, ( n-1 ) / 2
               temp = a( jd+1, jd )
               a( jd+1, jd ) = a( n+1-jd, n-jd )
               a( n+1-jd, n-jd ) = temp
  280       CONTINUE
         END IF
*
*        If ISIGN = 2, and no subdiagonals already, then apply
*        random rotations to make 2x2 blocks.
*
         IF( isign.EQ.2 .AND. itype.NE.2 .AND. itype.NE.3 ) THEN
            safmin = dlamch( 'S' )
            DO 290 jd = kbeg, kend - 1, 2
               IF( dlaran( iseed ).GT.half ) THEN
*
*                 Rotation on left.
*
                  cl = two*dlaran( iseed ) - one
                  sl = two*dlaran( iseed ) - one
                  temp = one / max( safmin, sqrt( cl**2+sl**2 ) )
                  cl = cl*temp
                  sl = sl*temp
*
*                 Rotation on right.
*
                  cr = two*dlaran( iseed ) - one
                  sr = two*dlaran( iseed ) - one
                  temp = one / max( safmin, sqrt( cr**2+sr**2 ) )
                  cr = cr*temp
                  sr = sr*temp
*
*                 Apply
*
                  sv1 = a( jd, jd )
                  sv2 = a( jd+1, jd+1 )
                  a( jd, jd ) = cl*cr*sv1 + sl*sr*sv2
                  a( jd+1, jd ) = -sl*cr*sv1 + cl*sr*sv2
                  a( jd, jd+1 ) = -cl*sr*sv1 + sl*cr*sv2
                  a( jd+1, jd+1 ) = sl*sr*sv1 + cl*cr*sv2
               END IF
  290       CONTINUE
         END IF
*
      END IF
*
*     Fill in upper triangle (except for 2x2 blocks)
*
      IF( triang.NE.zero ) THEN
         IF( isign.NE.2 .OR. itype.EQ.2 .OR. itype.EQ.3 ) THEN
            ioff = 1
         ELSE
            ioff = 2
            DO 300 jr = 1, n - 1
               IF( a( jr+1, jr ).EQ.zero )
     $            a( jr, jr+1 ) = triang*dlarnd( idist, iseed )
  300       CONTINUE
         END IF
*
         DO 320 jc = 2, n
            DO 310 jr = 1, jc - ioff
               a( jr, jc ) = triang*dlarnd( idist, iseed )
  310       CONTINUE
  320    CONTINUE
      END IF
*
      RETURN
*
*     End of DLATM4
*
      END