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