clatsy.f

Go to the documentation of this file.

*> \brief \b CLATSY
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE CLATSY( UPLO, N, X, LDX, ISEED )
* 
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            LDX, N
*       ..
*       .. Array Arguments ..
*       INTEGER            ISEED( * )
*       COMPLEX            X( LDX, * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CLATSY generates a special test matrix for the complex symmetric
*> (indefinite) factorization.  The pivot blocks of the generated matrix
*> will be in the following order:
*>    2x2 pivot block, non diagonalizable
*>    1x1 pivot block
*>    2x2 pivot block, diagonalizable
*>    (cycle repeats)
*> A row interchange is required for each non-diagonalizable 2x2 block.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER
*>          Specifies whether the generated matrix is to be upper or
*>          lower triangular.
*>          = 'U':  Upper triangular
*>          = 'L':  Lower triangular
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The dimension of the matrix to be generated.
*> \endverbatim
*>
*> \param[out] X
*> \verbatim
*>          X is COMPLEX array, dimension (LDX,N)
*>          The generated matrix, consisting of 3x3 and 2x2 diagonal
*>          blocks which result in the pivot sequence given above.
*>          The matrix outside of these diagonal blocks is zero.
*> \endverbatim
*>
*> \param[in] LDX
*> \verbatim
*>          LDX is INTEGER
*>          The leading dimension of the array X.
*> \endverbatim
*>
*> \param[in,out] ISEED
*> \verbatim
*>          ISEED is INTEGER array, dimension (4)
*>          On entry, the seed for the random number generator.  The last
*>          of the four integers must be odd.  (modified on exit)
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup complex_lin
*
*  =====================================================================
      SUBROUTINE clatsy( UPLO, N, X, LDX, ISEED )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            LDX, N
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( * )
      COMPLEX            X( ldx, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX            EYE
      parameter                ( eye = ( 0.0, 1.0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J, N5
      REAL               ALPHA, ALPHA3, BETA
      COMPLEX            A, B, C, R
*     ..
*     .. External Functions ..
      COMPLEX            CLARND
      EXTERNAL           clarnd
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, sqrt
*     ..
*     .. Executable Statements ..
*
*     Initialize constants
*
      alpha = ( 1.+sqrt( 17. ) ) / 8.
      beta = alpha - 1. / 1000.
      alpha3 = alpha*alpha*alpha
*
*     UPLO = 'U':  Upper triangular storage
*
      IF( uplo.EQ.'U' ) THEN
*
*        Fill the upper triangle of the matrix with zeros.
*
         DO 20 j = 1, n
            DO 10 i = 1, j
               x( i, j ) = 0.0
   10       CONTINUE
   20    CONTINUE
         n5 = n / 5
         n5 = n - 5*n5 + 1
*
         DO 30 i = n, n5, -5
            a = alpha3*clarnd( 5, iseed )
            b = clarnd( 5, iseed ) / alpha
            c = a - 2.*b*eye
            r = c / beta
            x( i, i ) = a
            x( i-2, i ) = b
            x( i-2, i-1 ) = r
            x( i-2, i-2 ) = c
            x( i-1, i-1 ) = clarnd( 2, iseed )
            x( i-3, i-3 ) = clarnd( 2, iseed )
            x( i-4, i-4 ) = clarnd( 2, iseed )
            IF( abs( x( i-3, i-3 ) ).GT.abs( x( i-4, i-4 ) ) ) THEN
               x( i-4, i-3 ) = 2.0*x( i-3, i-3 )
            ELSE
               x( i-4, i-3 ) = 2.0*x( i-4, i-4 )
            END IF
   30    CONTINUE
*
*        Clean-up for N not a multiple of 5.
*
         i = n5 - 1
         IF( i.GT.2 ) THEN
            a = alpha3*clarnd( 5, iseed )
            b = clarnd( 5, iseed ) / alpha
            c = a - 2.*b*eye
            r = c / beta
            x( i, i ) = a
            x( i-2, i ) = b
            x( i-2, i-1 ) = r
            x( i-2, i-2 ) = c
            x( i-1, i-1 ) = clarnd( 2, iseed )
            i = i - 3
         END IF
         IF( i.GT.1 ) THEN
            x( i, i ) = clarnd( 2, iseed )
            x( i-1, i-1 ) = clarnd( 2, iseed )
            IF( abs( x( i, i ) ).GT.abs( x( i-1, i-1 ) ) ) THEN
               x( i-1, i ) = 2.0*x( i, i )
            ELSE
               x( i-1, i ) = 2.0*x( i-1, i-1 )
            END IF
            i = i - 2
         ELSE IF( i.EQ.1 ) THEN
            x( i, i ) = clarnd( 2, iseed )
            i = i - 1
         END IF
*
*     UPLO = 'L':  Lower triangular storage
*
      ELSE
*
*        Fill the lower triangle of the matrix with zeros.
*
         DO 50 j = 1, n
            DO 40 i = j, n
               x( i, j ) = 0.0
   40       CONTINUE
   50    CONTINUE
         n5 = n / 5
         n5 = n5*5
*
         DO 60 i = 1, n5, 5
            a = alpha3*clarnd( 5, iseed )
            b = clarnd( 5, iseed ) / alpha
            c = a - 2.*b*eye
            r = c / beta
            x( i, i ) = a
            x( i+2, i ) = b
            x( i+2, i+1 ) = r
            x( i+2, i+2 ) = c
            x( i+1, i+1 ) = clarnd( 2, iseed )
            x( i+3, i+3 ) = clarnd( 2, iseed )
            x( i+4, i+4 ) = clarnd( 2, iseed )
            IF( abs( x( i+3, i+3 ) ).GT.abs( x( i+4, i+4 ) ) ) THEN
               x( i+4, i+3 ) = 2.0*x( i+3, i+3 )
            ELSE
               x( i+4, i+3 ) = 2.0*x( i+4, i+4 )
            END IF
   60    CONTINUE
*
*        Clean-up for N not a multiple of 5.
*
         i = n5 + 1
         IF( i.LT.n-1 ) THEN
            a = alpha3*clarnd( 5, iseed )
            b = clarnd( 5, iseed ) / alpha
            c = a - 2.*b*eye
            r = c / beta
            x( i, i ) = a
            x( i+2, i ) = b
            x( i+2, i+1 ) = r
            x( i+2, i+2 ) = c
            x( i+1, i+1 ) = clarnd( 2, iseed )
            i = i + 3
         END IF
         IF( i.LT.n ) THEN
            x( i, i ) = clarnd( 2, iseed )
            x( i+1, i+1 ) = clarnd( 2, iseed )
            IF( abs( x( i, i ) ).GT.abs( x( i+1, i+1 ) ) ) THEN
               x( i+1, i ) = 2.0*x( i, i )
            ELSE
               x( i+1, i ) = 2.0*x( i+1, i+1 )
            END IF
            i = i + 2
         ELSE IF( i.EQ.n ) THEN
            x( i, i ) = clarnd( 2, iseed )
            i = i + 1
         END IF
      END IF
*
      RETURN
*
*     End of CLATSY
*
      END