LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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