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