|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLAGGE( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO ) 00002 * 00003 * -- LAPACK auxiliary test routine (version 3.1) 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * November 2006 00006 * 00007 * .. Scalar Arguments .. 00008 INTEGER INFO, KL, KU, LDA, M, N 00009 * .. 00010 * .. Array Arguments .. 00011 INTEGER ISEED( 4 ) 00012 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * ) 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * DLAGGE generates a real general m by n matrix A, by pre- and post- 00019 * multiplying a real diagonal matrix D with random orthogonal matrices: 00020 * A = U*D*V. The lower and upper bandwidths may then be reduced to 00021 * kl and ku by additional orthogonal transformations. 00022 * 00023 * Arguments 00024 * ========= 00025 * 00026 * M (input) INTEGER 00027 * The number of rows of the matrix A. M >= 0. 00028 * 00029 * N (input) INTEGER 00030 * The number of columns of the matrix A. N >= 0. 00031 * 00032 * KL (input) INTEGER 00033 * The number of nonzero subdiagonals within the band of A. 00034 * 0 <= KL <= M-1. 00035 * 00036 * KU (input) INTEGER 00037 * The number of nonzero superdiagonals within the band of A. 00038 * 0 <= KU <= N-1. 00039 * 00040 * D (input) DOUBLE PRECISION array, dimension (min(M,N)) 00041 * The diagonal elements of the diagonal matrix D. 00042 * 00043 * A (output) DOUBLE PRECISION array, dimension (LDA,N) 00044 * The generated m by n matrix A. 00045 * 00046 * LDA (input) INTEGER 00047 * The leading dimension of the array A. LDA >= M. 00048 * 00049 * ISEED (input/output) INTEGER array, dimension (4) 00050 * On entry, the seed of the random number generator; the array 00051 * elements must be between 0 and 4095, and ISEED(4) must be 00052 * odd. 00053 * On exit, the seed is updated. 00054 * 00055 * WORK (workspace) DOUBLE PRECISION array, dimension (M+N) 00056 * 00057 * INFO (output) INTEGER 00058 * = 0: successful exit 00059 * < 0: if INFO = -i, the i-th argument had an illegal value 00060 * 00061 * ===================================================================== 00062 * 00063 * .. Parameters .. 00064 DOUBLE PRECISION ZERO, ONE 00065 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00066 * .. 00067 * .. Local Scalars .. 00068 INTEGER I, J 00069 DOUBLE PRECISION TAU, WA, WB, WN 00070 * .. 00071 * .. External Subroutines .. 00072 EXTERNAL DGEMV, DGER, DLARNV, DSCAL, XERBLA 00073 * .. 00074 * .. Intrinsic Functions .. 00075 INTRINSIC MAX, MIN, SIGN 00076 * .. 00077 * .. External Functions .. 00078 DOUBLE PRECISION DNRM2 00079 EXTERNAL DNRM2 00080 * .. 00081 * .. Executable Statements .. 00082 * 00083 * Test the input arguments 00084 * 00085 INFO = 0 00086 IF( M.LT.0 ) THEN 00087 INFO = -1 00088 ELSE IF( N.LT.0 ) THEN 00089 INFO = -2 00090 ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN 00091 INFO = -3 00092 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 00093 INFO = -4 00094 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00095 INFO = -7 00096 END IF 00097 IF( INFO.LT.0 ) THEN 00098 CALL XERBLA( 'DLAGGE', -INFO ) 00099 RETURN 00100 END IF 00101 * 00102 * initialize A to diagonal matrix 00103 * 00104 DO 20 J = 1, N 00105 DO 10 I = 1, M 00106 A( I, J ) = ZERO 00107 10 CONTINUE 00108 20 CONTINUE 00109 DO 30 I = 1, MIN( M, N ) 00110 A( I, I ) = D( I ) 00111 30 CONTINUE 00112 * 00113 * pre- and post-multiply A by random orthogonal matrices 00114 * 00115 DO 40 I = MIN( M, N ), 1, -1 00116 IF( I.LT.M ) THEN 00117 * 00118 * generate random reflection 00119 * 00120 CALL DLARNV( 3, ISEED, M-I+1, WORK ) 00121 WN = DNRM2( M-I+1, WORK, 1 ) 00122 WA = SIGN( WN, WORK( 1 ) ) 00123 IF( WN.EQ.ZERO ) THEN 00124 TAU = ZERO 00125 ELSE 00126 WB = WORK( 1 ) + WA 00127 CALL DSCAL( M-I, ONE / WB, WORK( 2 ), 1 ) 00128 WORK( 1 ) = ONE 00129 TAU = WB / WA 00130 END IF 00131 * 00132 * multiply A(i:m,i:n) by random reflection from the left 00133 * 00134 CALL DGEMV( 'Transpose', M-I+1, N-I+1, ONE, A( I, I ), LDA, 00135 $ WORK, 1, ZERO, WORK( M+1 ), 1 ) 00136 CALL DGER( M-I+1, N-I+1, -TAU, WORK, 1, WORK( M+1 ), 1, 00137 $ A( I, I ), LDA ) 00138 END IF 00139 IF( I.LT.N ) THEN 00140 * 00141 * generate random reflection 00142 * 00143 CALL DLARNV( 3, ISEED, N-I+1, WORK ) 00144 WN = DNRM2( N-I+1, WORK, 1 ) 00145 WA = SIGN( WN, WORK( 1 ) ) 00146 IF( WN.EQ.ZERO ) THEN 00147 TAU = ZERO 00148 ELSE 00149 WB = WORK( 1 ) + WA 00150 CALL DSCAL( N-I, ONE / WB, WORK( 2 ), 1 ) 00151 WORK( 1 ) = ONE 00152 TAU = WB / WA 00153 END IF 00154 * 00155 * multiply A(i:m,i:n) by random reflection from the right 00156 * 00157 CALL DGEMV( 'No transpose', M-I+1, N-I+1, ONE, A( I, I ), 00158 $ LDA, WORK, 1, ZERO, WORK( N+1 ), 1 ) 00159 CALL DGER( M-I+1, N-I+1, -TAU, WORK( N+1 ), 1, WORK, 1, 00160 $ A( I, I ), LDA ) 00161 END IF 00162 40 CONTINUE 00163 * 00164 * Reduce number of subdiagonals to KL and number of superdiagonals 00165 * to KU 00166 * 00167 DO 70 I = 1, MAX( M-1-KL, N-1-KU ) 00168 IF( KL.LE.KU ) THEN 00169 * 00170 * annihilate subdiagonal elements first (necessary if KL = 0) 00171 * 00172 IF( I.LE.MIN( M-1-KL, N ) ) THEN 00173 * 00174 * generate reflection to annihilate A(kl+i+1:m,i) 00175 * 00176 WN = DNRM2( M-KL-I+1, A( KL+I, I ), 1 ) 00177 WA = SIGN( WN, A( KL+I, I ) ) 00178 IF( WN.EQ.ZERO ) THEN 00179 TAU = ZERO 00180 ELSE 00181 WB = A( KL+I, I ) + WA 00182 CALL DSCAL( M-KL-I, ONE / WB, A( KL+I+1, I ), 1 ) 00183 A( KL+I, I ) = ONE 00184 TAU = WB / WA 00185 END IF 00186 * 00187 * apply reflection to A(kl+i:m,i+1:n) from the left 00188 * 00189 CALL DGEMV( 'Transpose', M-KL-I+1, N-I, ONE, 00190 $ A( KL+I, I+1 ), LDA, A( KL+I, I ), 1, ZERO, 00191 $ WORK, 1 ) 00192 CALL DGER( M-KL-I+1, N-I, -TAU, A( KL+I, I ), 1, WORK, 1, 00193 $ A( KL+I, I+1 ), LDA ) 00194 A( KL+I, I ) = -WA 00195 END IF 00196 * 00197 IF( I.LE.MIN( N-1-KU, M ) ) THEN 00198 * 00199 * generate reflection to annihilate A(i,ku+i+1:n) 00200 * 00201 WN = DNRM2( N-KU-I+1, A( I, KU+I ), LDA ) 00202 WA = SIGN( WN, A( I, KU+I ) ) 00203 IF( WN.EQ.ZERO ) THEN 00204 TAU = ZERO 00205 ELSE 00206 WB = A( I, KU+I ) + WA 00207 CALL DSCAL( N-KU-I, ONE / WB, A( I, KU+I+1 ), LDA ) 00208 A( I, KU+I ) = ONE 00209 TAU = WB / WA 00210 END IF 00211 * 00212 * apply reflection to A(i+1:m,ku+i:n) from the right 00213 * 00214 CALL DGEMV( 'No transpose', M-I, N-KU-I+1, ONE, 00215 $ A( I+1, KU+I ), LDA, A( I, KU+I ), LDA, ZERO, 00216 $ WORK, 1 ) 00217 CALL DGER( M-I, N-KU-I+1, -TAU, WORK, 1, A( I, KU+I ), 00218 $ LDA, A( I+1, KU+I ), LDA ) 00219 A( I, KU+I ) = -WA 00220 END IF 00221 ELSE 00222 * 00223 * annihilate superdiagonal elements first (necessary if 00224 * KU = 0) 00225 * 00226 IF( I.LE.MIN( N-1-KU, M ) ) THEN 00227 * 00228 * generate reflection to annihilate A(i,ku+i+1:n) 00229 * 00230 WN = DNRM2( N-KU-I+1, A( I, KU+I ), LDA ) 00231 WA = SIGN( WN, A( I, KU+I ) ) 00232 IF( WN.EQ.ZERO ) THEN 00233 TAU = ZERO 00234 ELSE 00235 WB = A( I, KU+I ) + WA 00236 CALL DSCAL( N-KU-I, ONE / WB, A( I, KU+I+1 ), LDA ) 00237 A( I, KU+I ) = ONE 00238 TAU = WB / WA 00239 END IF 00240 * 00241 * apply reflection to A(i+1:m,ku+i:n) from the right 00242 * 00243 CALL DGEMV( 'No transpose', M-I, N-KU-I+1, ONE, 00244 $ A( I+1, KU+I ), LDA, A( I, KU+I ), LDA, ZERO, 00245 $ WORK, 1 ) 00246 CALL DGER( M-I, N-KU-I+1, -TAU, WORK, 1, A( I, KU+I ), 00247 $ LDA, A( I+1, KU+I ), LDA ) 00248 A( I, KU+I ) = -WA 00249 END IF 00250 * 00251 IF( I.LE.MIN( M-1-KL, N ) ) THEN 00252 * 00253 * generate reflection to annihilate A(kl+i+1:m,i) 00254 * 00255 WN = DNRM2( M-KL-I+1, A( KL+I, I ), 1 ) 00256 WA = SIGN( WN, A( KL+I, I ) ) 00257 IF( WN.EQ.ZERO ) THEN 00258 TAU = ZERO 00259 ELSE 00260 WB = A( KL+I, I ) + WA 00261 CALL DSCAL( M-KL-I, ONE / WB, A( KL+I+1, I ), 1 ) 00262 A( KL+I, I ) = ONE 00263 TAU = WB / WA 00264 END IF 00265 * 00266 * apply reflection to A(kl+i:m,i+1:n) from the left 00267 * 00268 CALL DGEMV( 'Transpose', M-KL-I+1, N-I, ONE, 00269 $ A( KL+I, I+1 ), LDA, A( KL+I, I ), 1, ZERO, 00270 $ WORK, 1 ) 00271 CALL DGER( M-KL-I+1, N-I, -TAU, A( KL+I, I ), 1, WORK, 1, 00272 $ A( KL+I, I+1 ), LDA ) 00273 A( KL+I, I ) = -WA 00274 END IF 00275 END IF 00276 * 00277 DO 50 J = KL + I + 1, M 00278 A( J, I ) = ZERO 00279 50 CONTINUE 00280 * 00281 DO 60 J = KU + I + 1, N 00282 A( I, J ) = ZERO 00283 60 CONTINUE 00284 70 CONTINUE 00285 RETURN 00286 * 00287 * End of DLAGGE 00288 * 00289 END
1.7.3 
