1*> \brief \b ZLAGHE 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 ZLAGHE( 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* DOUBLE PRECISION D( * ) 19* COMPLEX*16 A( LDA, * ), WORK( * ) 20* .. 21* 22* 23*> \par Purpose: 24* ============= 25*> 26*> \verbatim 27*> 28*> ZLAGHE 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 DOUBLE PRECISION 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*16 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*16 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*> \ingroup complex16_matgen 99* 100* ===================================================================== 101 SUBROUTINE ZLAGHE( N, K, D, A, LDA, ISEED, WORK, INFO ) 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 DOUBLE PRECISION D( * ) 113 COMPLEX*16 A( LDA, * ), WORK( * ) 114* .. 115* 116* ===================================================================== 117* 118* .. Parameters .. 119 COMPLEX*16 ZERO, ONE, HALF 120 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 121 $ ONE = ( 1.0D+0, 0.0D+0 ), 122 $ HALF = ( 0.5D+0, 0.0D+0 ) ) 123* .. 124* .. Local Scalars .. 125 INTEGER I, J 126 DOUBLE PRECISION WN 127 COMPLEX*16 ALPHA, TAU, WA, WB 128* .. 129* .. External Subroutines .. 130 EXTERNAL XERBLA, ZAXPY, ZGEMV, ZGERC, ZHEMV, ZHER2, 131 $ ZLARNV, ZSCAL 132* .. 133* .. External Functions .. 134 DOUBLE PRECISION DZNRM2 135 COMPLEX*16 ZDOTC 136 EXTERNAL DZNRM2, ZDOTC 137* .. 138* .. Intrinsic Functions .. 139 INTRINSIC ABS, DBLE, DCONJG, MAX 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( 'ZLAGHE', -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 ZLARNV( 3, ISEED, N-I+1, WORK ) 176 WN = DZNRM2( 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 ZSCAL( N-I, ONE / WB, WORK( 2 ), 1 ) 183 WORK( 1 ) = ONE 184 TAU = DBLE( 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 ZHEMV( '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*ZDOTC( N-I+1, WORK( N+1 ), 1, WORK, 1 ) 198 CALL ZAXPY( 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 ZHER2( '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 = DZNRM2( 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 ZSCAL( N-K-I, ONE / WB, A( K+I+1, I ), 1 ) 219 A( K+I, I ) = ONE 220 TAU = DBLE( 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 ZGEMV( '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 ZGERC( 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 ZHEMV( '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*ZDOTC( N-K-I+1, WORK, 1, A( K+I, I ), 1 ) 240 CALL ZAXPY( 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 ZHER2( '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 ) = DCONJG( A( I, J ) ) 258 70 CONTINUE 259 80 CONTINUE 260 RETURN 261* 262* End of ZLAGHE 263* 264 END 265