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