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*> \date December 2016 99* 100*> \ingroup complex16_matgen 101* 102* ===================================================================== 103 SUBROUTINE ZLAGHE( N, K, D, A, LDA, ISEED, WORK, INFO ) 104* 105* -- LAPACK auxiliary routine (version 3.7.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* December 2016 109* 110* .. Scalar Arguments .. 111 INTEGER INFO, K, LDA, N 112* .. 113* .. Array Arguments .. 114 INTEGER ISEED( 4 ) 115 DOUBLE PRECISION D( * ) 116 COMPLEX*16 A( LDA, * ), WORK( * ) 117* .. 118* 119* ===================================================================== 120* 121* .. Parameters .. 122 COMPLEX*16 ZERO, ONE, HALF 123 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 124 $ ONE = ( 1.0D+0, 0.0D+0 ), 125 $ HALF = ( 0.5D+0, 0.0D+0 ) ) 126* .. 127* .. Local Scalars .. 128 INTEGER I, J 129 DOUBLE PRECISION WN 130 COMPLEX*16 ALPHA, TAU, WA, WB 131* .. 132* .. External Subroutines .. 133 EXTERNAL XERBLA, ZAXPY, ZGEMV, ZGERC, ZHEMV, ZHER2, 134 $ ZLARNV, ZSCAL 135* .. 136* .. External Functions .. 137 DOUBLE PRECISION DZNRM2 138 COMPLEX*16 ZDOTC 139 EXTERNAL DZNRM2, ZDOTC 140* .. 141* .. Intrinsic Functions .. 142 INTRINSIC ABS, DBLE, DCONJG, MAX 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( 'ZLAGHE', -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 ZLARNV( 3, ISEED, N-I+1, WORK ) 179 WN = DZNRM2( 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 ZSCAL( N-I, ONE / WB, WORK( 2 ), 1 ) 186 WORK( 1 ) = ONE 187 TAU = DBLE( 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 ZHEMV( '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*ZDOTC( N-I+1, WORK( N+1 ), 1, WORK, 1 ) 201 CALL ZAXPY( 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 ZHER2( '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 = DZNRM2( 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 ZSCAL( N-K-I, ONE / WB, A( K+I+1, I ), 1 ) 222 A( K+I, I ) = ONE 223 TAU = DBLE( 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 ZGEMV( '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 ZGERC( 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 ZHEMV( '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*ZDOTC( N-K-I+1, WORK, 1, A( K+I, I ), 1 ) 243 CALL ZAXPY( 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 ZHER2( '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 ) = DCONJG( A( I, J ) ) 261 70 CONTINUE 262 80 CONTINUE 263 RETURN 264* 265* End of ZLAGHE 266* 267 END 268