00001 SUBROUTINE CLAROT( LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT, 00002 $ XRIGHT ) 00003 * 00004 * -- LAPACK auxiliary test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 LOGICAL LLEFT, LRIGHT, LROWS 00010 INTEGER LDA, NL 00011 COMPLEX C, S, XLEFT, XRIGHT 00012 * .. 00013 * .. Array Arguments .. 00014 COMPLEX A( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CLAROT applies a (Givens) rotation to two adjacent rows or 00021 * columns, where one element of the first and/or last column/row 00022 * for use on matrices stored in some format other than GE, so 00023 * that elements of the matrix may be used or modified for which 00024 * no array element is provided. 00025 * 00026 * One example is a symmetric matrix in SB format (bandwidth=4), for 00027 * which UPLO='L': Two adjacent rows will have the format: 00028 * 00029 * row j: * * * * * . . . . 00030 * row j+1: * * * * * . . . . 00031 * 00032 * '*' indicates elements for which storage is provided, 00033 * '.' indicates elements for which no storage is provided, but 00034 * are not necessarily zero; their values are determined by 00035 * symmetry. ' ' indicates elements which are necessarily zero, 00036 * and have no storage provided. 00037 * 00038 * Those columns which have two '*'s can be handled by SROT. 00039 * Those columns which have no '*'s can be ignored, since as long 00040 * as the Givens rotations are carefully applied to preserve 00041 * symmetry, their values are determined. 00042 * Those columns which have one '*' have to be handled separately, 00043 * by using separate variables "p" and "q": 00044 * 00045 * row j: * * * * * p . . . 00046 * row j+1: q * * * * * . . . . 00047 * 00048 * The element p would have to be set correctly, then that column 00049 * is rotated, setting p to its new value. The next call to 00050 * CLAROT would rotate columns j and j+1, using p, and restore 00051 * symmetry. The element q would start out being zero, and be 00052 * made non-zero by the rotation. Later, rotations would presumably 00053 * be chosen to zero q out. 00054 * 00055 * Typical Calling Sequences: rotating the i-th and (i+1)-st rows. 00056 * ------- ------- --------- 00057 * 00058 * General dense matrix: 00059 * 00060 * CALL CLAROT(.TRUE.,.FALSE.,.FALSE., N, C,S, 00061 * A(i,1),LDA, DUMMY, DUMMY) 00062 * 00063 * General banded matrix in GB format: 00064 * 00065 * j = MAX(1, i-KL ) 00066 * NL = MIN( N, i+KU+1 ) + 1-j 00067 * CALL CLAROT( .TRUE., i-KL.GE.1, i+KU.LT.N, NL, C,S, 00068 * A(KU+i+1-j,j),LDA-1, XLEFT, XRIGHT ) 00069 * 00070 * [ note that i+1-j is just MIN(i,KL+1) ] 00071 * 00072 * Symmetric banded matrix in SY format, bandwidth K, 00073 * lower triangle only: 00074 * 00075 * j = MAX(1, i-K ) 00076 * NL = MIN( K+1, i ) + 1 00077 * CALL CLAROT( .TRUE., i-K.GE.1, .TRUE., NL, C,S, 00078 * A(i,j), LDA, XLEFT, XRIGHT ) 00079 * 00080 * Same, but upper triangle only: 00081 * 00082 * NL = MIN( K+1, N-i ) + 1 00083 * CALL CLAROT( .TRUE., .TRUE., i+K.LT.N, NL, C,S, 00084 * A(i,i), LDA, XLEFT, XRIGHT ) 00085 * 00086 * Symmetric banded matrix in SB format, bandwidth K, 00087 * lower triangle only: 00088 * 00089 * [ same as for SY, except:] 00090 * . . . . 00091 * A(i+1-j,j), LDA-1, XLEFT, XRIGHT ) 00092 * 00093 * [ note that i+1-j is just MIN(i,K+1) ] 00094 * 00095 * Same, but upper triangle only: 00096 * . . . 00097 * A(K+1,i), LDA-1, XLEFT, XRIGHT ) 00098 * 00099 * Rotating columns is just the transpose of rotating rows, except 00100 * for GB and SB: (rotating columns i and i+1) 00101 * 00102 * GB: 00103 * j = MAX(1, i-KU ) 00104 * NL = MIN( N, i+KL+1 ) + 1-j 00105 * CALL CLAROT( .TRUE., i-KU.GE.1, i+KL.LT.N, NL, C,S, 00106 * A(KU+j+1-i,i),LDA-1, XTOP, XBOTTM ) 00107 * 00108 * [note that KU+j+1-i is just MAX(1,KU+2-i)] 00109 * 00110 * SB: (upper triangle) 00111 * 00112 * . . . . . . 00113 * A(K+j+1-i,i),LDA-1, XTOP, XBOTTM ) 00114 * 00115 * SB: (lower triangle) 00116 * 00117 * . . . . . . 00118 * A(1,i),LDA-1, XTOP, XBOTTM ) 00119 * 00120 * Arguments 00121 * ========= 00122 * 00123 * LROWS - LOGICAL 00124 * If .TRUE., then CLAROT will rotate two rows. If .FALSE., 00125 * then it will rotate two columns. 00126 * Not modified. 00127 * 00128 * LLEFT - LOGICAL 00129 * If .TRUE., then XLEFT will be used instead of the 00130 * corresponding element of A for the first element in the 00131 * second row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.) 00132 * If .FALSE., then the corresponding element of A will be 00133 * used. 00134 * Not modified. 00135 * 00136 * LRIGHT - LOGICAL 00137 * If .TRUE., then XRIGHT will be used instead of the 00138 * corresponding element of A for the last element in the 00139 * first row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.) If 00140 * .FALSE., then the corresponding element of A will be used. 00141 * Not modified. 00142 * 00143 * NL - INTEGER 00144 * The length of the rows (if LROWS=.TRUE.) or columns (if 00145 * LROWS=.FALSE.) to be rotated. If XLEFT and/or XRIGHT are 00146 * used, the columns/rows they are in should be included in 00147 * NL, e.g., if LLEFT = LRIGHT = .TRUE., then NL must be at 00148 * least 2. The number of rows/columns to be rotated 00149 * exclusive of those involving XLEFT and/or XRIGHT may 00150 * not be negative, i.e., NL minus how many of LLEFT and 00151 * LRIGHT are .TRUE. must be at least zero; if not, XERBLA 00152 * will be called. 00153 * Not modified. 00154 * 00155 * C, S - COMPLEX 00156 * Specify the Givens rotation to be applied. If LROWS is 00157 * true, then the matrix ( c s ) 00158 * ( _ _ ) 00159 * (-s c ) is applied from the left; 00160 * if false, then the transpose (not conjugated) thereof is 00161 * applied from the right. Note that in contrast to the 00162 * output of CROTG or to most versions of CROT, both C and S 00163 * are complex. For a Givens rotation, |C|**2 + |S|**2 should 00164 * be 1, but this is not checked. 00165 * Not modified. 00166 * 00167 * A - COMPLEX array. 00168 * The array containing the rows/columns to be rotated. The 00169 * first element of A should be the upper left element to 00170 * be rotated. 00171 * Read and modified. 00172 * 00173 * LDA - INTEGER 00174 * The "effective" leading dimension of A. If A contains 00175 * a matrix stored in GE, HE, or SY format, then this is just 00176 * the leading dimension of A as dimensioned in the calling 00177 * routine. If A contains a matrix stored in band (GB, HB, or 00178 * SB) format, then this should be *one less* than the leading 00179 * dimension used in the calling routine. Thus, if A were 00180 * dimensioned A(LDA,*) in CLAROT, then A(1,j) would be the 00181 * j-th element in the first of the two rows to be rotated, 00182 * and A(2,j) would be the j-th in the second, regardless of 00183 * how the array may be stored in the calling routine. [A 00184 * cannot, however, actually be dimensioned thus, since for 00185 * band format, the row number may exceed LDA, which is not 00186 * legal FORTRAN.] 00187 * If LROWS=.TRUE., then LDA must be at least 1, otherwise 00188 * it must be at least NL minus the number of .TRUE. values 00189 * in XLEFT and XRIGHT. 00190 * Not modified. 00191 * 00192 * XLEFT - COMPLEX 00193 * If LLEFT is .TRUE., then XLEFT will be used and modified 00194 * instead of A(2,1) (if LROWS=.TRUE.) or A(1,2) 00195 * (if LROWS=.FALSE.). 00196 * Read and modified. 00197 * 00198 * XRIGHT - COMPLEX 00199 * If LRIGHT is .TRUE., then XRIGHT will be used and modified 00200 * instead of A(1,NL) (if LROWS=.TRUE.) or A(NL,1) 00201 * (if LROWS=.FALSE.). 00202 * Read and modified. 00203 * 00204 * ===================================================================== 00205 * 00206 * .. Local Scalars .. 00207 INTEGER IINC, INEXT, IX, IY, IYT, J, NT 00208 COMPLEX TEMPX 00209 * .. 00210 * .. Local Arrays .. 00211 COMPLEX XT( 2 ), YT( 2 ) 00212 * .. 00213 * .. External Subroutines .. 00214 EXTERNAL XERBLA 00215 * .. 00216 * .. Intrinsic Functions .. 00217 INTRINSIC CONJG 00218 * .. 00219 * .. Executable Statements .. 00220 * 00221 * Set up indices, arrays for ends 00222 * 00223 IF( LROWS ) THEN 00224 IINC = LDA 00225 INEXT = 1 00226 ELSE 00227 IINC = 1 00228 INEXT = LDA 00229 END IF 00230 * 00231 IF( LLEFT ) THEN 00232 NT = 1 00233 IX = 1 + IINC 00234 IY = 2 + LDA 00235 XT( 1 ) = A( 1 ) 00236 YT( 1 ) = XLEFT 00237 ELSE 00238 NT = 0 00239 IX = 1 00240 IY = 1 + INEXT 00241 END IF 00242 * 00243 IF( LRIGHT ) THEN 00244 IYT = 1 + INEXT + ( NL-1 )*IINC 00245 NT = NT + 1 00246 XT( NT ) = XRIGHT 00247 YT( NT ) = A( IYT ) 00248 END IF 00249 * 00250 * Check for errors 00251 * 00252 IF( NL.LT.NT ) THEN 00253 CALL XERBLA( 'CLAROT', 4 ) 00254 RETURN 00255 END IF 00256 IF( LDA.LE.0 .OR. ( .NOT.LROWS .AND. LDA.LT.NL-NT ) ) THEN 00257 CALL XERBLA( 'CLAROT', 8 ) 00258 RETURN 00259 END IF 00260 * 00261 * Rotate 00262 * 00263 * CROT( NL-NT, A(IX),IINC, A(IY),IINC, C, S ) with complex C, S 00264 * 00265 DO 10 J = 0, NL - NT - 1 00266 TEMPX = C*A( IX+J*IINC ) + S*A( IY+J*IINC ) 00267 A( IY+J*IINC ) = -CONJG( S )*A( IX+J*IINC ) + 00268 $ CONJG( C )*A( IY+J*IINC ) 00269 A( IX+J*IINC ) = TEMPX 00270 10 CONTINUE 00271 * 00272 * CROT( NT, XT,1, YT,1, C, S ) with complex C, S 00273 * 00274 DO 20 J = 1, NT 00275 TEMPX = C*XT( J ) + S*YT( J ) 00276 YT( J ) = -CONJG( S )*XT( J ) + CONJG( C )*YT( J ) 00277 XT( J ) = TEMPX 00278 20 CONTINUE 00279 * 00280 * Stuff values back into XLEFT, XRIGHT, etc. 00281 * 00282 IF( LLEFT ) THEN 00283 A( 1 ) = XT( 1 ) 00284 XLEFT = YT( 1 ) 00285 END IF 00286 * 00287 IF( LRIGHT ) THEN 00288 XRIGHT = XT( NT ) 00289 A( IYT ) = YT( NT ) 00290 END IF 00291 * 00292 RETURN 00293 * 00294 * End of CLAROT 00295 * 00296 END
1.7.2