zlatsp.f

Go to the documentation of this file.

*> \brief \b ZLATSP
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZLATSP( UPLO, N, X, ISEED )
* 
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            N
*       ..
*       .. Array Arguments ..
*       INTEGER            ISEED( * )
*       COMPLEX*16         X( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZLATSP 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*16 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 complex16_lin
*
*  =====================================================================
      SUBROUTINE zlatsp( 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*16         x( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX*16         eye
      parameter( eye = ( 0.0d0, 1.0d0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            j, jj, n5
      DOUBLE PRECISION   alpha, alpha3, beta
      COMPLEX*16         a, b, c, r
*     ..
*     .. External Functions ..
      COMPLEX*16         zlarnd
      EXTERNAL           zlarnd
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, sqrt
*     ..
*     .. Executable Statements ..
*
*     Initialize constants
*
      alpha = ( 1.d0+sqrt( 17.d0 ) ) / 8.d0
      beta = alpha - 1.d0 / 1000.d0
      alpha3 = alpha*alpha*alpha
*
*     Fill the matrix with zeros.
*
      DO 10 j = 1, n*( n+1 ) / 2
         x( j ) = 0.0d0
   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*zlarnd( 5, iseed )
            b = zlarnd( 5, iseed ) / alpha
            c = a - 2.d0*b*eye
            r = c / beta
            x( jj ) = a
            x( jj-2 ) = b
            jj = jj - j
            x( jj ) = zlarnd( 2, iseed )
            x( jj-1 ) = r
            jj = jj - ( j-1 )
            x( jj ) = c
            jj = jj - ( j-2 )
            x( jj ) = zlarnd( 2, iseed )
            jj = jj - ( j-3 )
            x( jj ) = zlarnd( 2, iseed )
            IF( abs( x( jj+( j-3 ) ) ).GT.abs( x( jj ) ) ) THEN
               x( jj+( j-4 ) ) = 2.0d0*x( jj+( j-3 ) )
            ELSE
               x( jj+( j-4 ) ) = 2.0d0*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*zlarnd( 5, iseed )
            b = zlarnd( 5, iseed ) / alpha
            c = a - 2.d0*b*eye
            r = c / beta
            x( jj ) = a
            x( jj-2 ) = b
            jj = jj - j
            x( jj ) = zlarnd( 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 ) = zlarnd( 2, iseed )
            x( jj-j ) = zlarnd( 2, iseed )
            IF( abs( x( jj ) ).GT.abs( x( jj-j ) ) ) THEN
               x( jj-1 ) = 2.0d0*x( jj )
            ELSE
               x( jj-1 ) = 2.0d0*x( jj-j )
            END IF
            jj = jj - j - ( j-1 )
            j = j - 2
         ELSE IF( j.EQ.1 ) THEN
            x( jj ) = zlarnd( 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*zlarnd( 5, iseed )
            b = zlarnd( 5, iseed ) / alpha
            c = a - 2.d0*b*eye
            r = c / beta
            x( jj ) = a
            x( jj+2 ) = b
            jj = jj + ( n-j+1 )
            x( jj ) = zlarnd( 2, iseed )
            x( jj+1 ) = r
            jj = jj + ( n-j )
            x( jj ) = c
            jj = jj + ( n-j-1 )
            x( jj ) = zlarnd( 2, iseed )
            jj = jj + ( n-j-2 )
            x( jj ) = zlarnd( 2, iseed )
            IF( abs( x( jj-( n-j-2 ) ) ).GT.abs( x( jj ) ) ) THEN
               x( jj-( n-j-2 )+1 ) = 2.0d0*x( jj-( n-j-2 ) )
            ELSE
               x( jj-( n-j-2 )+1 ) = 2.0d0*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*zlarnd( 5, iseed )
            b = zlarnd( 5, iseed ) / alpha
            c = a - 2.d0*b*eye
            r = c / beta
            x( jj ) = a
            x( jj+2 ) = b
            jj = jj + ( n-j+1 )
            x( jj ) = zlarnd( 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 ) = zlarnd( 2, iseed )
            x( jj+( n-j+1 ) ) = zlarnd( 2, iseed )
            IF( abs( x( jj ) ).GT.abs( x( jj+( n-j+1 ) ) ) ) THEN
               x( jj+1 ) = 2.0d0*x( jj )
            ELSE
               x( jj+1 ) = 2.0d0*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 ) = zlarnd( 2, iseed )
            jj = jj + ( n-j+1 )
            j = j + 1
         END IF
      END IF
*
      return
*
*     End of ZLATSP
*
      END