![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLAGGE 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE ZLAGGE( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER INFO, KL, KU, LDA, M, N 00015 * .. 00016 * .. Array Arguments .. 00017 * INTEGER ISEED( 4 ) 00018 * DOUBLE PRECISION D( * ) 00019 * COMPLEX*16 A( LDA, * ), WORK( * ) 00020 * .. 00021 * 00022 * 00023 *> \par Purpose: 00024 * ============= 00025 *> 00026 *> \verbatim 00027 *> 00028 *> ZLAGGE generates a complex general m by n matrix A, by pre- and post- 00029 *> multiplying a real diagonal matrix D with random unitary matrices: 00030 *> A = U*D*V. The lower and upper bandwidths may then be reduced to 00031 *> kl and ku by additional unitary transformations. 00032 *> \endverbatim 00033 * 00034 * Arguments: 00035 * ========== 00036 * 00037 *> \param[in] M 00038 *> \verbatim 00039 *> M is INTEGER 00040 *> The number of rows of the matrix A. M >= 0. 00041 *> \endverbatim 00042 *> 00043 *> \param[in] N 00044 *> \verbatim 00045 *> N is INTEGER 00046 *> The number of columns of the matrix A. N >= 0. 00047 *> \endverbatim 00048 *> 00049 *> \param[in] KL 00050 *> \verbatim 00051 *> KL is INTEGER 00052 *> The number of nonzero subdiagonals within the band of A. 00053 *> 0 <= KL <= M-1. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] KU 00057 *> \verbatim 00058 *> KU is INTEGER 00059 *> The number of nonzero superdiagonals within the band of A. 00060 *> 0 <= KU <= N-1. 00061 *> \endverbatim 00062 *> 00063 *> \param[in] D 00064 *> \verbatim 00065 *> D is DOUBLE PRECISION array, dimension (min(M,N)) 00066 *> The diagonal elements of the diagonal matrix D. 00067 *> \endverbatim 00068 *> 00069 *> \param[out] A 00070 *> \verbatim 00071 *> A is COMPLEX*16 array, dimension (LDA,N) 00072 *> The generated m by n matrix A. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] LDA 00076 *> \verbatim 00077 *> LDA is INTEGER 00078 *> The leading dimension of the array A. LDA >= M. 00079 *> \endverbatim 00080 *> 00081 *> \param[in,out] ISEED 00082 *> \verbatim 00083 *> ISEED is INTEGER array, dimension (4) 00084 *> On entry, the seed of the random number generator; the array 00085 *> elements must be between 0 and 4095, and ISEED(4) must be 00086 *> odd. 00087 *> On exit, the seed is updated. 00088 *> \endverbatim 00089 *> 00090 *> \param[out] WORK 00091 *> \verbatim 00092 *> WORK is COMPLEX*16 array, dimension (M+N) 00093 *> \endverbatim 00094 *> 00095 *> \param[out] INFO 00096 *> \verbatim 00097 *> INFO is INTEGER 00098 *> = 0: successful exit 00099 *> < 0: if INFO = -i, the i-th argument had an illegal value 00100 *> \endverbatim 00101 * 00102 * Authors: 00103 * ======== 00104 * 00105 *> \author Univ. of Tennessee 00106 *> \author Univ. of California Berkeley 00107 *> \author Univ. of Colorado Denver 00108 *> \author NAG Ltd. 00109 * 00110 *> \date November 2011 00111 * 00112 *> \ingroup complex16_matgen 00113 * 00114 * ===================================================================== 00115 SUBROUTINE ZLAGGE( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO ) 00116 * 00117 * -- LAPACK auxiliary routine (version 3.4.0) -- 00118 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00120 * November 2011 00121 * 00122 * .. Scalar Arguments .. 00123 INTEGER INFO, KL, KU, LDA, M, N 00124 * .. 00125 * .. Array Arguments .. 00126 INTEGER ISEED( 4 ) 00127 DOUBLE PRECISION D( * ) 00128 COMPLEX*16 A( LDA, * ), WORK( * ) 00129 * .. 00130 * 00131 * ===================================================================== 00132 * 00133 * .. Parameters .. 00134 COMPLEX*16 ZERO, ONE 00135 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 00136 $ ONE = ( 1.0D+0, 0.0D+0 ) ) 00137 * .. 00138 * .. Local Scalars .. 00139 INTEGER I, J 00140 DOUBLE PRECISION WN 00141 COMPLEX*16 TAU, WA, WB 00142 * .. 00143 * .. External Subroutines .. 00144 EXTERNAL XERBLA, ZGEMV, ZGERC, ZLACGV, ZLARNV, ZSCAL 00145 * .. 00146 * .. Intrinsic Functions .. 00147 INTRINSIC ABS, DBLE, MAX, MIN 00148 * .. 00149 * .. External Functions .. 00150 DOUBLE PRECISION DZNRM2 00151 EXTERNAL DZNRM2 00152 * .. 00153 * .. Executable Statements .. 00154 * 00155 * Test the input arguments 00156 * 00157 INFO = 0 00158 IF( M.LT.0 ) THEN 00159 INFO = -1 00160 ELSE IF( N.LT.0 ) THEN 00161 INFO = -2 00162 ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN 00163 INFO = -3 00164 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 00165 INFO = -4 00166 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00167 INFO = -7 00168 END IF 00169 IF( INFO.LT.0 ) THEN 00170 CALL XERBLA( 'ZLAGGE', -INFO ) 00171 RETURN 00172 END IF 00173 * 00174 * initialize A to diagonal matrix 00175 * 00176 DO 20 J = 1, N 00177 DO 10 I = 1, M 00178 A( I, J ) = ZERO 00179 10 CONTINUE 00180 20 CONTINUE 00181 DO 30 I = 1, MIN( M, N ) 00182 A( I, I ) = D( I ) 00183 30 CONTINUE 00184 * 00185 * pre- and post-multiply A by random unitary matrices 00186 * 00187 DO 40 I = MIN( M, N ), 1, -1 00188 IF( I.LT.M ) THEN 00189 * 00190 * generate random reflection 00191 * 00192 CALL ZLARNV( 3, ISEED, M-I+1, WORK ) 00193 WN = DZNRM2( M-I+1, WORK, 1 ) 00194 WA = ( WN / ABS( WORK( 1 ) ) )*WORK( 1 ) 00195 IF( WN.EQ.ZERO ) THEN 00196 TAU = ZERO 00197 ELSE 00198 WB = WORK( 1 ) + WA 00199 CALL ZSCAL( M-I, ONE / WB, WORK( 2 ), 1 ) 00200 WORK( 1 ) = ONE 00201 TAU = DBLE( WB / WA ) 00202 END IF 00203 * 00204 * multiply A(i:m,i:n) by random reflection from the left 00205 * 00206 CALL ZGEMV( 'Conjugate transpose', M-I+1, N-I+1, ONE, 00207 $ A( I, I ), LDA, WORK, 1, ZERO, WORK( M+1 ), 1 ) 00208 CALL ZGERC( M-I+1, N-I+1, -TAU, WORK, 1, WORK( M+1 ), 1, 00209 $ A( I, I ), LDA ) 00210 END IF 00211 IF( I.LT.N ) THEN 00212 * 00213 * generate random reflection 00214 * 00215 CALL ZLARNV( 3, ISEED, N-I+1, WORK ) 00216 WN = DZNRM2( N-I+1, WORK, 1 ) 00217 WA = ( WN / ABS( WORK( 1 ) ) )*WORK( 1 ) 00218 IF( WN.EQ.ZERO ) THEN 00219 TAU = ZERO 00220 ELSE 00221 WB = WORK( 1 ) + WA 00222 CALL ZSCAL( N-I, ONE / WB, WORK( 2 ), 1 ) 00223 WORK( 1 ) = ONE 00224 TAU = DBLE( WB / WA ) 00225 END IF 00226 * 00227 * multiply A(i:m,i:n) by random reflection from the right 00228 * 00229 CALL ZGEMV( 'No transpose', M-I+1, N-I+1, ONE, A( I, I ), 00230 $ LDA, WORK, 1, ZERO, WORK( N+1 ), 1 ) 00231 CALL ZGERC( M-I+1, N-I+1, -TAU, WORK( N+1 ), 1, WORK, 1, 00232 $ A( I, I ), LDA ) 00233 END IF 00234 40 CONTINUE 00235 * 00236 * Reduce number of subdiagonals to KL and number of superdiagonals 00237 * to KU 00238 * 00239 DO 70 I = 1, MAX( M-1-KL, N-1-KU ) 00240 IF( KL.LE.KU ) THEN 00241 * 00242 * annihilate subdiagonal elements first (necessary if KL = 0) 00243 * 00244 IF( I.LE.MIN( M-1-KL, N ) ) THEN 00245 * 00246 * generate reflection to annihilate A(kl+i+1:m,i) 00247 * 00248 WN = DZNRM2( M-KL-I+1, A( KL+I, I ), 1 ) 00249 WA = ( WN / ABS( A( KL+I, I ) ) )*A( KL+I, I ) 00250 IF( WN.EQ.ZERO ) THEN 00251 TAU = ZERO 00252 ELSE 00253 WB = A( KL+I, I ) + WA 00254 CALL ZSCAL( M-KL-I, ONE / WB, A( KL+I+1, I ), 1 ) 00255 A( KL+I, I ) = ONE 00256 TAU = DBLE( WB / WA ) 00257 END IF 00258 * 00259 * apply reflection to A(kl+i:m,i+1:n) from the left 00260 * 00261 CALL ZGEMV( 'Conjugate transpose', M-KL-I+1, N-I, ONE, 00262 $ A( KL+I, I+1 ), LDA, A( KL+I, I ), 1, ZERO, 00263 $ WORK, 1 ) 00264 CALL ZGERC( M-KL-I+1, N-I, -TAU, A( KL+I, I ), 1, WORK, 00265 $ 1, A( KL+I, I+1 ), LDA ) 00266 A( KL+I, I ) = -WA 00267 END IF 00268 * 00269 IF( I.LE.MIN( N-1-KU, M ) ) THEN 00270 * 00271 * generate reflection to annihilate A(i,ku+i+1:n) 00272 * 00273 WN = DZNRM2( N-KU-I+1, A( I, KU+I ), LDA ) 00274 WA = ( WN / ABS( A( I, KU+I ) ) )*A( I, KU+I ) 00275 IF( WN.EQ.ZERO ) THEN 00276 TAU = ZERO 00277 ELSE 00278 WB = A( I, KU+I ) + WA 00279 CALL ZSCAL( N-KU-I, ONE / WB, A( I, KU+I+1 ), LDA ) 00280 A( I, KU+I ) = ONE 00281 TAU = DBLE( WB / WA ) 00282 END IF 00283 * 00284 * apply reflection to A(i+1:m,ku+i:n) from the right 00285 * 00286 CALL ZLACGV( N-KU-I+1, A( I, KU+I ), LDA ) 00287 CALL ZGEMV( 'No transpose', M-I, N-KU-I+1, ONE, 00288 $ A( I+1, KU+I ), LDA, A( I, KU+I ), LDA, ZERO, 00289 $ WORK, 1 ) 00290 CALL ZGERC( M-I, N-KU-I+1, -TAU, WORK, 1, A( I, KU+I ), 00291 $ LDA, A( I+1, KU+I ), LDA ) 00292 A( I, KU+I ) = -WA 00293 END IF 00294 ELSE 00295 * 00296 * annihilate superdiagonal elements first (necessary if 00297 * KU = 0) 00298 * 00299 IF( I.LE.MIN( N-1-KU, M ) ) THEN 00300 * 00301 * generate reflection to annihilate A(i,ku+i+1:n) 00302 * 00303 WN = DZNRM2( N-KU-I+1, A( I, KU+I ), LDA ) 00304 WA = ( WN / ABS( A( I, KU+I ) ) )*A( I, KU+I ) 00305 IF( WN.EQ.ZERO ) THEN 00306 TAU = ZERO 00307 ELSE 00308 WB = A( I, KU+I ) + WA 00309 CALL ZSCAL( N-KU-I, ONE / WB, A( I, KU+I+1 ), LDA ) 00310 A( I, KU+I ) = ONE 00311 TAU = DBLE( WB / WA ) 00312 END IF 00313 * 00314 * apply reflection to A(i+1:m,ku+i:n) from the right 00315 * 00316 CALL ZLACGV( N-KU-I+1, A( I, KU+I ), LDA ) 00317 CALL ZGEMV( 'No transpose', M-I, N-KU-I+1, ONE, 00318 $ A( I+1, KU+I ), LDA, A( I, KU+I ), LDA, ZERO, 00319 $ WORK, 1 ) 00320 CALL ZGERC( M-I, N-KU-I+1, -TAU, WORK, 1, A( I, KU+I ), 00321 $ LDA, A( I+1, KU+I ), LDA ) 00322 A( I, KU+I ) = -WA 00323 END IF 00324 * 00325 IF( I.LE.MIN( M-1-KL, N ) ) THEN 00326 * 00327 * generate reflection to annihilate A(kl+i+1:m,i) 00328 * 00329 WN = DZNRM2( M-KL-I+1, A( KL+I, I ), 1 ) 00330 WA = ( WN / ABS( A( KL+I, I ) ) )*A( KL+I, I ) 00331 IF( WN.EQ.ZERO ) THEN 00332 TAU = ZERO 00333 ELSE 00334 WB = A( KL+I, I ) + WA 00335 CALL ZSCAL( M-KL-I, ONE / WB, A( KL+I+1, I ), 1 ) 00336 A( KL+I, I ) = ONE 00337 TAU = DBLE( WB / WA ) 00338 END IF 00339 * 00340 * apply reflection to A(kl+i:m,i+1:n) from the left 00341 * 00342 CALL ZGEMV( 'Conjugate transpose', M-KL-I+1, N-I, ONE, 00343 $ A( KL+I, I+1 ), LDA, A( KL+I, I ), 1, ZERO, 00344 $ WORK, 1 ) 00345 CALL ZGERC( M-KL-I+1, N-I, -TAU, A( KL+I, I ), 1, WORK, 00346 $ 1, A( KL+I, I+1 ), LDA ) 00347 A( KL+I, I ) = -WA 00348 END IF 00349 END IF 00350 * 00351 DO 50 J = KL + I + 1, M 00352 A( J, I ) = ZERO 00353 50 CONTINUE 00354 * 00355 DO 60 J = KU + I + 1, N 00356 A( I, J ) = ZERO 00357 60 CONTINUE 00358 70 CONTINUE 00359 RETURN 00360 * 00361 * End of ZLAGGE 00362 * 00363 END