1*> \brief \b ZLAGSY 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 ZLAGSY( 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*> ZLAGSY generates a complex symmetric matrix A, by pre- and post- 29*> multiplying a real diagonal matrix D with a random unitary matrix: 30*> A = U*D*U**T. The semi-bandwidth may then be reduced to k by 31*> additional 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 symmetric 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 November 2011 99* 100*> \ingroup complex16_matgen 101* 102* ===================================================================== 103 SUBROUTINE ZLAGSY( 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 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, II, J, JJ 129 DOUBLE PRECISION WN 130 COMPLEX*16 ALPHA, TAU, WA, WB 131* .. 132* .. External Subroutines .. 133 EXTERNAL XERBLA, ZAXPY, ZGEMV, ZGERC, ZLACGV, ZLARNV, 134 $ ZSCAL, ZSYMV 135* .. 136* .. External Functions .. 137 DOUBLE PRECISION DZNRM2 138 COMPLEX*16 ZDOTC 139 EXTERNAL DZNRM2, ZDOTC 140* .. 141* .. Intrinsic Functions .. 142 INTRINSIC ABS, DBLE, 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( 'ZLAGSY', -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 symmetric matrix 173* 174 DO 60 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 * conjg(u) 194* 195 CALL ZLACGV( N-I+1, WORK, 1 ) 196 CALL ZSYMV( 'Lower', N-I+1, TAU, A( I, I ), LDA, WORK, 1, ZERO, 197 $ WORK( N+1 ), 1 ) 198 CALL ZLACGV( N-I+1, WORK, 1 ) 199* 200* compute v := y - 1/2 * tau * ( u, y ) * u 201* 202 ALPHA = -HALF*TAU*ZDOTC( N-I+1, WORK, 1, WORK( N+1 ), 1 ) 203 CALL ZAXPY( N-I+1, ALPHA, WORK, 1, WORK( N+1 ), 1 ) 204* 205* apply the transformation as a rank-2 update to A(i:n,i:n) 206* 207* CALL ZSYR2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1, 208* $ A( I, I ), LDA ) 209* 210 DO 50 JJ = I, N 211 DO 40 II = JJ, N 212 A( II, JJ ) = A( II, JJ ) - 213 $ WORK( II-I+1 )*WORK( N+JJ-I+1 ) - 214 $ WORK( N+II-I+1 )*WORK( JJ-I+1 ) 215 40 CONTINUE 216 50 CONTINUE 217 60 CONTINUE 218* 219* Reduce number of subdiagonals to K 220* 221 DO 100 I = 1, N - 1 - K 222* 223* generate reflection to annihilate A(k+i+1:n,i) 224* 225 WN = DZNRM2( N-K-I+1, A( K+I, I ), 1 ) 226 WA = ( WN / ABS( A( K+I, I ) ) )*A( K+I, I ) 227 IF( WN.EQ.ZERO ) THEN 228 TAU = ZERO 229 ELSE 230 WB = A( K+I, I ) + WA 231 CALL ZSCAL( N-K-I, ONE / WB, A( K+I+1, I ), 1 ) 232 A( K+I, I ) = ONE 233 TAU = DBLE( WB / WA ) 234 END IF 235* 236* apply reflection to A(k+i:n,i+1:k+i-1) from the left 237* 238 CALL ZGEMV( 'Conjugate transpose', N-K-I+1, K-1, ONE, 239 $ A( K+I, I+1 ), LDA, A( K+I, I ), 1, ZERO, WORK, 1 ) 240 CALL ZGERC( N-K-I+1, K-1, -TAU, A( K+I, I ), 1, WORK, 1, 241 $ A( K+I, I+1 ), LDA ) 242* 243* apply reflection to A(k+i:n,k+i:n) from the left and the right 244* 245* compute y := tau * A * conjg(u) 246* 247 CALL ZLACGV( N-K-I+1, A( K+I, I ), 1 ) 248 CALL ZSYMV( 'Lower', N-K-I+1, TAU, A( K+I, K+I ), LDA, 249 $ A( K+I, I ), 1, ZERO, WORK, 1 ) 250 CALL ZLACGV( N-K-I+1, A( K+I, I ), 1 ) 251* 252* compute v := y - 1/2 * tau * ( u, y ) * u 253* 254 ALPHA = -HALF*TAU*ZDOTC( N-K-I+1, A( K+I, I ), 1, WORK, 1 ) 255 CALL ZAXPY( N-K-I+1, ALPHA, A( K+I, I ), 1, WORK, 1 ) 256* 257* apply symmetric rank-2 update to A(k+i:n,k+i:n) 258* 259* CALL ZSYR2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1, 260* $ A( K+I, K+I ), LDA ) 261* 262 DO 80 JJ = K + I, N 263 DO 70 II = JJ, N 264 A( II, JJ ) = A( II, JJ ) - A( II, I )*WORK( JJ-K-I+1 ) - 265 $ WORK( II-K-I+1 )*A( JJ, I ) 266 70 CONTINUE 267 80 CONTINUE 268* 269 A( K+I, I ) = -WA 270 DO 90 J = K + I + 1, N 271 A( J, I ) = ZERO 272 90 CONTINUE 273 100 CONTINUE 274* 275* Store full symmetric matrix 276* 277 DO 120 J = 1, N 278 DO 110 I = J + 1, N 279 A( J, I ) = A( I, J ) 280 110 CONTINUE 281 120 CONTINUE 282 RETURN 283* 284* End of ZLAGSY 285* 286 END 287