 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLASR( SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER DIRECT, PIVOT, SIDE 00010 INTEGER LDA, M, N 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION A( LDA, * ), C( * ), S( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DLASR applies a sequence of plane rotations to a real matrix A, 00020 * from either the left or the right. 00021 * 00022 * When SIDE = 'L', the transformation takes the form 00023 * 00024 * A := P*A 00025 * 00026 * and when SIDE = 'R', the transformation takes the form 00027 * 00028 * A := A*P**T 00029 * 00030 * where P is an orthogonal matrix consisting of a sequence of z plane 00031 * rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R', 00032 * and P**T is the transpose of P. 00033 * 00034 * When DIRECT = 'F' (Forward sequence), then 00035 * 00036 * P = P(z-1) * ... * P(2) * P(1) 00037 * 00038 * and when DIRECT = 'B' (Backward sequence), then 00039 * 00040 * P = P(1) * P(2) * ... * P(z-1) 00041 * 00042 * where P(k) is a plane rotation matrix defined by the 2-by-2 rotation 00043 * 00044 * R(k) = ( c(k) s(k) ) 00045 * = ( -s(k) c(k) ). 00046 * 00047 * When PIVOT = 'V' (Variable pivot), the rotation is performed 00048 * for the plane (k,k+1), i.e., P(k) has the form 00049 * 00050 * P(k) = ( 1 ) 00051 * ( ... ) 00052 * ( 1 ) 00053 * ( c(k) s(k) ) 00054 * ( -s(k) c(k) ) 00055 * ( 1 ) 00056 * ( ... ) 00057 * ( 1 ) 00058 * 00059 * where R(k) appears as a rank-2 modification to the identity matrix in 00060 * rows and columns k and k+1. 00061 * 00062 * When PIVOT = 'T' (Top pivot), the rotation is performed for the 00063 * plane (1,k+1), so P(k) has the form 00064 * 00065 * P(k) = ( c(k) s(k) ) 00066 * ( 1 ) 00067 * ( ... ) 00068 * ( 1 ) 00069 * ( -s(k) c(k) ) 00070 * ( 1 ) 00071 * ( ... ) 00072 * ( 1 ) 00073 * 00074 * where R(k) appears in rows and columns 1 and k+1. 00075 * 00076 * Similarly, when PIVOT = 'B' (Bottom pivot), the rotation is 00077 * performed for the plane (k,z), giving P(k) the form 00078 * 00079 * P(k) = ( 1 ) 00080 * ( ... ) 00081 * ( 1 ) 00082 * ( c(k) s(k) ) 00083 * ( 1 ) 00084 * ( ... ) 00085 * ( 1 ) 00086 * ( -s(k) c(k) ) 00087 * 00088 * where R(k) appears in rows and columns k and z. The rotations are 00089 * performed without ever forming P(k) explicitly. 00090 * 00091 * Arguments 00092 * ========= 00093 * 00094 * SIDE (input) CHARACTER*1 00095 * Specifies whether the plane rotation matrix P is applied to 00096 * A on the left or the right. 00097 * = 'L': Left, compute A := P*A 00098 * = 'R': Right, compute A:= A*P**T 00099 * 00100 * PIVOT (input) CHARACTER*1 00101 * Specifies the plane for which P(k) is a plane rotation 00102 * matrix. 00103 * = 'V': Variable pivot, the plane (k,k+1) 00104 * = 'T': Top pivot, the plane (1,k+1) 00105 * = 'B': Bottom pivot, the plane (k,z) 00106 * 00107 * DIRECT (input) CHARACTER*1 00108 * Specifies whether P is a forward or backward sequence of 00109 * plane rotations. 00110 * = 'F': Forward, P = P(z-1)*...*P(2)*P(1) 00111 * = 'B': Backward, P = P(1)*P(2)*...*P(z-1) 00112 * 00113 * M (input) INTEGER 00114 * The number of rows of the matrix A. If m <= 1, an immediate 00115 * return is effected. 00116 * 00117 * N (input) INTEGER 00118 * The number of columns of the matrix A. If n <= 1, an 00119 * immediate return is effected. 00120 * 00121 * C (input) DOUBLE PRECISION array, dimension 00122 * (M-1) if SIDE = 'L' 00123 * (N-1) if SIDE = 'R' 00124 * The cosines c(k) of the plane rotations. 00125 * 00126 * S (input) DOUBLE PRECISION array, dimension 00127 * (M-1) if SIDE = 'L' 00128 * (N-1) if SIDE = 'R' 00129 * The sines s(k) of the plane rotations. The 2-by-2 plane 00130 * rotation part of the matrix P(k), R(k), has the form 00131 * R(k) = ( c(k) s(k) ) 00132 * ( -s(k) c(k) ). 00133 * 00134 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00135 * The M-by-N matrix A. On exit, A is overwritten by P*A if 00136 * SIDE = 'R' or by A*P**T if SIDE = 'L'. 00137 * 00138 * LDA (input) INTEGER 00139 * The leading dimension of the array A. LDA >= max(1,M). 00140 * 00141 * ===================================================================== 00142 * 00143 * .. Parameters .. 00144 DOUBLE PRECISION ONE, ZERO 00145 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00146 * .. 00147 * .. Local Scalars .. 00148 INTEGER I, INFO, J 00149 DOUBLE PRECISION CTEMP, STEMP, TEMP 00150 * .. 00151 * .. External Functions .. 00152 LOGICAL LSAME 00153 EXTERNAL LSAME 00154 * .. 00155 * .. External Subroutines .. 00156 EXTERNAL XERBLA 00157 * .. 00158 * .. Intrinsic Functions .. 00159 INTRINSIC MAX 00160 * .. 00161 * .. Executable Statements .. 00162 * 00163 * Test the input parameters 00164 * 00165 INFO = 0 00166 IF( .NOT.( LSAME( SIDE, 'L' ) .OR. LSAME( SIDE, 'R' ) ) ) THEN 00167 INFO = 1 00168 ELSE IF( .NOT.( LSAME( PIVOT, 'V' ) .OR. LSAME( PIVOT, 00169 $ 'T' ) .OR. LSAME( PIVOT, 'B' ) ) ) THEN 00170 INFO = 2 00171 ELSE IF( .NOT.( LSAME( DIRECT, 'F' ) .OR. LSAME( DIRECT, 'B' ) ) ) 00172 $ THEN 00173 INFO = 3 00174 ELSE IF( M.LT.0 ) THEN 00175 INFO = 4 00176 ELSE IF( N.LT.0 ) THEN 00177 INFO = 5 00178 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00179 INFO = 9 00180 END IF 00181 IF( INFO.NE.0 ) THEN 00182 CALL XERBLA( 'DLASR ', INFO ) 00183 RETURN 00184 END IF 00185 * 00186 * Quick return if possible 00187 * 00188 IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) 00189 $ RETURN 00190 IF( LSAME( SIDE, 'L' ) ) THEN 00191 * 00192 * Form P * A 00193 * 00194 IF( LSAME( PIVOT, 'V' ) ) THEN 00195 IF( LSAME( DIRECT, 'F' ) ) THEN 00196 DO 20 J = 1, M - 1 00197 CTEMP = C( J ) 00198 STEMP = S( J ) 00199 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00200 DO 10 I = 1, N 00201 TEMP = A( J+1, I ) 00202 A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I ) 00203 A( J, I ) = STEMP*TEMP + CTEMP*A( J, I ) 00204 10 CONTINUE 00205 END IF 00206 20 CONTINUE 00207 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 00208 DO 40 J = M - 1, 1, -1 00209 CTEMP = C( J ) 00210 STEMP = S( J ) 00211 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00212 DO 30 I = 1, N 00213 TEMP = A( J+1, I ) 00214 A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I ) 00215 A( J, I ) = STEMP*TEMP + CTEMP*A( J, I ) 00216 30 CONTINUE 00217 END IF 00218 40 CONTINUE 00219 END IF 00220 ELSE IF( LSAME( PIVOT, 'T' ) ) THEN 00221 IF( LSAME( DIRECT, 'F' ) ) THEN 00222 DO 60 J = 2, M 00223 CTEMP = C( J-1 ) 00224 STEMP = S( J-1 ) 00225 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00226 DO 50 I = 1, N 00227 TEMP = A( J, I ) 00228 A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I ) 00229 A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I ) 00230 50 CONTINUE 00231 END IF 00232 60 CONTINUE 00233 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 00234 DO 80 J = M, 2, -1 00235 CTEMP = C( J-1 ) 00236 STEMP = S( J-1 ) 00237 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00238 DO 70 I = 1, N 00239 TEMP = A( J, I ) 00240 A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I ) 00241 A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I ) 00242 70 CONTINUE 00243 END IF 00244 80 CONTINUE 00245 END IF 00246 ELSE IF( LSAME( PIVOT, 'B' ) ) THEN 00247 IF( LSAME( DIRECT, 'F' ) ) THEN 00248 DO 100 J = 1, M - 1 00249 CTEMP = C( J ) 00250 STEMP = S( J ) 00251 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00252 DO 90 I = 1, N 00253 TEMP = A( J, I ) 00254 A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP 00255 A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP 00256 90 CONTINUE 00257 END IF 00258 100 CONTINUE 00259 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 00260 DO 120 J = M - 1, 1, -1 00261 CTEMP = C( J ) 00262 STEMP = S( J ) 00263 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00264 DO 110 I = 1, N 00265 TEMP = A( J, I ) 00266 A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP 00267 A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP 00268 110 CONTINUE 00269 END IF 00270 120 CONTINUE 00271 END IF 00272 END IF 00273 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00274 * 00275 * Form A * P**T 00276 * 00277 IF( LSAME( PIVOT, 'V' ) ) THEN 00278 IF( LSAME( DIRECT, 'F' ) ) THEN 00279 DO 140 J = 1, N - 1 00280 CTEMP = C( J ) 00281 STEMP = S( J ) 00282 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00283 DO 130 I = 1, M 00284 TEMP = A( I, J+1 ) 00285 A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J ) 00286 A( I, J ) = STEMP*TEMP + CTEMP*A( I, J ) 00287 130 CONTINUE 00288 END IF 00289 140 CONTINUE 00290 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 00291 DO 160 J = N - 1, 1, -1 00292 CTEMP = C( J ) 00293 STEMP = S( J ) 00294 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00295 DO 150 I = 1, M 00296 TEMP = A( I, J+1 ) 00297 A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J ) 00298 A( I, J ) = STEMP*TEMP + CTEMP*A( I, J ) 00299 150 CONTINUE 00300 END IF 00301 160 CONTINUE 00302 END IF 00303 ELSE IF( LSAME( PIVOT, 'T' ) ) THEN 00304 IF( LSAME( DIRECT, 'F' ) ) THEN 00305 DO 180 J = 2, N 00306 CTEMP = C( J-1 ) 00307 STEMP = S( J-1 ) 00308 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00309 DO 170 I = 1, M 00310 TEMP = A( I, J ) 00311 A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 ) 00312 A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 ) 00313 170 CONTINUE 00314 END IF 00315 180 CONTINUE 00316 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 00317 DO 200 J = N, 2, -1 00318 CTEMP = C( J-1 ) 00319 STEMP = S( J-1 ) 00320 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00321 DO 190 I = 1, M 00322 TEMP = A( I, J ) 00323 A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 ) 00324 A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 ) 00325 190 CONTINUE 00326 END IF 00327 200 CONTINUE 00328 END IF 00329 ELSE IF( LSAME( PIVOT, 'B' ) ) THEN 00330 IF( LSAME( DIRECT, 'F' ) ) THEN 00331 DO 220 J = 1, N - 1 00332 CTEMP = C( J ) 00333 STEMP = S( J ) 00334 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00335 DO 210 I = 1, M 00336 TEMP = A( I, J ) 00337 A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP 00338 A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP 00339 210 CONTINUE 00340 END IF 00341 220 CONTINUE 00342 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 00343 DO 240 J = N - 1, 1, -1 00344 CTEMP = C( J ) 00345 STEMP = S( J ) 00346 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 00347 DO 230 I = 1, M 00348 TEMP = A( I, J ) 00349 A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP 00350 A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP 00351 230 CONTINUE 00352 END IF 00353 240 CONTINUE 00354 END IF 00355 END IF 00356 END IF 00357 * 00358 RETURN 00359 * 00360 * End of DLASR 00361 * 00362 END
1.7.3 
