![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLAROR 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 ZLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO ) 00012 * 00013 * .. Scalar Arguments .. 00014 * CHARACTER INIT, SIDE 00015 * INTEGER INFO, LDA, M, N 00016 * .. 00017 * .. Array Arguments .. 00018 * INTEGER ISEED( 4 ) 00019 * COMPLEX*16 A( LDA, * ), X( * ) 00020 * .. 00021 * 00022 * 00023 *> \par Purpose: 00024 * ============= 00025 *> 00026 *> \verbatim 00027 *> 00028 *> ZLAROR pre- or post-multiplies an M by N matrix A by a random 00029 *> unitary matrix U, overwriting A. A may optionally be 00030 *> initialized to the identity matrix before multiplying by U. 00031 *> U is generated using the method of G.W. Stewart 00032 *> ( SIAM J. Numer. Anal. 17, 1980, pp. 403-409 ). 00033 *> (BLAS-2 version) 00034 *> \endverbatim 00035 * 00036 * Arguments: 00037 * ========== 00038 * 00039 *> \param[in] SIDE 00040 *> \verbatim 00041 *> SIDE is CHARACTER*1 00042 *> SIDE specifies whether A is multiplied on the left or right 00043 *> by U. 00044 *> SIDE = 'L' Multiply A on the left (premultiply) by U 00045 *> SIDE = 'R' Multiply A on the right (postmultiply) by UC> SIDE = 'C' Multiply A on the left by U and the right by UC> SIDE = 'T' Multiply A on the left by U and the right by U' 00046 *> Not modified. 00047 *> \endverbatim 00048 *> 00049 *> \param[in] INIT 00050 *> \verbatim 00051 *> INIT is CHARACTER*1 00052 *> INIT specifies whether or not A should be initialized to 00053 *> the identity matrix. 00054 *> INIT = 'I' Initialize A to (a section of) the 00055 *> identity matrix before applying U. 00056 *> INIT = 'N' No initialization. Apply U to the 00057 *> input matrix A. 00058 *> 00059 *> INIT = 'I' may be used to generate square (i.e., unitary) 00060 *> or rectangular orthogonal matrices (orthogonality being 00061 *> in the sense of ZDOTC): 00062 *> 00063 *> For square matrices, M=N, and SIDE many be either 'L' or 00064 *> 'R'; the rows will be orthogonal to each other, as will the 00065 *> columns. 00066 *> For rectangular matrices where M < N, SIDE = 'R' will 00067 *> produce a dense matrix whose rows will be orthogonal and 00068 *> whose columns will not, while SIDE = 'L' will produce a 00069 *> matrix whose rows will be orthogonal, and whose first M 00070 *> columns will be orthogonal, the remaining columns being 00071 *> zero. 00072 *> For matrices where M > N, just use the previous 00073 *> explanation, interchanging 'L' and 'R' and "rows" and 00074 *> "columns". 00075 *> 00076 *> Not modified. 00077 *> \endverbatim 00078 *> 00079 *> \param[in] M 00080 *> \verbatim 00081 *> M is INTEGER 00082 *> Number of rows of A. Not modified. 00083 *> \endverbatim 00084 *> 00085 *> \param[in] N 00086 *> \verbatim 00087 *> N is INTEGER 00088 *> Number of columns of A. Not modified. 00089 *> \endverbatim 00090 *> 00091 *> \param[in,out] A 00092 *> \verbatim 00093 *> A is COMPLEX*16 array, dimension ( LDA, N ) 00094 *> Input and output array. Overwritten by U A ( if SIDE = 'L' ) 00095 *> or by A U ( if SIDE = 'R' ) 00096 *> or by U A U* ( if SIDE = 'C') 00097 *> or by U A U' ( if SIDE = 'T') on exit. 00098 *> \endverbatim 00099 *> 00100 *> \param[in] LDA 00101 *> \verbatim 00102 *> LDA is INTEGER 00103 *> Leading dimension of A. Must be at least MAX ( 1, M ). 00104 *> Not modified. 00105 *> \endverbatim 00106 *> 00107 *> \param[in,out] ISEED 00108 *> \verbatim 00109 *> ISEED is INTEGER array, dimension ( 4 ) 00110 *> On entry ISEED specifies the seed of the random number 00111 *> generator. The array elements should be between 0 and 4095; 00112 *> if not they will be reduced mod 4096. Also, ISEED(4) must 00113 *> be odd. The random number generator uses a linear 00114 *> congruential sequence limited to small integers, and so 00115 *> should produce machine independent random numbers. The 00116 *> values of ISEED are changed on exit, and can be used in the 00117 *> next call to ZLAROR to continue the same random number 00118 *> sequence. 00119 *> Modified. 00120 *> \endverbatim 00121 *> 00122 *> \param[out] X 00123 *> \verbatim 00124 *> X is COMPLEX*16 array, dimension ( 3*MAX( M, N ) ) 00125 *> Workspace. Of length: 00126 *> 2*M + N if SIDE = 'L', 00127 *> 2*N + M if SIDE = 'R', 00128 *> 3*N if SIDE = 'C' or 'T'. 00129 *> Modified. 00130 *> \endverbatim 00131 *> 00132 *> \param[out] INFO 00133 *> \verbatim 00134 *> INFO is INTEGER 00135 *> An error flag. It is set to: 00136 *> 0 if no error. 00137 *> 1 if ZLARND returned a bad random number (installation 00138 *> problem) 00139 *> -1 if SIDE is not L, R, C, or T. 00140 *> -3 if M is negative. 00141 *> -4 if N is negative or if SIDE is C or T and N is not equal 00142 *> to M. 00143 *> -6 if LDA is less than M. 00144 *> \endverbatim 00145 * 00146 * Authors: 00147 * ======== 00148 * 00149 *> \author Univ. of Tennessee 00150 *> \author Univ. of California Berkeley 00151 *> \author Univ. of Colorado Denver 00152 *> \author NAG Ltd. 00153 * 00154 *> \date November 2011 00155 * 00156 *> \ingroup complex16_matgen 00157 * 00158 * ===================================================================== 00159 SUBROUTINE ZLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO ) 00160 * 00161 * -- LAPACK auxiliary routine (version 3.4.0) -- 00162 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00163 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00164 * November 2011 00165 * 00166 * .. Scalar Arguments .. 00167 CHARACTER INIT, SIDE 00168 INTEGER INFO, LDA, M, N 00169 * .. 00170 * .. Array Arguments .. 00171 INTEGER ISEED( 4 ) 00172 COMPLEX*16 A( LDA, * ), X( * ) 00173 * .. 00174 * 00175 * ===================================================================== 00176 * 00177 * .. Parameters .. 00178 DOUBLE PRECISION ZERO, ONE, TOOSML 00179 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, 00180 $ TOOSML = 1.0D-20 ) 00181 COMPLEX*16 CZERO, CONE 00182 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00183 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00184 * .. 00185 * .. Local Scalars .. 00186 INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM 00187 DOUBLE PRECISION FACTOR, XABS, XNORM 00188 COMPLEX*16 CSIGN, XNORMS 00189 * .. 00190 * .. External Functions .. 00191 LOGICAL LSAME 00192 DOUBLE PRECISION DZNRM2 00193 COMPLEX*16 ZLARND 00194 EXTERNAL LSAME, DZNRM2, ZLARND 00195 * .. 00196 * .. External Subroutines .. 00197 EXTERNAL XERBLA, ZGEMV, ZGERC, ZLACGV, ZLASET, ZSCAL 00198 * .. 00199 * .. Intrinsic Functions .. 00200 INTRINSIC ABS, DCMPLX, DCONJG 00201 * .. 00202 * .. Executable Statements .. 00203 * 00204 INFO = 0 00205 IF( N.EQ.0 .OR. M.EQ.0 ) 00206 $ RETURN 00207 * 00208 ITYPE = 0 00209 IF( LSAME( SIDE, 'L' ) ) THEN 00210 ITYPE = 1 00211 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00212 ITYPE = 2 00213 ELSE IF( LSAME( SIDE, 'C' ) ) THEN 00214 ITYPE = 3 00215 ELSE IF( LSAME( SIDE, 'T' ) ) THEN 00216 ITYPE = 4 00217 END IF 00218 * 00219 * Check for argument errors. 00220 * 00221 IF( ITYPE.EQ.0 ) THEN 00222 INFO = -1 00223 ELSE IF( M.LT.0 ) THEN 00224 INFO = -3 00225 ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN 00226 INFO = -4 00227 ELSE IF( LDA.LT.M ) THEN 00228 INFO = -6 00229 END IF 00230 IF( INFO.NE.0 ) THEN 00231 CALL XERBLA( 'ZLAROR', -INFO ) 00232 RETURN 00233 END IF 00234 * 00235 IF( ITYPE.EQ.1 ) THEN 00236 NXFRM = M 00237 ELSE 00238 NXFRM = N 00239 END IF 00240 * 00241 * Initialize A to the identity matrix if desired 00242 * 00243 IF( LSAME( INIT, 'I' ) ) 00244 $ CALL ZLASET( 'Full', M, N, CZERO, CONE, A, LDA ) 00245 * 00246 * If no rotation possible, still multiply by 00247 * a random complex number from the circle |x| = 1 00248 * 00249 * 2) Compute Rotation by computing Householder 00250 * Transformations H(2), H(3), ..., H(n). Note that the 00251 * order in which they are computed is irrelevant. 00252 * 00253 DO 10 J = 1, NXFRM 00254 X( J ) = CZERO 00255 10 CONTINUE 00256 * 00257 DO 30 IXFRM = 2, NXFRM 00258 KBEG = NXFRM - IXFRM + 1 00259 * 00260 * Generate independent normal( 0, 1 ) random numbers 00261 * 00262 DO 20 J = KBEG, NXFRM 00263 X( J ) = ZLARND( 3, ISEED ) 00264 20 CONTINUE 00265 * 00266 * Generate a Householder transformation from the random vector X 00267 * 00268 XNORM = DZNRM2( IXFRM, X( KBEG ), 1 ) 00269 XABS = ABS( X( KBEG ) ) 00270 IF( XABS.NE.CZERO ) THEN 00271 CSIGN = X( KBEG ) / XABS 00272 ELSE 00273 CSIGN = CONE 00274 END IF 00275 XNORMS = CSIGN*XNORM 00276 X( NXFRM+KBEG ) = -CSIGN 00277 FACTOR = XNORM*( XNORM+XABS ) 00278 IF( ABS( FACTOR ).LT.TOOSML ) THEN 00279 INFO = 1 00280 CALL XERBLA( 'ZLAROR', -INFO ) 00281 RETURN 00282 ELSE 00283 FACTOR = ONE / FACTOR 00284 END IF 00285 X( KBEG ) = X( KBEG ) + XNORMS 00286 * 00287 * Apply Householder transformation to A 00288 * 00289 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 .OR. ITYPE.EQ.4 ) THEN 00290 * 00291 * Apply H(k) on the left of A 00292 * 00293 CALL ZGEMV( 'C', IXFRM, N, CONE, A( KBEG, 1 ), LDA, 00294 $ X( KBEG ), 1, CZERO, X( 2*NXFRM+1 ), 1 ) 00295 CALL ZGERC( IXFRM, N, -DCMPLX( FACTOR ), X( KBEG ), 1, 00296 $ X( 2*NXFRM+1 ), 1, A( KBEG, 1 ), LDA ) 00297 * 00298 END IF 00299 * 00300 IF( ITYPE.GE.2 .AND. ITYPE.LE.4 ) THEN 00301 * 00302 * Apply H(k)* (or H(k)') on the right of A 00303 * 00304 IF( ITYPE.EQ.4 ) THEN 00305 CALL ZLACGV( IXFRM, X( KBEG ), 1 ) 00306 END IF 00307 * 00308 CALL ZGEMV( 'N', M, IXFRM, CONE, A( 1, KBEG ), LDA, 00309 $ X( KBEG ), 1, CZERO, X( 2*NXFRM+1 ), 1 ) 00310 CALL ZGERC( M, IXFRM, -DCMPLX( FACTOR ), X( 2*NXFRM+1 ), 1, 00311 $ X( KBEG ), 1, A( 1, KBEG ), LDA ) 00312 * 00313 END IF 00314 30 CONTINUE 00315 * 00316 X( 1 ) = ZLARND( 3, ISEED ) 00317 XABS = ABS( X( 1 ) ) 00318 IF( XABS.NE.ZERO ) THEN 00319 CSIGN = X( 1 ) / XABS 00320 ELSE 00321 CSIGN = CONE 00322 END IF 00323 X( 2*NXFRM ) = CSIGN 00324 * 00325 * Scale the matrix A by D. 00326 * 00327 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 .OR. ITYPE.EQ.4 ) THEN 00328 DO 40 IROW = 1, M 00329 CALL ZSCAL( N, DCONJG( X( NXFRM+IROW ) ), A( IROW, 1 ), 00330 $ LDA ) 00331 40 CONTINUE 00332 END IF 00333 * 00334 IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN 00335 DO 50 JCOL = 1, N 00336 CALL ZSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 ) 00337 50 CONTINUE 00338 END IF 00339 * 00340 IF( ITYPE.EQ.4 ) THEN 00341 DO 60 JCOL = 1, N 00342 CALL ZSCAL( M, DCONJG( X( NXFRM+JCOL ) ), A( 1, JCOL ), 1 ) 00343 60 CONTINUE 00344 END IF 00345 RETURN 00346 * 00347 * End of ZLAROR 00348 * 00349 END