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