LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ claghe()

subroutine claghe ( integer  N,
integer  K,
real, dimension( * )  D,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( 4 )  ISEED,
complex, dimension( * )  WORK,
integer  INFO 
)

CLAGHE

Purpose:
 CLAGHE generates a complex hermitian matrix A, by pre- and post-
 multiplying a real diagonal matrix D with a random unitary matrix:
 A = U*D*U'. The semi-bandwidth may then be reduced to k by additional
 unitary transformations.
Parameters
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]K
          K is INTEGER
          The number of nonzero subdiagonals within the band of A.
          0 <= K <= N-1.
[in]D
          D is REAL array, dimension (N)
          The diagonal elements of the diagonal matrix D.
[out]A
          A is COMPLEX array, dimension (LDA,N)
          The generated n by n hermitian matrix A (the full matrix is
          stored).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= N.
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry, the seed of the random number generator; the array
          elements must be between 0 and 4095, and ISEED(4) must be
          odd.
          On exit, the seed is updated.
[out]WORK
          WORK is COMPLEX array, dimension (2*N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 101 of file claghe.f.

102 *
103 * -- LAPACK auxiliary routine --
104 * -- LAPACK is a software package provided by Univ. of Tennessee, --
105 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
106 *
107 * .. Scalar Arguments ..
108  INTEGER INFO, K, LDA, N
109 * ..
110 * .. Array Arguments ..
111  INTEGER ISEED( 4 )
112  REAL D( * )
113  COMPLEX A( LDA, * ), WORK( * )
114 * ..
115 *
116 * =====================================================================
117 *
118 * .. Parameters ..
119  COMPLEX ZERO, ONE, HALF
120  parameter( zero = ( 0.0e+0, 0.0e+0 ),
121  $ one = ( 1.0e+0, 0.0e+0 ),
122  $ half = ( 0.5e+0, 0.0e+0 ) )
123 * ..
124 * .. Local Scalars ..
125  INTEGER I, J
126  REAL WN
127  COMPLEX ALPHA, TAU, WA, WB
128 * ..
129 * .. External Subroutines ..
130  EXTERNAL caxpy, cgemv, cgerc, chemv, cher2, clarnv,
131  $ cscal, xerbla
132 * ..
133 * .. External Functions ..
134  REAL SCNRM2
135  COMPLEX CDOTC
136  EXTERNAL scnrm2, cdotc
137 * ..
138 * .. Intrinsic Functions ..
139  INTRINSIC abs, conjg, max, real
140 * ..
141 * .. Executable Statements ..
142 *
143 * Test the input arguments
144 *
145  info = 0
146  IF( n.LT.0 ) THEN
147  info = -1
148  ELSE IF( k.LT.0 .OR. k.GT.n-1 ) THEN
149  info = -2
150  ELSE IF( lda.LT.max( 1, n ) ) THEN
151  info = -5
152  END IF
153  IF( info.LT.0 ) THEN
154  CALL xerbla( 'CLAGHE', -info )
155  RETURN
156  END IF
157 *
158 * initialize lower triangle of A to diagonal matrix
159 *
160  DO 20 j = 1, n
161  DO 10 i = j + 1, n
162  a( i, j ) = zero
163  10 CONTINUE
164  20 CONTINUE
165  DO 30 i = 1, n
166  a( i, i ) = d( i )
167  30 CONTINUE
168 *
169 * Generate lower triangle of hermitian matrix
170 *
171  DO 40 i = n - 1, 1, -1
172 *
173 * generate random reflection
174 *
175  CALL clarnv( 3, iseed, n-i+1, work )
176  wn = scnrm2( n-i+1, work, 1 )
177  wa = ( wn / abs( work( 1 ) ) )*work( 1 )
178  IF( wn.EQ.zero ) THEN
179  tau = zero
180  ELSE
181  wb = work( 1 ) + wa
182  CALL cscal( n-i, one / wb, work( 2 ), 1 )
183  work( 1 ) = one
184  tau = real( wb / wa )
185  END IF
186 *
187 * apply random reflection to A(i:n,i:n) from the left
188 * and the right
189 *
190 * compute y := tau * A * u
191 *
192  CALL chemv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
193  $ work( n+1 ), 1 )
194 *
195 * compute v := y - 1/2 * tau * ( y, u ) * u
196 *
197  alpha = -half*tau*cdotc( n-i+1, work( n+1 ), 1, work, 1 )
198  CALL caxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
199 *
200 * apply the transformation as a rank-2 update to A(i:n,i:n)
201 *
202  CALL cher2( 'Lower', n-i+1, -one, work, 1, work( n+1 ), 1,
203  $ a( i, i ), lda )
204  40 CONTINUE
205 *
206 * Reduce number of subdiagonals to K
207 *
208  DO 60 i = 1, n - 1 - k
209 *
210 * generate reflection to annihilate A(k+i+1:n,i)
211 *
212  wn = scnrm2( n-k-i+1, a( k+i, i ), 1 )
213  wa = ( wn / abs( a( k+i, i ) ) )*a( k+i, i )
214  IF( wn.EQ.zero ) THEN
215  tau = zero
216  ELSE
217  wb = a( k+i, i ) + wa
218  CALL cscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
219  a( k+i, i ) = one
220  tau = real( wb / wa )
221  END IF
222 *
223 * apply reflection to A(k+i:n,i+1:k+i-1) from the left
224 *
225  CALL cgemv( 'Conjugate transpose', n-k-i+1, k-1, one,
226  $ a( k+i, i+1 ), lda, a( k+i, i ), 1, zero, work, 1 )
227  CALL cgerc( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
228  $ a( k+i, i+1 ), lda )
229 *
230 * apply reflection to A(k+i:n,k+i:n) from the left and the right
231 *
232 * compute y := tau * A * u
233 *
234  CALL chemv( 'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
235  $ a( k+i, i ), 1, zero, work, 1 )
236 *
237 * compute v := y - 1/2 * tau * ( y, u ) * u
238 *
239  alpha = -half*tau*cdotc( n-k-i+1, work, 1, a( k+i, i ), 1 )
240  CALL caxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
241 *
242 * apply hermitian rank-2 update to A(k+i:n,k+i:n)
243 *
244  CALL cher2( 'Lower', n-k-i+1, -one, a( k+i, i ), 1, work, 1,
245  $ a( k+i, k+i ), lda )
246 *
247  a( k+i, i ) = -wa
248  DO 50 j = k + i + 1, n
249  a( j, i ) = zero
250  50 CONTINUE
251  60 CONTINUE
252 *
253 * Store full hermitian matrix
254 *
255  DO 80 j = 1, n
256  DO 70 i = j + 1, n
257  a( j, i ) = conjg( a( i, j ) )
258  70 CONTINUE
259  80 CONTINUE
260  RETURN
261 *
262 * End of CLAGHE
263 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:83
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERC
Definition: cgerc.f:130
subroutine cher2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CHER2
Definition: cher2.f:150
subroutine chemv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CHEMV
Definition: chemv.f:154
subroutine clarnv(IDIST, ISEED, N, X)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: clarnv.f:99
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition: scnrm2.f90:90
Here is the call graph for this function:
Here is the caller graph for this function: