![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLAROR 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 DLAROR( 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 * DOUBLE PRECISION A( LDA, * ), X( * ) 00020 * .. 00021 * 00022 * 00023 *> \par Purpose: 00024 * ============= 00025 *> 00026 *> \verbatim 00027 *> 00028 *> DLAROR pre- or post-multiplies an M by N matrix A by a random 00029 *> orthogonal matrix U, overwriting A. A may optionally be initialized 00030 *> to the identity matrix before multiplying by U. U is generated using 00031 *> the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409). 00032 *> \endverbatim 00033 * 00034 * Arguments: 00035 * ========== 00036 * 00037 *> \param[in] SIDE 00038 *> \verbatim 00039 *> SIDE is CHARACTER*1 00040 *> Specifies whether A is multiplied on the left or right by U. 00041 *> = 'L': Multiply A on the left (premultiply) by U 00042 *> = 'R': Multiply A on the right (postmultiply) by U' 00043 *> = 'C' or 'T': Multiply A on the left by U and the right 00044 *> by U' (Here, U' means U-transpose.) 00045 *> \endverbatim 00046 *> 00047 *> \param[in] INIT 00048 *> \verbatim 00049 *> INIT is CHARACTER*1 00050 *> Specifies whether or not A should be initialized to the 00051 *> identity matrix. 00052 *> = 'I': Initialize A to (a section of) the identity matrix 00053 *> before applying U. 00054 *> = 'N': No initialization. Apply U to the input matrix A. 00055 *> 00056 *> INIT = 'I' may be used to generate square or rectangular 00057 *> orthogonal matrices: 00058 *> 00059 *> For M = N and SIDE = 'L' or 'R', the rows will be orthogonal 00060 *> to each other, as will the columns. 00061 *> 00062 *> If M < N, SIDE = 'R' produces a dense matrix whose rows are 00063 *> orthogonal and whose columns are not, while SIDE = 'L' 00064 *> produces a matrix whose rows are orthogonal, and whose first 00065 *> M columns are orthogonal, and whose remaining columns are 00066 *> zero. 00067 *> 00068 *> If M > N, SIDE = 'L' produces a dense matrix whose columns 00069 *> are orthogonal and whose rows are not, while SIDE = 'R' 00070 *> produces a matrix whose columns are orthogonal, and whose 00071 *> first M rows are orthogonal, and whose remaining rows are 00072 *> zero. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] M 00076 *> \verbatim 00077 *> M is INTEGER 00078 *> The number of rows of A. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] N 00082 *> \verbatim 00083 *> N is INTEGER 00084 *> The number of columns of A. 00085 *> \endverbatim 00086 *> 00087 *> \param[in,out] A 00088 *> \verbatim 00089 *> A is DOUBLE PRECISION array, dimension (LDA, N) 00090 *> On entry, the array A. 00091 *> On exit, overwritten by U A ( if SIDE = 'L' ), 00092 *> or by A U ( if SIDE = 'R' ), 00093 *> or by U A U' ( if SIDE = 'C' or 'T'). 00094 *> \endverbatim 00095 *> 00096 *> \param[in] LDA 00097 *> \verbatim 00098 *> LDA is INTEGER 00099 *> The leading dimension of the array A. LDA >= max(1,M). 00100 *> \endverbatim 00101 *> 00102 *> \param[in,out] ISEED 00103 *> \verbatim 00104 *> ISEED is INTEGER array, dimension (4) 00105 *> On entry ISEED specifies the seed of the random number 00106 *> generator. The array elements should be between 0 and 4095; 00107 *> if not they will be reduced mod 4096. Also, ISEED(4) must 00108 *> be odd. The random number generator uses a linear 00109 *> congruential sequence limited to small integers, and so 00110 *> should produce machine independent random numbers. The 00111 *> values of ISEED are changed on exit, and can be used in the 00112 *> next call to DLAROR to continue the same random number 00113 *> sequence. 00114 *> \endverbatim 00115 *> 00116 *> \param[out] X 00117 *> \verbatim 00118 *> X is DOUBLE PRECISION array, dimension (3*MAX( M, N )) 00119 *> Workspace of length 00120 *> 2*M + N if SIDE = 'L', 00121 *> 2*N + M if SIDE = 'R', 00122 *> 3*N if SIDE = 'C' or 'T'. 00123 *> \endverbatim 00124 *> 00125 *> \param[out] INFO 00126 *> \verbatim 00127 *> INFO is INTEGER 00128 *> An error flag. It is set to: 00129 *> = 0: normal return 00130 *> < 0: if INFO = -k, the k-th argument had an illegal value 00131 *> = 1: if the random numbers generated by DLARND are bad. 00132 *> \endverbatim 00133 * 00134 * Authors: 00135 * ======== 00136 * 00137 *> \author Univ. of Tennessee 00138 *> \author Univ. of California Berkeley 00139 *> \author Univ. of Colorado Denver 00140 *> \author NAG Ltd. 00141 * 00142 *> \date November 2011 00143 * 00144 *> \ingroup double_matgen 00145 * 00146 * ===================================================================== 00147 SUBROUTINE DLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO ) 00148 * 00149 * -- LAPACK auxiliary routine (version 3.4.0) -- 00150 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00151 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00152 * November 2011 00153 * 00154 * .. Scalar Arguments .. 00155 CHARACTER INIT, SIDE 00156 INTEGER INFO, LDA, M, N 00157 * .. 00158 * .. Array Arguments .. 00159 INTEGER ISEED( 4 ) 00160 DOUBLE PRECISION A( LDA, * ), X( * ) 00161 * .. 00162 * 00163 * ===================================================================== 00164 * 00165 * .. Parameters .. 00166 DOUBLE PRECISION ZERO, ONE, TOOSML 00167 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, 00168 $ TOOSML = 1.0D-20 ) 00169 * .. 00170 * .. Local Scalars .. 00171 INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM 00172 DOUBLE PRECISION FACTOR, XNORM, XNORMS 00173 * .. 00174 * .. External Functions .. 00175 LOGICAL LSAME 00176 DOUBLE PRECISION DLARND, DNRM2 00177 EXTERNAL LSAME, DLARND, DNRM2 00178 * .. 00179 * .. External Subroutines .. 00180 EXTERNAL DGEMV, DGER, DLASET, DSCAL, XERBLA 00181 * .. 00182 * .. Intrinsic Functions .. 00183 INTRINSIC ABS, SIGN 00184 * .. 00185 * .. Executable Statements .. 00186 * 00187 INFO = 0 00188 IF( N.EQ.0 .OR. M.EQ.0 ) 00189 $ RETURN 00190 * 00191 ITYPE = 0 00192 IF( LSAME( SIDE, 'L' ) ) THEN 00193 ITYPE = 1 00194 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00195 ITYPE = 2 00196 ELSE IF( LSAME( SIDE, 'C' ) .OR. LSAME( SIDE, 'T' ) ) THEN 00197 ITYPE = 3 00198 END IF 00199 * 00200 * Check for argument errors. 00201 * 00202 IF( ITYPE.EQ.0 ) THEN 00203 INFO = -1 00204 ELSE IF( M.LT.0 ) THEN 00205 INFO = -3 00206 ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN 00207 INFO = -4 00208 ELSE IF( LDA.LT.M ) THEN 00209 INFO = -6 00210 END IF 00211 IF( INFO.NE.0 ) THEN 00212 CALL XERBLA( 'DLAROR', -INFO ) 00213 RETURN 00214 END IF 00215 * 00216 IF( ITYPE.EQ.1 ) THEN 00217 NXFRM = M 00218 ELSE 00219 NXFRM = N 00220 END IF 00221 * 00222 * Initialize A to the identity matrix if desired 00223 * 00224 IF( LSAME( INIT, 'I' ) ) 00225 $ CALL DLASET( 'Full', M, N, ZERO, ONE, A, LDA ) 00226 * 00227 * If no rotation possible, multiply by random +/-1 00228 * 00229 * Compute rotation by computing Householder transformations 00230 * H(2), H(3), ..., H(nhouse) 00231 * 00232 DO 10 J = 1, NXFRM 00233 X( J ) = ZERO 00234 10 CONTINUE 00235 * 00236 DO 30 IXFRM = 2, NXFRM 00237 KBEG = NXFRM - IXFRM + 1 00238 * 00239 * Generate independent normal( 0, 1 ) random numbers 00240 * 00241 DO 20 J = KBEG, NXFRM 00242 X( J ) = DLARND( 3, ISEED ) 00243 20 CONTINUE 00244 * 00245 * Generate a Householder transformation from the random vector X 00246 * 00247 XNORM = DNRM2( IXFRM, X( KBEG ), 1 ) 00248 XNORMS = SIGN( XNORM, X( KBEG ) ) 00249 X( KBEG+NXFRM ) = SIGN( ONE, -X( KBEG ) ) 00250 FACTOR = XNORMS*( XNORMS+X( KBEG ) ) 00251 IF( ABS( FACTOR ).LT.TOOSML ) THEN 00252 INFO = 1 00253 CALL XERBLA( 'DLAROR', INFO ) 00254 RETURN 00255 ELSE 00256 FACTOR = ONE / FACTOR 00257 END IF 00258 X( KBEG ) = X( KBEG ) + XNORMS 00259 * 00260 * Apply Householder transformation to A 00261 * 00262 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN 00263 * 00264 * Apply H(k) from the left. 00265 * 00266 CALL DGEMV( 'T', IXFRM, N, ONE, A( KBEG, 1 ), LDA, 00267 $ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 ) 00268 CALL DGER( IXFRM, N, -FACTOR, X( KBEG ), 1, X( 2*NXFRM+1 ), 00269 $ 1, A( KBEG, 1 ), LDA ) 00270 * 00271 END IF 00272 * 00273 IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN 00274 * 00275 * Apply H(k) from the right. 00276 * 00277 CALL DGEMV( 'N', M, IXFRM, ONE, A( 1, KBEG ), LDA, 00278 $ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 ) 00279 CALL DGER( M, IXFRM, -FACTOR, X( 2*NXFRM+1 ), 1, X( KBEG ), 00280 $ 1, A( 1, KBEG ), LDA ) 00281 * 00282 END IF 00283 30 CONTINUE 00284 * 00285 X( 2*NXFRM ) = SIGN( ONE, DLARND( 3, ISEED ) ) 00286 * 00287 * Scale the matrix A by D. 00288 * 00289 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN 00290 DO 40 IROW = 1, M 00291 CALL DSCAL( N, X( NXFRM+IROW ), A( IROW, 1 ), LDA ) 00292 40 CONTINUE 00293 END IF 00294 * 00295 IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN 00296 DO 50 JCOL = 1, N 00297 CALL DSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 ) 00298 50 CONTINUE 00299 END IF 00300 RETURN 00301 * 00302 * End of DLAROR 00303 * 00304 END