1*> \brief \b ZLAGGE 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 ZLAGGE( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO ) 12* 13* .. Scalar Arguments .. 14* INTEGER INFO, KL, KU, LDA, M, 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*> ZLAGGE generates a complex general m by n matrix A, by pre- and post- 29*> multiplying a real diagonal matrix D with random unitary matrices: 30*> A = U*D*V. The lower and upper bandwidths may then be reduced to 31*> kl and ku by additional unitary transformations. 32*> \endverbatim 33* 34* Arguments: 35* ========== 36* 37*> \param[in] M 38*> \verbatim 39*> M is INTEGER 40*> The number of rows of the matrix A. M >= 0. 41*> \endverbatim 42*> 43*> \param[in] N 44*> \verbatim 45*> N is INTEGER 46*> The number of columns of the matrix A. N >= 0. 47*> \endverbatim 48*> 49*> \param[in] KL 50*> \verbatim 51*> KL is INTEGER 52*> The number of nonzero subdiagonals within the band of A. 53*> 0 <= KL <= M-1. 54*> \endverbatim 55*> 56*> \param[in] KU 57*> \verbatim 58*> KU is INTEGER 59*> The number of nonzero superdiagonals within the band of A. 60*> 0 <= KU <= N-1. 61*> \endverbatim 62*> 63*> \param[in] D 64*> \verbatim 65*> D is DOUBLE PRECISION array, dimension (min(M,N)) 66*> The diagonal elements of the diagonal matrix D. 67*> \endverbatim 68*> 69*> \param[out] A 70*> \verbatim 71*> A is COMPLEX*16 array, dimension (LDA,N) 72*> The generated m by n matrix A. 73*> \endverbatim 74*> 75*> \param[in] LDA 76*> \verbatim 77*> LDA is INTEGER 78*> The leading dimension of the array A. LDA >= M. 79*> \endverbatim 80*> 81*> \param[in,out] ISEED 82*> \verbatim 83*> ISEED is INTEGER array, dimension (4) 84*> On entry, the seed of the random number generator; the array 85*> elements must be between 0 and 4095, and ISEED(4) must be 86*> odd. 87*> On exit, the seed is updated. 88*> \endverbatim 89*> 90*> \param[out] WORK 91*> \verbatim 92*> WORK is COMPLEX*16 array, dimension (M+N) 93*> \endverbatim 94*> 95*> \param[out] INFO 96*> \verbatim 97*> INFO is INTEGER 98*> = 0: successful exit 99*> < 0: if INFO = -i, the i-th argument had an illegal value 100*> \endverbatim 101* 102* Authors: 103* ======== 104* 105*> \author Univ. of Tennessee 106*> \author Univ. of California Berkeley 107*> \author Univ. of Colorado Denver 108*> \author NAG Ltd. 109* 110*> \ingroup complex16_matgen 111* 112* ===================================================================== 113 SUBROUTINE ZLAGGE( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO ) 114* 115* -- LAPACK auxiliary routine -- 116* -- LAPACK is a software package provided by Univ. of Tennessee, -- 117* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 118* 119* .. Scalar Arguments .. 120 INTEGER INFO, KL, KU, LDA, M, N 121* .. 122* .. Array Arguments .. 123 INTEGER ISEED( 4 ) 124 DOUBLE PRECISION D( * ) 125 COMPLEX*16 A( LDA, * ), WORK( * ) 126* .. 127* 128* ===================================================================== 129* 130* .. Parameters .. 131 COMPLEX*16 ZERO, ONE 132 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 133 $ ONE = ( 1.0D+0, 0.0D+0 ) ) 134* .. 135* .. Local Scalars .. 136 INTEGER I, J 137 DOUBLE PRECISION WN 138 COMPLEX*16 TAU, WA, WB 139* .. 140* .. External Subroutines .. 141 EXTERNAL XERBLA, ZGEMV, ZGERC, ZLACGV, ZLARNV, ZSCAL 142* .. 143* .. Intrinsic Functions .. 144 INTRINSIC ABS, DBLE, MAX, MIN 145* .. 146* .. External Functions .. 147 DOUBLE PRECISION DZNRM2 148 EXTERNAL DZNRM2 149* .. 150* .. Executable Statements .. 151* 152* Test the input arguments 153* 154 INFO = 0 155 IF( M.LT.0 ) THEN 156 INFO = -1 157 ELSE IF( N.LT.0 ) THEN 158 INFO = -2 159 ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN 160 INFO = -3 161 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 162 INFO = -4 163 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 164 INFO = -7 165 END IF 166 IF( INFO.LT.0 ) THEN 167 CALL XERBLA( 'ZLAGGE', -INFO ) 168 RETURN 169 END IF 170* 171* initialize A to diagonal matrix 172* 173 DO 20 J = 1, N 174 DO 10 I = 1, M 175 A( I, J ) = ZERO 176 10 CONTINUE 177 20 CONTINUE 178 DO 30 I = 1, MIN( M, N ) 179 A( I, I ) = D( I ) 180 30 CONTINUE 181* 182* Quick exit if the user wants a diagonal matrix 183* 184 IF(( KL .EQ. 0 ).AND.( KU .EQ. 0)) RETURN 185* 186* pre- and post-multiply A by random unitary matrices 187* 188 DO 40 I = MIN( M, N ), 1, -1 189 IF( I.LT.M ) THEN 190* 191* generate random reflection 192* 193 CALL ZLARNV( 3, ISEED, M-I+1, WORK ) 194 WN = DZNRM2( M-I+1, WORK, 1 ) 195 WA = ( WN / ABS( WORK( 1 ) ) )*WORK( 1 ) 196 IF( WN.EQ.ZERO ) THEN 197 TAU = ZERO 198 ELSE 199 WB = WORK( 1 ) + WA 200 CALL ZSCAL( M-I, ONE / WB, WORK( 2 ), 1 ) 201 WORK( 1 ) = ONE 202 TAU = DBLE( WB / WA ) 203 END IF 204* 205* multiply A(i:m,i:n) by random reflection from the left 206* 207 CALL ZGEMV( 'Conjugate transpose', M-I+1, N-I+1, ONE, 208 $ A( I, I ), LDA, WORK, 1, ZERO, WORK( M+1 ), 1 ) 209 CALL ZGERC( M-I+1, N-I+1, -TAU, WORK, 1, WORK( M+1 ), 1, 210 $ A( I, I ), LDA ) 211 END IF 212 IF( I.LT.N ) THEN 213* 214* generate random reflection 215* 216 CALL ZLARNV( 3, ISEED, N-I+1, WORK ) 217 WN = DZNRM2( N-I+1, WORK, 1 ) 218 WA = ( WN / ABS( WORK( 1 ) ) )*WORK( 1 ) 219 IF( WN.EQ.ZERO ) THEN 220 TAU = ZERO 221 ELSE 222 WB = WORK( 1 ) + WA 223 CALL ZSCAL( N-I, ONE / WB, WORK( 2 ), 1 ) 224 WORK( 1 ) = ONE 225 TAU = DBLE( WB / WA ) 226 END IF 227* 228* multiply A(i:m,i:n) by random reflection from the right 229* 230 CALL ZGEMV( 'No transpose', M-I+1, N-I+1, ONE, A( I, I ), 231 $ LDA, WORK, 1, ZERO, WORK( N+1 ), 1 ) 232 CALL ZGERC( M-I+1, N-I+1, -TAU, WORK( N+1 ), 1, WORK, 1, 233 $ A( I, I ), LDA ) 234 END IF 235 40 CONTINUE 236* 237* Reduce number of subdiagonals to KL and number of superdiagonals 238* to KU 239* 240 DO 70 I = 1, MAX( M-1-KL, N-1-KU ) 241 IF( KL.LE.KU ) THEN 242* 243* annihilate subdiagonal elements first (necessary if KL = 0) 244* 245 IF( I.LE.MIN( M-1-KL, N ) ) THEN 246* 247* generate reflection to annihilate A(kl+i+1:m,i) 248* 249 WN = DZNRM2( M-KL-I+1, A( KL+I, I ), 1 ) 250 WA = ( WN / ABS( A( KL+I, I ) ) )*A( KL+I, I ) 251 IF( WN.EQ.ZERO ) THEN 252 TAU = ZERO 253 ELSE 254 WB = A( KL+I, I ) + WA 255 CALL ZSCAL( M-KL-I, ONE / WB, A( KL+I+1, I ), 1 ) 256 A( KL+I, I ) = ONE 257 TAU = DBLE( WB / WA ) 258 END IF 259* 260* apply reflection to A(kl+i:m,i+1:n) from the left 261* 262 CALL ZGEMV( 'Conjugate transpose', M-KL-I+1, N-I, ONE, 263 $ A( KL+I, I+1 ), LDA, A( KL+I, I ), 1, ZERO, 264 $ WORK, 1 ) 265 CALL ZGERC( M-KL-I+1, N-I, -TAU, A( KL+I, I ), 1, WORK, 266 $ 1, A( KL+I, I+1 ), LDA ) 267 A( KL+I, I ) = -WA 268 END IF 269* 270 IF( I.LE.MIN( N-1-KU, M ) ) THEN 271* 272* generate reflection to annihilate A(i,ku+i+1:n) 273* 274 WN = DZNRM2( N-KU-I+1, A( I, KU+I ), LDA ) 275 WA = ( WN / ABS( A( I, KU+I ) ) )*A( I, KU+I ) 276 IF( WN.EQ.ZERO ) THEN 277 TAU = ZERO 278 ELSE 279 WB = A( I, KU+I ) + WA 280 CALL ZSCAL( N-KU-I, ONE / WB, A( I, KU+I+1 ), LDA ) 281 A( I, KU+I ) = ONE 282 TAU = DBLE( WB / WA ) 283 END IF 284* 285* apply reflection to A(i+1:m,ku+i:n) from the right 286* 287 CALL ZLACGV( N-KU-I+1, A( I, KU+I ), LDA ) 288 CALL ZGEMV( 'No transpose', M-I, N-KU-I+1, ONE, 289 $ A( I+1, KU+I ), LDA, A( I, KU+I ), LDA, ZERO, 290 $ WORK, 1 ) 291 CALL ZGERC( M-I, N-KU-I+1, -TAU, WORK, 1, A( I, KU+I ), 292 $ LDA, A( I+1, KU+I ), LDA ) 293 A( I, KU+I ) = -WA 294 END IF 295 ELSE 296* 297* annihilate superdiagonal elements first (necessary if 298* KU = 0) 299* 300 IF( I.LE.MIN( N-1-KU, M ) ) THEN 301* 302* generate reflection to annihilate A(i,ku+i+1:n) 303* 304 WN = DZNRM2( N-KU-I+1, A( I, KU+I ), LDA ) 305 WA = ( WN / ABS( A( I, KU+I ) ) )*A( I, KU+I ) 306 IF( WN.EQ.ZERO ) THEN 307 TAU = ZERO 308 ELSE 309 WB = A( I, KU+I ) + WA 310 CALL ZSCAL( N-KU-I, ONE / WB, A( I, KU+I+1 ), LDA ) 311 A( I, KU+I ) = ONE 312 TAU = DBLE( WB / WA ) 313 END IF 314* 315* apply reflection to A(i+1:m,ku+i:n) from the right 316* 317 CALL ZLACGV( N-KU-I+1, A( I, KU+I ), LDA ) 318 CALL ZGEMV( 'No transpose', M-I, N-KU-I+1, ONE, 319 $ A( I+1, KU+I ), LDA, A( I, KU+I ), LDA, ZERO, 320 $ WORK, 1 ) 321 CALL ZGERC( M-I, N-KU-I+1, -TAU, WORK, 1, A( I, KU+I ), 322 $ LDA, A( I+1, KU+I ), LDA ) 323 A( I, KU+I ) = -WA 324 END IF 325* 326 IF( I.LE.MIN( M-1-KL, N ) ) THEN 327* 328* generate reflection to annihilate A(kl+i+1:m,i) 329* 330 WN = DZNRM2( M-KL-I+1, A( KL+I, I ), 1 ) 331 WA = ( WN / ABS( A( KL+I, I ) ) )*A( KL+I, I ) 332 IF( WN.EQ.ZERO ) THEN 333 TAU = ZERO 334 ELSE 335 WB = A( KL+I, I ) + WA 336 CALL ZSCAL( M-KL-I, ONE / WB, A( KL+I+1, I ), 1 ) 337 A( KL+I, I ) = ONE 338 TAU = DBLE( WB / WA ) 339 END IF 340* 341* apply reflection to A(kl+i:m,i+1:n) from the left 342* 343 CALL ZGEMV( 'Conjugate transpose', M-KL-I+1, N-I, ONE, 344 $ A( KL+I, I+1 ), LDA, A( KL+I, I ), 1, ZERO, 345 $ WORK, 1 ) 346 CALL ZGERC( M-KL-I+1, N-I, -TAU, A( KL+I, I ), 1, WORK, 347 $ 1, A( KL+I, I+1 ), LDA ) 348 A( KL+I, I ) = -WA 349 END IF 350 END IF 351* 352 IF (I .LE. N) THEN 353 DO 50 J = KL + I + 1, M 354 A( J, I ) = ZERO 355 50 CONTINUE 356 END IF 357* 358 IF (I .LE. M) THEN 359 DO 60 J = KU + I + 1, N 360 A( I, J ) = ZERO 361 60 CONTINUE 362 END IF 363 70 CONTINUE 364 RETURN 365* 366* End of ZLAGGE 367* 368 END 369