LAPACK  3.8.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 *> \date December 2016
81 *
82 *> \ingroup complex_lin
83 *
84 * =====================================================================
85  SUBROUTINE clatsp( UPLO, N, X, ISEED )
86 *
87 * -- LAPACK test routine (version 3.7.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 * December 2016
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
subroutine clatsp(UPLO, N, X, ISEED)
CLATSP
Definition: clatsp.f:86