LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
claghe.f
Go to the documentation of this file.
1 *> \brief \b CLAGHE
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 CLAGHE( N, K, D, A, LDA, ISEED, WORK, INFO )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER INFO, K, LDA, N
15 * ..
16 * .. Array Arguments ..
17 * INTEGER ISEED( 4 )
18 * REAL D( * )
19 * COMPLEX A( LDA, * ), WORK( * )
20 * ..
21 *
22 *
23 *> \par Purpose:
24 * =============
25 *>
26 *> \verbatim
27 *>
28 *> CLAGHE generates a complex hermitian matrix A, by pre- and post-
29 *> multiplying a real diagonal matrix D with a random unitary matrix:
30 *> A = U*D*U'. The semi-bandwidth may then be reduced to k by additional
31 *> unitary transformations.
32 *> \endverbatim
33 *
34 * Arguments:
35 * ==========
36 *
37 *> \param[in] N
38 *> \verbatim
39 *> N is INTEGER
40 *> The order of the matrix A. N >= 0.
41 *> \endverbatim
42 *>
43 *> \param[in] K
44 *> \verbatim
45 *> K is INTEGER
46 *> The number of nonzero subdiagonals within the band of A.
47 *> 0 <= K <= N-1.
48 *> \endverbatim
49 *>
50 *> \param[in] D
51 *> \verbatim
52 *> D is REAL array, dimension (N)
53 *> The diagonal elements of the diagonal matrix D.
54 *> \endverbatim
55 *>
56 *> \param[out] A
57 *> \verbatim
58 *> A is COMPLEX array, dimension (LDA,N)
59 *> The generated n by n hermitian matrix A (the full matrix is
60 *> stored).
61 *> \endverbatim
62 *>
63 *> \param[in] LDA
64 *> \verbatim
65 *> LDA is INTEGER
66 *> The leading dimension of the array A. LDA >= N.
67 *> \endverbatim
68 *>
69 *> \param[in,out] ISEED
70 *> \verbatim
71 *> ISEED is INTEGER array, dimension (4)
72 *> On entry, the seed of the random number generator; the array
73 *> elements must be between 0 and 4095, and ISEED(4) must be
74 *> odd.
75 *> On exit, the seed is updated.
76 *> \endverbatim
77 *>
78 *> \param[out] WORK
79 *> \verbatim
80 *> WORK is COMPLEX array, dimension (2*N)
81 *> \endverbatim
82 *>
83 *> \param[out] INFO
84 *> \verbatim
85 *> INFO is INTEGER
86 *> = 0: successful exit
87 *> < 0: if INFO = -i, the i-th argument had an illegal value
88 *> \endverbatim
89 *
90 * Authors:
91 * ========
92 *
93 *> \author Univ. of Tennessee
94 *> \author Univ. of California Berkeley
95 *> \author Univ. of Colorado Denver
96 *> \author NAG Ltd.
97 *
98 *> \date November 2011
99 *
100 *> \ingroup complex_matgen
101 *
102 * =====================================================================
103  SUBROUTINE claghe( N, K, D, A, LDA, ISEED, WORK, INFO )
104 *
105 * -- LAPACK auxiliary routine (version 3.4.0) --
106 * -- LAPACK is a software package provided by Univ. of Tennessee, --
107 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
108 * November 2011
109 *
110 * .. Scalar Arguments ..
111  INTEGER info, k, lda, n
112 * ..
113 * .. Array Arguments ..
114  INTEGER iseed( 4 )
115  REAL d( * )
116  COMPLEX a( lda, * ), work( * )
117 * ..
118 *
119 * =====================================================================
120 *
121 * .. Parameters ..
122  COMPLEX zero, one, half
123  parameter( zero = ( 0.0e+0, 0.0e+0 ),
124  $ one = ( 1.0e+0, 0.0e+0 ),
125  $ half = ( 0.5e+0, 0.0e+0 ) )
126 * ..
127 * .. Local Scalars ..
128  INTEGER i, j
129  REAL wn
130  COMPLEX alpha, tau, wa, wb
131 * ..
132 * .. External Subroutines ..
133  EXTERNAL caxpy, cgemv, cgerc, chemv, cher2, clarnv,
134  $ cscal, xerbla
135 * ..
136 * .. External Functions ..
137  REAL scnrm2
138  COMPLEX cdotc
139  EXTERNAL scnrm2, cdotc
140 * ..
141 * .. Intrinsic Functions ..
142  INTRINSIC abs, conjg, max, real
143 * ..
144 * .. Executable Statements ..
145 *
146 * Test the input arguments
147 *
148  info = 0
149  IF( n.LT.0 ) THEN
150  info = -1
151  ELSE IF( k.LT.0 .OR. k.GT.n-1 ) THEN
152  info = -2
153  ELSE IF( lda.LT.max( 1, n ) ) THEN
154  info = -5
155  END IF
156  IF( info.LT.0 ) THEN
157  CALL xerbla( 'CLAGHE', -info )
158  return
159  END IF
160 *
161 * initialize lower triangle of A to diagonal matrix
162 *
163  DO 20 j = 1, n
164  DO 10 i = j + 1, n
165  a( i, j ) = zero
166  10 continue
167  20 continue
168  DO 30 i = 1, n
169  a( i, i ) = d( i )
170  30 continue
171 *
172 * Generate lower triangle of hermitian matrix
173 *
174  DO 40 i = n - 1, 1, -1
175 *
176 * generate random reflection
177 *
178  CALL clarnv( 3, iseed, n-i+1, work )
179  wn = scnrm2( n-i+1, work, 1 )
180  wa = ( wn / abs( work( 1 ) ) )*work( 1 )
181  IF( wn.EQ.zero ) THEN
182  tau = zero
183  ELSE
184  wb = work( 1 ) + wa
185  CALL cscal( n-i, one / wb, work( 2 ), 1 )
186  work( 1 ) = one
187  tau = REAL( wb / wa )
188  END IF
189 *
190 * apply random reflection to A(i:n,i:n) from the left
191 * and the right
192 *
193 * compute y := tau * A * u
194 *
195  CALL chemv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
196  $ work( n+1 ), 1 )
197 *
198 * compute v := y - 1/2 * tau * ( y, u ) * u
199 *
200  alpha = -half*tau*cdotc( n-i+1, work( n+1 ), 1, work, 1 )
201  CALL caxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
202 *
203 * apply the transformation as a rank-2 update to A(i:n,i:n)
204 *
205  CALL cher2( 'Lower', n-i+1, -one, work, 1, work( n+1 ), 1,
206  $ a( i, i ), lda )
207  40 continue
208 *
209 * Reduce number of subdiagonals to K
210 *
211  DO 60 i = 1, n - 1 - k
212 *
213 * generate reflection to annihilate A(k+i+1:n,i)
214 *
215  wn = scnrm2( n-k-i+1, a( k+i, i ), 1 )
216  wa = ( wn / abs( a( k+i, i ) ) )*a( k+i, i )
217  IF( wn.EQ.zero ) THEN
218  tau = zero
219  ELSE
220  wb = a( k+i, i ) + wa
221  CALL cscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
222  a( k+i, i ) = one
223  tau = REAL( wb / wa )
224  END IF
225 *
226 * apply reflection to A(k+i:n,i+1:k+i-1) from the left
227 *
228  CALL cgemv( 'Conjugate transpose', n-k-i+1, k-1, one,
229  $ a( k+i, i+1 ), lda, a( k+i, i ), 1, zero, work, 1 )
230  CALL cgerc( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
231  $ a( k+i, i+1 ), lda )
232 *
233 * apply reflection to A(k+i:n,k+i:n) from the left and the right
234 *
235 * compute y := tau * A * u
236 *
237  CALL chemv( 'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
238  $ a( k+i, i ), 1, zero, work, 1 )
239 *
240 * compute v := y - 1/2 * tau * ( y, u ) * u
241 *
242  alpha = -half*tau*cdotc( n-k-i+1, work, 1, a( k+i, i ), 1 )
243  CALL caxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
244 *
245 * apply hermitian rank-2 update to A(k+i:n,k+i:n)
246 *
247  CALL cher2( 'Lower', n-k-i+1, -one, a( k+i, i ), 1, work, 1,
248  $ a( k+i, k+i ), lda )
249 *
250  a( k+i, i ) = -wa
251  DO 50 j = k + i + 1, n
252  a( j, i ) = zero
253  50 continue
254  60 continue
255 *
256 * Store full hermitian matrix
257 *
258  DO 80 j = 1, n
259  DO 70 i = j + 1, n
260  a( j, i ) = conjg( a( i, j ) )
261  70 continue
262  80 continue
263  return
264 *
265 * End of CLAGHE
266 *
267  END