clatsp.f

Go to the documentation of this file.

*> \brief \b CLATSP
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE CLATSP( UPLO, N, X, ISEED )
* 
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            N
*       ..
*       .. Array Arguments ..
*       INTEGER            ISEED( * )
*       COMPLEX            X( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CLATSP generates a special test matrix for the complex symmetric
*> (indefinite) factorization for packed matrices.  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 (N*(N+1)/2)
*>          The generated matrix in packed storage format.  The matrix
*>          consists of 3x3 and 2x2 diagonal blocks which result in the
*>          pivot sequence given above.  The matrix outside these
*>          diagonal blocks is zero.
*> \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 clatsp( UPLO, N, X, 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            N
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( * )
      COMPLEX            X( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX            EYE
      parameter                ( eye = ( 0.0, 1.0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            J, JJ, 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
*
*     Fill the matrix with zeros.
*
      DO 10 j = 1, n*( n+1 ) / 2
         x( j ) = 0.0
   10 CONTINUE
*
*     UPLO = 'U':  Upper triangular storage
*
      IF( uplo.EQ.'U' ) THEN
         n5 = n / 5
         n5 = n - 5*n5 + 1
*
         jj = n*( n+1 ) / 2
         DO 20 j = n, n5, -5
            a = alpha3*clarnd( 5, iseed )
            b = clarnd( 5, iseed ) / alpha
            c = a - 2.*b*eye
            r = c / beta
            x( jj ) = a
            x( jj-2 ) = b
            jj = jj - j
            x( jj ) = clarnd( 2, iseed )
            x( jj-1 ) = r
            jj = jj - ( j-1 )
            x( jj ) = c
            jj = jj - ( j-2 )
            x( jj ) = clarnd( 2, iseed )
            jj = jj - ( j-3 )
            x( jj ) = clarnd( 2, iseed )
            IF( abs( x( jj+( j-3 ) ) ).GT.abs( x( jj ) ) ) THEN
               x( jj+( j-4 ) ) = 2.0*x( jj+( j-3 ) )
            ELSE
               x( jj+( j-4 ) ) = 2.0*x( jj )
            END IF
            jj = jj - ( j-4 )
   20    CONTINUE
*
*        Clean-up for N not a multiple of 5.
*
         j = n5 - 1
         IF( j.GT.2 ) THEN
            a = alpha3*clarnd( 5, iseed )
            b = clarnd( 5, iseed ) / alpha
            c = a - 2.*b*eye
            r = c / beta
            x( jj ) = a
            x( jj-2 ) = b
            jj = jj - j
            x( jj ) = clarnd( 2, iseed )
            x( jj-1 ) = r
            jj = jj - ( j-1 )
            x( jj ) = c
            jj = jj - ( j-2 )
            j = j - 3
         END IF
         IF( j.GT.1 ) THEN
            x( jj ) = clarnd( 2, iseed )
            x( jj-j ) = clarnd( 2, iseed )
            IF( abs( x( jj ) ).GT.abs( x( jj-j ) ) ) THEN
               x( jj-1 ) = 2.0*x( jj )
            ELSE
               x( jj-1 ) = 2.0*x( jj-j )
            END IF
            jj = jj - j - ( j-1 )
            j = j - 2
         ELSE IF( j.EQ.1 ) THEN
            x( jj ) = clarnd( 2, iseed )
            j = j - 1
         END IF
*
*     UPLO = 'L':  Lower triangular storage
*
      ELSE
         n5 = n / 5
         n5 = n5*5
*
         jj = 1
         DO 30 j = 1, n5, 5
            a = alpha3*clarnd( 5, iseed )
            b = clarnd( 5, iseed ) / alpha
            c = a - 2.*b*eye
            r = c / beta
            x( jj ) = a
            x( jj+2 ) = b
            jj = jj + ( n-j+1 )
            x( jj ) = clarnd( 2, iseed )
            x( jj+1 ) = r
            jj = jj + ( n-j )
            x( jj ) = c
            jj = jj + ( n-j-1 )
            x( jj ) = clarnd( 2, iseed )
            jj = jj + ( n-j-2 )
            x( jj ) = clarnd( 2, iseed )
            IF( abs( x( jj-( n-j-2 ) ) ).GT.abs( x( jj ) ) ) THEN
               x( jj-( n-j-2 )+1 ) = 2.0*x( jj-( n-j-2 ) )
            ELSE
               x( jj-( n-j-2 )+1 ) = 2.0*x( jj )
            END IF
            jj = jj + ( n-j-3 )
   30    CONTINUE
*
*        Clean-up for N not a multiple of 5.
*
         j = n5 + 1
         IF( j.LT.n-1 ) THEN
            a = alpha3*clarnd( 5, iseed )
            b = clarnd( 5, iseed ) / alpha
            c = a - 2.*b*eye
            r = c / beta
            x( jj ) = a
            x( jj+2 ) = b
            jj = jj + ( n-j+1 )
            x( jj ) = clarnd( 2, iseed )
            x( jj+1 ) = r
            jj = jj + ( n-j )
            x( jj ) = c
            jj = jj + ( n-j-1 )
            j = j + 3
         END IF
         IF( j.LT.n ) THEN
            x( jj ) = clarnd( 2, iseed )
            x( jj+( n-j+1 ) ) = clarnd( 2, iseed )
            IF( abs( x( jj ) ).GT.abs( x( jj+( n-j+1 ) ) ) ) THEN
               x( jj+1 ) = 2.0*x( jj )
            ELSE
               x( jj+1 ) = 2.0*x( jj+( n-j+1 ) )
            END IF
            jj = jj + ( n-j+1 ) + ( n-j )
            j = j + 2
         ELSE IF( j.EQ.n ) THEN
            x( jj ) = clarnd( 2, iseed )
            jj = jj + ( n-j+1 )
            j = j + 1
         END IF
      END IF
*
      RETURN
*
*     End of CLATSP
*
      END