LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
clatsp.f
Go to the documentation of this file.
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*> \ingroup complex_lin
81*
82* =====================================================================
83 SUBROUTINE clatsp( UPLO, N, X, ISEED )
84*
85* -- LAPACK test routine --
86* -- LAPACK is a software package provided by Univ. of Tennessee, --
87* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
88*
89* .. Scalar Arguments ..
90 CHARACTER UPLO
91 INTEGER N
92* ..
93* .. Array Arguments ..
94 INTEGER ISEED( * )
95 COMPLEX X( * )
96* ..
97*
98* =====================================================================
99*
100* .. Parameters ..
101 COMPLEX EYE
102 parameter( eye = ( 0.0, 1.0 ) )
103* ..
104* .. Local Scalars ..
105 INTEGER J, JJ, N5
106 REAL ALPHA, ALPHA3, BETA
107 COMPLEX A, B, C, R
108* ..
109* .. External Functions ..
110 COMPLEX CLARND
111 EXTERNAL clarnd
112* ..
113* .. Intrinsic Functions ..
114 INTRINSIC abs, sqrt
115* ..
116* .. Executable Statements ..
117*
118* Initialize constants
119*
120 alpha = ( 1.+sqrt( 17. ) ) / 8.
121 beta = alpha - 1. / 1000.
122 alpha3 = alpha*alpha*alpha
123*
124* Fill the matrix with zeros.
125*
126 DO 10 j = 1, n*( n+1 ) / 2
127 x( j ) = 0.0
128 10 CONTINUE
129*
130* UPLO = 'U': Upper triangular storage
131*
132 IF( uplo.EQ.'U' ) THEN
133 n5 = n / 5
134 n5 = n - 5*n5 + 1
135*
136 jj = n*( n+1 ) / 2
137 DO 20 j = n, n5, -5
138 a = alpha3*clarnd( 5, iseed )
139 b = clarnd( 5, iseed ) / alpha
140 c = a - 2.*b*eye
141 r = c / beta
142 x( jj ) = a
143 x( jj-2 ) = b
144 jj = jj - j
145 x( jj ) = clarnd( 2, iseed )
146 x( jj-1 ) = r
147 jj = jj - ( j-1 )
148 x( jj ) = c
149 jj = jj - ( j-2 )
150 x( jj ) = clarnd( 2, iseed )
151 jj = jj - ( j-3 )
152 x( jj ) = clarnd( 2, iseed )
153 IF( abs( x( jj+( j-3 ) ) ).GT.abs( x( jj ) ) ) THEN
154 x( jj+( j-4 ) ) = 2.0*x( jj+( j-3 ) )
155 ELSE
156 x( jj+( j-4 ) ) = 2.0*x( jj )
157 END IF
158 jj = jj - ( j-4 )
159 20 CONTINUE
160*
161* Clean-up for N not a multiple of 5.
162*
163 j = n5 - 1
164 IF( j.GT.2 ) THEN
165 a = alpha3*clarnd( 5, iseed )
166 b = clarnd( 5, iseed ) / alpha
167 c = a - 2.*b*eye
168 r = c / beta
169 x( jj ) = a
170 x( jj-2 ) = b
171 jj = jj - j
172 x( jj ) = clarnd( 2, iseed )
173 x( jj-1 ) = r
174 jj = jj - ( j-1 )
175 x( jj ) = c
176 jj = jj - ( j-2 )
177 j = j - 3
178 END IF
179 IF( j.GT.1 ) THEN
180 x( jj ) = clarnd( 2, iseed )
181 x( jj-j ) = clarnd( 2, iseed )
182 IF( abs( x( jj ) ).GT.abs( x( jj-j ) ) ) THEN
183 x( jj-1 ) = 2.0*x( jj )
184 ELSE
185 x( jj-1 ) = 2.0*x( jj-j )
186 END IF
187 jj = jj - j - ( j-1 )
188 j = j - 2
189 ELSE IF( j.EQ.1 ) THEN
190 x( jj ) = clarnd( 2, iseed )
191 j = j - 1
192 END IF
193*
194* UPLO = 'L': Lower triangular storage
195*
196 ELSE
197 n5 = n / 5
198 n5 = n5*5
199*
200 jj = 1
201 DO 30 j = 1, n5, 5
202 a = alpha3*clarnd( 5, iseed )
203 b = clarnd( 5, iseed ) / alpha
204 c = a - 2.*b*eye
205 r = c / beta
206 x( jj ) = a
207 x( jj+2 ) = b
208 jj = jj + ( n-j+1 )
209 x( jj ) = clarnd( 2, iseed )
210 x( jj+1 ) = r
211 jj = jj + ( n-j )
212 x( jj ) = c
213 jj = jj + ( n-j-1 )
214 x( jj ) = clarnd( 2, iseed )
215 jj = jj + ( n-j-2 )
216 x( jj ) = clarnd( 2, iseed )
217 IF( abs( x( jj-( n-j-2 ) ) ).GT.abs( x( jj ) ) ) THEN
218 x( jj-( n-j-2 )+1 ) = 2.0*x( jj-( n-j-2 ) )
219 ELSE
220 x( jj-( n-j-2 )+1 ) = 2.0*x( jj )
221 END IF
222 jj = jj + ( n-j-3 )
223 30 CONTINUE
224*
225* Clean-up for N not a multiple of 5.
226*
227 j = n5 + 1
228 IF( j.LT.n-1 ) THEN
229 a = alpha3*clarnd( 5, iseed )
230 b = clarnd( 5, iseed ) / alpha
231 c = a - 2.*b*eye
232 r = c / beta
233 x( jj ) = a
234 x( jj+2 ) = b
235 jj = jj + ( n-j+1 )
236 x( jj ) = clarnd( 2, iseed )
237 x( jj+1 ) = r
238 jj = jj + ( n-j )
239 x( jj ) = c
240 jj = jj + ( n-j-1 )
241 j = j + 3
242 END IF
243 IF( j.LT.n ) THEN
244 x( jj ) = clarnd( 2, iseed )
245 x( jj+( n-j+1 ) ) = clarnd( 2, iseed )
246 IF( abs( x( jj ) ).GT.abs( x( jj+( n-j+1 ) ) ) ) THEN
247 x( jj+1 ) = 2.0*x( jj )
248 ELSE
249 x( jj+1 ) = 2.0*x( jj+( n-j+1 ) )
250 END IF
251 jj = jj + ( n-j+1 ) + ( n-j )
252 j = j + 2
253 ELSE IF( j.EQ.n ) THEN
254 x( jj ) = clarnd( 2, iseed )
255 jj = jj + ( n-j+1 )
256 j = j + 1
257 END IF
258 END IF
259*
260 RETURN
261*
262* End of CLATSP
263*
264 END
subroutine clatsp(UPLO, N, X, ISEED)
CLATSP
Definition: clatsp.f:84