|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 |
SUBROUTINE ZLATSP( UPLO, N, X, ISEED )
* * -- LAPACK auxiliary test routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER N * .. * .. Array Arguments .. INTEGER ISEED( * ) COMPLEX*16 X( * ) * .. * * Purpose * ======= * * 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. * * Arguments * ========= * * UPLO (input) CHARACTER * Specifies whether the generated matrix is to be upper or * lower triangular. * = 'U': Upper triangular * = 'L': Lower triangular * * N (input) INTEGER * The dimension of the matrix to be generated. * * X (output) 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. * * ISEED (input/output) 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) * * ===================================================================== * * .. 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 |