MAGMA  1.2.0 MatrixAlgebraonGPUandMulticoreArchitectures
zlatsp.f
Go to the documentation of this file.
1  SUBROUTINE zlatsp( UPLO, N, X, ISEED )
2 *
3 * -- LAPACK auxiliary test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8  CHARACTER uplo
9  INTEGER n
10 * ..
11 * .. Array Arguments ..
12  INTEGER iseed( * )
13  COMPLEX*16 x( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZLATSP generates a special test matrix for the complex symmetric
20 * (indefinite) factorization for packed matrices. The pivot blocks of
21 * the generated matrix will be in the following order:
22 * 2x2 pivot block, non diagonalizable
23 * 1x1 pivot block
24 * 2x2 pivot block, diagonalizable
25 * (cycle repeats)
26 * A row interchange is required for each non-diagonalizable 2x2 block.
27 *
28 * Arguments
29 * =========
30 *
31 * UPLO (input) CHARACTER
32 * Specifies whether the generated matrix is to be upper or
33 * lower triangular.
34 * = 'U': Upper triangular
35 * = 'L': Lower triangular
36 *
37 * N (input) INTEGER
38 * The dimension of the matrix to be generated.
39 *
40 * X (output) COMPLEX*16 array, dimension (N*(N+1)/2)
41 * The generated matrix in packed storage format. The matrix
42 * consists of 3x3 and 2x2 diagonal blocks which result in the
43 * pivot sequence given above. The matrix outside these
44 * diagonal blocks is zero.
45 *
46 * ISEED (input/output) INTEGER array, dimension (4)
47 * On entry, the seed for the random number generator. The last
48 * of the four integers must be odd. (modified on exit)
49 *
50 * =====================================================================
51 *
52 * .. Parameters ..
53  COMPLEX*16 eye
54  parameter( eye = ( 0.0d0, 1.0d0 ) )
55 * ..
56 * .. Local Scalars ..
57  INTEGER j, jj, n5
58  DOUBLE PRECISION alpha, alpha3, beta
59  COMPLEX*16 a, b, c, r
60 * ..
61 * .. External Functions ..
62  COMPLEX*16 zlarnd
63  EXTERNAL zlarnd
64 * ..
65 * .. Intrinsic Functions ..
66  INTRINSIC abs, sqrt
67 * ..
68 * .. Executable Statements ..
69 *
70 * Initialize constants
71 *
72  alpha = ( 1.d0+sqrt( 17.d0 ) ) / 8.d0
73  beta = alpha - 1.d0 / 1000.d0
74  alpha3 = alpha*alpha*alpha
75 *
76 * Fill the matrix with zeros.
77 *
78  DO 10 j = 1, n*( n+1 ) / 2
79  x( j ) = 0.0d0
80  10 continue
81 *
82 * UPLO = 'U': Upper triangular storage
83 *
84  IF( uplo.EQ.'U' ) THEN
85  n5 = n / 5
86  n5 = n - 5*n5 + 1
87 *
88  jj = n*( n+1 ) / 2
89  DO 20 j = n, n5, -5
90  a = alpha3*zlarnd( 5, iseed )
91  b = zlarnd( 5, iseed ) / alpha
92  c = a - 2.d0*b*eye
93  r = c / beta
94  x( jj ) = a
95  x( jj-2 ) = b
96  jj = jj - j
97  x( jj ) = zlarnd( 2, iseed )
98  x( jj-1 ) = r
99  jj = jj - ( j-1 )
100  x( jj ) = c
101  jj = jj - ( j-2 )
102  x( jj ) = zlarnd( 2, iseed )
103  jj = jj - ( j-3 )
104  x( jj ) = zlarnd( 2, iseed )
105  IF( abs( x( jj+( j-3 ) ) ).GT.abs( x( jj ) ) ) THEN
106  x( jj+( j-4 ) ) = 2.0d0*x( jj+( j-3 ) )
107  ELSE
108  x( jj+( j-4 ) ) = 2.0d0*x( jj )
109  END IF
110  jj = jj - ( j-4 )
111  20 continue
112 *
113 * Clean-up for N not a multiple of 5.
114 *
115  j = n5 - 1
116  IF( j.GT.2 ) THEN
117  a = alpha3*zlarnd( 5, iseed )
118  b = zlarnd( 5, iseed ) / alpha
119  c = a - 2.d0*b*eye
120  r = c / beta
121  x( jj ) = a
122  x( jj-2 ) = b
123  jj = jj - j
124  x( jj ) = zlarnd( 2, iseed )
125  x( jj-1 ) = r
126  jj = jj - ( j-1 )
127  x( jj ) = c
128  jj = jj - ( j-2 )
129  j = j - 3
130  END IF
131  IF( j.GT.1 ) THEN
132  x( jj ) = zlarnd( 2, iseed )
133  x( jj-j ) = zlarnd( 2, iseed )
134  IF( abs( x( jj ) ).GT.abs( x( jj-j ) ) ) THEN
135  x( jj-1 ) = 2.0d0*x( jj )
136  ELSE
137  x( jj-1 ) = 2.0d0*x( jj-j )
138  END IF
139  jj = jj - j - ( j-1 )
140  j = j - 2
141  ELSE IF( j.EQ.1 ) THEN
142  x( jj ) = zlarnd( 2, iseed )
143  j = j - 1
144  END IF
145 *
146 * UPLO = 'L': Lower triangular storage
147 *
148  ELSE
149  n5 = n / 5
150  n5 = n5*5
151 *
152  jj = 1
153  DO 30 j = 1, n5, 5
154  a = alpha3*zlarnd( 5, iseed )
155  b = zlarnd( 5, iseed ) / alpha
156  c = a - 2.d0*b*eye
157  r = c / beta
158  x( jj ) = a
159  x( jj+2 ) = b
160  jj = jj + ( n-j+1 )
161  x( jj ) = zlarnd( 2, iseed )
162  x( jj+1 ) = r
163  jj = jj + ( n-j )
164  x( jj ) = c
165  jj = jj + ( n-j-1 )
166  x( jj ) = zlarnd( 2, iseed )
167  jj = jj + ( n-j-2 )
168  x( jj ) = zlarnd( 2, iseed )
169  IF( abs( x( jj-( n-j-2 ) ) ).GT.abs( x( jj ) ) ) THEN
170  x( jj-( n-j-2 )+1 ) = 2.0d0*x( jj-( n-j-2 ) )
171  ELSE
172  x( jj-( n-j-2 )+1 ) = 2.0d0*x( jj )
173  END IF
174  jj = jj + ( n-j-3 )
175  30 continue
176 *
177 * Clean-up for N not a multiple of 5.
178 *
179  j = n5 + 1
180  IF( j.LT.n-1 ) THEN
181  a = alpha3*zlarnd( 5, iseed )
182  b = zlarnd( 5, iseed ) / alpha
183  c = a - 2.d0*b*eye
184  r = c / beta
185  x( jj ) = a
186  x( jj+2 ) = b
187  jj = jj + ( n-j+1 )
188  x( jj ) = zlarnd( 2, iseed )
189  x( jj+1 ) = r
190  jj = jj + ( n-j )
191  x( jj ) = c
192  jj = jj + ( n-j-1 )
193  j = j + 3
194  END IF
195  IF( j.LT.n ) THEN
196  x( jj ) = zlarnd( 2, iseed )
197  x( jj+( n-j+1 ) ) = zlarnd( 2, iseed )
198  IF( abs( x( jj ) ).GT.abs( x( jj+( n-j+1 ) ) ) ) THEN
199  x( jj+1 ) = 2.0d0*x( jj )
200  ELSE
201  x( jj+1 ) = 2.0d0*x( jj+( n-j+1 ) )
202  END IF
203  jj = jj + ( n-j+1 ) + ( n-j )
204  j = j + 2
205  ELSE IF( j.EQ.n ) THEN
206  x( jj ) = zlarnd( 2, iseed )
207  jj = jj + ( n-j+1 )
208  j = j + 1
209  END IF
210  END IF
211 *
212  return
213 *
214 * End of ZLATSP
215 *
216  END