TESTING/LIN/clatsp.f

   1 *> \brief \b CLATSP
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE CLATSP( UPLO, N, X, ISEED )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       CHARACTER          UPLO
  15 *       INTEGER            N
  16 *       ..
  17 *       .. Array Arguments ..
  18 *       INTEGER            ISEED( * )
  19 *       COMPLEX            X( * )
  20 *       ..
  21 *
  22 *
  23 *> \par Purpose:
  24 *  =============
  25 *>
  26 *> \verbatim
  27 *>
  28 *> CLATSP generates a special test matrix for the complex symmetric
  29 *> (indefinite) factorization for packed matrices.  The pivot blocks of
  30 *> the generated matrix will be in the following order:
  31 *>    2x2 pivot block, non diagonalizable
  32 *>    1x1 pivot block
  33 *>    2x2 pivot block, diagonalizable
  34 *>    (cycle repeats)
  35 *> A row interchange is required for each non-diagonalizable 2x2 block.
  36 *> \endverbatim
  37 *
  38 *  Arguments:
  39 *  ==========
  40 *
  41 *> \param[in] UPLO
  42 *> \verbatim
  43 *>          UPLO is CHARACTER
  44 *>          Specifies whether the generated matrix is to be upper or
  45 *>          lower triangular.
  46 *>          = 'U':  Upper triangular
  47 *>          = 'L':  Lower triangular
  48 *> \endverbatim
  49 *>
  50 *> \param[in] N
  51 *> \verbatim
  52 *>          N is INTEGER
  53 *>          The dimension of the matrix to be generated.
  54 *> \endverbatim
  55 *>
  56 *> \param[out] X
  57 *> \verbatim
  58 *>          X is COMPLEX array, dimension (N*(N+1)/2)
  59 *>          The generated matrix in packed storage format.  The matrix
  60 *>          consists of 3x3 and 2x2 diagonal blocks which result in the
  61 *>          pivot sequence given above.  The matrix outside these
  62 *>          diagonal blocks is zero.
  63 *> \endverbatim
  64 *>
  65 *> \param[in,out] ISEED
  66 *> \verbatim
  67 *>          ISEED is INTEGER array, dimension (4)
  68 *>          On entry, the seed for the random number generator.  The last
  69 *>          of the four integers must be odd.  (modified on exit)
  70 *> \endverbatim
  71 *
  72 *  Authors:
  73 *  ========
  74 *
  75 *> \author Univ. of Tennessee
  76 *> \author Univ. of California Berkeley
  77 *> \author Univ. of Colorado Denver
  78 *> \author NAG Ltd.
  79 *
  80 *> \date November 2011
  81 *
  82 *> \ingroup complex_lin
  83 *
  84 *  =====================================================================
  85       SUBROUTINE CLATSP( UPLO, N, X, ISEED )
  86 *
  87 *  -- LAPACK test routine (version 3.4.0) --
  88 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  89 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  90 *     November 2011
  91 *
  92 *     .. Scalar Arguments ..
  93       CHARACTER          UPLO
  94       INTEGER            N
  95 *     ..
  96 *     .. Array Arguments ..
  97       INTEGER            ISEED( * )
  98       COMPLEX            X( * )
  99 *     ..
 100 *
 101 *  =====================================================================
 102 *
 103 *     .. Parameters ..
 104       COMPLEX            EYE
 105       PARAMETER          ( EYE = ( 0.0, 1.0 ) )
 106 *     ..
 107 *     .. Local Scalars ..
 108       INTEGER            J, JJ, N5
 109       REAL               ALPHA, ALPHA3, BETA
 110       COMPLEX            A, B, C, R
 111 *     ..
 112 *     .. External Functions ..
 113       COMPLEX            CLARND
 114       EXTERNAL           CLARND
 115 *     ..
 116 *     .. Intrinsic Functions ..
 117       INTRINSIC          ABS, SQRT
 118 *     ..
 119 *     .. Executable Statements ..
 120 *
 121 *     Initialize constants
 122 *
 123       ALPHA = ( 1.+SQRT( 17. ) ) / 8.
 124       BETA = ALPHA - 1. / 1000.
 125       ALPHA3 = ALPHA*ALPHA*ALPHA
 126 *
 127 *     Fill the matrix with zeros.
 128 *
 129       DO 10 J = 1, N*( N+1 ) / 2
 130          X( J ) = 0.0
 131    10 CONTINUE
 132 *
 133 *     UPLO = 'U':  Upper triangular storage
 134 *
 135       IF( UPLO.EQ.'U' ) THEN
 136          N5 = N / 5
 137          N5 = N - 5*N5 + 1
 138 *
 139          JJ = N*( N+1 ) / 2
 140          DO 20 J = N, N5, -5
 141             A = ALPHA3*CLARND( 5, ISEED )
 142             B = CLARND( 5, ISEED ) / ALPHA
 143             C = A - 2.*B*EYE
 144             R = C / BETA
 145             X( JJ ) = A
 146             X( JJ-2 ) = B
 147             JJ = JJ - J
 148             X( JJ ) = CLARND( 2, ISEED )
 149             X( JJ-1 ) = R
 150             JJ = JJ - ( J-1 )
 151             X( JJ ) = C
 152             JJ = JJ - ( J-2 )
 153             X( JJ ) = CLARND( 2, ISEED )
 154             JJ = JJ - ( J-3 )
 155             X( JJ ) = CLARND( 2, ISEED )
 156             IF( ABS( X( JJ+( J-3 ) ) ).GT.ABS( X( JJ ) ) ) THEN
 157                X( JJ+( J-4 ) ) = 2.0*X( JJ+( J-3 ) )
 158             ELSE
 159                X( JJ+( J-4 ) ) = 2.0*X( JJ )
 160             END IF
 161             JJ = JJ - ( J-4 )
 162    20    CONTINUE
 163 *
 164 *        Clean-up for N not a multiple of 5.
 165 *
 166          J = N5 - 1
 167          IF( J.GT.2 ) THEN
 168             A = ALPHA3*CLARND( 5, ISEED )
 169             B = CLARND( 5, ISEED ) / ALPHA
 170             C = A - 2.*B*EYE
 171             R = C / BETA
 172             X( JJ ) = A
 173             X( JJ-2 ) = B
 174             JJ = JJ - J
 175             X( JJ ) = CLARND( 2, ISEED )
 176             X( JJ-1 ) = R
 177             JJ = JJ - ( J-1 )
 178             X( JJ ) = C
 179             JJ = JJ - ( J-2 )
 180             J = J - 3
 181          END IF
 182          IF( J.GT.1 ) THEN
 183             X( JJ ) = CLARND( 2, ISEED )
 184             X( JJ-J ) = CLARND( 2, ISEED )
 185             IF( ABS( X( JJ ) ).GT.ABS( X( JJ-J ) ) ) THEN
 186                X( JJ-1 ) = 2.0*X( JJ )
 187             ELSE
 188                X( JJ-1 ) = 2.0*X( JJ-J )
 189             END IF
 190             JJ = JJ - J - ( J-1 )
 191             J = J - 2
 192          ELSE IF( J.EQ.1 ) THEN
 193             X( JJ ) = CLARND( 2, ISEED )
 194             J = J - 1
 195          END IF
 196 *
 197 *     UPLO = 'L':  Lower triangular storage
 198 *
 199       ELSE
 200          N5 = N / 5
 201          N5 = N5*5
 202 *
 203          JJ = 1
 204          DO 30 J = 1, N5, 5
 205             A = ALPHA3*CLARND( 5, ISEED )
 206             B = CLARND( 5, ISEED ) / ALPHA
 207             C = A - 2.*B*EYE
 208             R = C / BETA
 209             X( JJ ) = A
 210             X( JJ+2 ) = B
 211             JJ = JJ + ( N-J+1 )
 212             X( JJ ) = CLARND( 2, ISEED )
 213             X( JJ+1 ) = R
 214             JJ = JJ + ( N-J )
 215             X( JJ ) = C
 216             JJ = JJ + ( N-J-1 )
 217             X( JJ ) = CLARND( 2, ISEED )
 218             JJ = JJ + ( N-J-2 )
 219             X( JJ ) = CLARND( 2, ISEED )
 220             IF( ABS( X( JJ-( N-J-2 ) ) ).GT.ABS( X( JJ ) ) ) THEN
 221                X( JJ-( N-J-2 )+1 ) = 2.0*X( JJ-( N-J-2 ) )
 222             ELSE
 223                X( JJ-( N-J-2 )+1 ) = 2.0*X( JJ )
 224             END IF
 225             JJ = JJ + ( N-J-3 )
 226    30    CONTINUE
 227 *
 228 *        Clean-up for N not a multiple of 5.
 229 *
 230          J = N5 + 1
 231          IF( J.LT.N-1 ) THEN
 232             A = ALPHA3*CLARND( 5, ISEED )
 233             B = CLARND( 5, ISEED ) / ALPHA
 234             C = A - 2.*B*EYE
 235             R = C / BETA
 236             X( JJ ) = A
 237             X( JJ+2 ) = B
 238             JJ = JJ + ( N-J+1 )
 239             X( JJ ) = CLARND( 2, ISEED )
 240             X( JJ+1 ) = R
 241             JJ = JJ + ( N-J )
 242             X( JJ ) = C
 243             JJ = JJ + ( N-J-1 )
 244             J = J + 3
 245          END IF
 246          IF( J.LT.N ) THEN
 247             X( JJ ) = CLARND( 2, ISEED )
 248             X( JJ+( N-J+1 ) ) = CLARND( 2, ISEED )
 249             IF( ABS( X( JJ ) ).GT.ABS( X( JJ+( N-J+1 ) ) ) ) THEN
 250                X( JJ+1 ) = 2.0*X( JJ )
 251             ELSE
 252                X( JJ+1 ) = 2.0*X( JJ+( N-J+1 ) )
 253             END IF
 254             JJ = JJ + ( N-J+1 ) + ( N-J )
 255             J = J + 2
 256          ELSE IF( J.EQ.N ) THEN
 257             X( JJ ) = CLARND( 2, ISEED )
 258             JJ = JJ + ( N-J+1 )
 259             J = J + 1
 260          END IF
 261       END IF
 262 *
 263       RETURN
 264 *
 265 *     End of CLATSP
 266 *
 267       END