00001 /* clarot.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__4 = 4; 00019 static integer c__8 = 8; 00020 00021 /* Subroutine */ int clarot_(logical *lrows, logical *lleft, logical *lright, 00022 integer *nl, complex *c__, complex *s, complex *a, integer *lda, 00023 complex *xleft, complex *xright) 00024 { 00025 /* System generated locals */ 00026 integer i__1, i__2, i__3, i__4; 00027 complex q__1, q__2, q__3, q__4, q__5, q__6; 00028 00029 /* Builtin functions */ 00030 void r_cnjg(complex *, complex *); 00031 00032 /* Local variables */ 00033 integer j, ix, iy, nt; 00034 complex xt[2], yt[2]; 00035 integer iyt, iinc, inext; 00036 complex tempx; 00037 extern /* Subroutine */ int xerbla_(char *, integer *); 00038 00039 00040 /* -- LAPACK auxiliary test routine (version 3.1) -- */ 00041 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00042 /* November 2006 */ 00043 00044 /* .. Scalar Arguments .. */ 00045 /* .. */ 00046 /* .. Array Arguments .. */ 00047 /* .. */ 00048 00049 /* Purpose */ 00050 /* ======= */ 00051 00052 /* CLAROT applies a (Givens) rotation to two adjacent rows or */ 00053 /* columns, where one element of the first and/or last column/row */ 00054 /* for use on matrices stored in some format other than GE, so */ 00055 /* that elements of the matrix may be used or modified for which */ 00056 /* no array element is provided. */ 00057 00058 /* One example is a symmetric matrix in SB format (bandwidth=4), for */ 00059 /* which UPLO='L': Two adjacent rows will have the format: */ 00060 00061 /* row j: * * * * * . . . . */ 00062 /* row j+1: * * * * * . . . . */ 00063 00064 /* '*' indicates elements for which storage is provided, */ 00065 /* '.' indicates elements for which no storage is provided, but */ 00066 /* are not necessarily zero; their values are determined by */ 00067 /* symmetry. ' ' indicates elements which are necessarily zero, */ 00068 /* and have no storage provided. */ 00069 00070 /* Those columns which have two '*'s can be handled by SROT. */ 00071 /* Those columns which have no '*'s can be ignored, since as long */ 00072 /* as the Givens rotations are carefully applied to preserve */ 00073 /* symmetry, their values are determined. */ 00074 /* Those columns which have one '*' have to be handled separately, */ 00075 /* by using separate variables "p" and "q": */ 00076 00077 /* row j: * * * * * p . . . */ 00078 /* row j+1: q * * * * * . . . . */ 00079 00080 /* The element p would have to be set correctly, then that column */ 00081 /* is rotated, setting p to its new value. The next call to */ 00082 /* CLAROT would rotate columns j and j+1, using p, and restore */ 00083 /* symmetry. The element q would start out being zero, and be */ 00084 /* made non-zero by the rotation. Later, rotations would presumably */ 00085 /* be chosen to zero q out. */ 00086 00087 /* Typical Calling Sequences: rotating the i-th and (i+1)-st rows. */ 00088 /* ------- ------- --------- */ 00089 00090 /* General dense matrix: */ 00091 00092 /* CALL CLAROT(.TRUE.,.FALSE.,.FALSE., N, C,S, */ 00093 /* A(i,1),LDA, DUMMY, DUMMY) */ 00094 00095 /* General banded matrix in GB format: */ 00096 00097 /* j = MAX(1, i-KL ) */ 00098 /* NL = MIN( N, i+KU+1 ) + 1-j */ 00099 /* CALL CLAROT( .TRUE., i-KL.GE.1, i+KU.LT.N, NL, C,S, */ 00100 /* A(KU+i+1-j,j),LDA-1, XLEFT, XRIGHT ) */ 00101 00102 /* [ note that i+1-j is just MIN(i,KL+1) ] */ 00103 00104 /* Symmetric banded matrix in SY format, bandwidth K, */ 00105 /* lower triangle only: */ 00106 00107 /* j = MAX(1, i-K ) */ 00108 /* NL = MIN( K+1, i ) + 1 */ 00109 /* CALL CLAROT( .TRUE., i-K.GE.1, .TRUE., NL, C,S, */ 00110 /* A(i,j), LDA, XLEFT, XRIGHT ) */ 00111 00112 /* Same, but upper triangle only: */ 00113 00114 /* NL = MIN( K+1, N-i ) + 1 */ 00115 /* CALL CLAROT( .TRUE., .TRUE., i+K.LT.N, NL, C,S, */ 00116 /* A(i,i), LDA, XLEFT, XRIGHT ) */ 00117 00118 /* Symmetric banded matrix in SB format, bandwidth K, */ 00119 /* lower triangle only: */ 00120 00121 /* [ same as for SY, except:] */ 00122 /* . . . . */ 00123 /* A(i+1-j,j), LDA-1, XLEFT, XRIGHT ) */ 00124 00125 /* [ note that i+1-j is just MIN(i,K+1) ] */ 00126 00127 /* Same, but upper triangle only: */ 00128 /* . . . */ 00129 /* A(K+1,i), LDA-1, XLEFT, XRIGHT ) */ 00130 00131 /* Rotating columns is just the transpose of rotating rows, except */ 00132 /* for GB and SB: (rotating columns i and i+1) */ 00133 00134 /* GB: */ 00135 /* j = MAX(1, i-KU ) */ 00136 /* NL = MIN( N, i+KL+1 ) + 1-j */ 00137 /* CALL CLAROT( .TRUE., i-KU.GE.1, i+KL.LT.N, NL, C,S, */ 00138 /* A(KU+j+1-i,i),LDA-1, XTOP, XBOTTM ) */ 00139 00140 /* [note that KU+j+1-i is just MAX(1,KU+2-i)] */ 00141 00142 /* SB: (upper triangle) */ 00143 00144 /* . . . . . . */ 00145 /* A(K+j+1-i,i),LDA-1, XTOP, XBOTTM ) */ 00146 00147 /* SB: (lower triangle) */ 00148 00149 /* . . . . . . */ 00150 /* A(1,i),LDA-1, XTOP, XBOTTM ) */ 00151 00152 /* Arguments */ 00153 /* ========= */ 00154 00155 /* LROWS - LOGICAL */ 00156 /* If .TRUE., then CLAROT will rotate two rows. If .FALSE., */ 00157 /* then it will rotate two columns. */ 00158 /* Not modified. */ 00159 00160 /* LLEFT - LOGICAL */ 00161 /* If .TRUE., then XLEFT will be used instead of the */ 00162 /* corresponding element of A for the first element in the */ 00163 /* second row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.) */ 00164 /* If .FALSE., then the corresponding element of A will be */ 00165 /* used. */ 00166 /* Not modified. */ 00167 00168 /* LRIGHT - LOGICAL */ 00169 /* If .TRUE., then XRIGHT will be used instead of the */ 00170 /* corresponding element of A for the last element in the */ 00171 /* first row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.) If */ 00172 /* .FALSE., then the corresponding element of A will be used. */ 00173 /* Not modified. */ 00174 00175 /* NL - INTEGER */ 00176 /* The length of the rows (if LROWS=.TRUE.) or columns (if */ 00177 /* LROWS=.FALSE.) to be rotated. If XLEFT and/or XRIGHT are */ 00178 /* used, the columns/rows they are in should be included in */ 00179 /* NL, e.g., if LLEFT = LRIGHT = .TRUE., then NL must be at */ 00180 /* least 2. The number of rows/columns to be rotated */ 00181 /* exclusive of those involving XLEFT and/or XRIGHT may */ 00182 /* not be negative, i.e., NL minus how many of LLEFT and */ 00183 /* LRIGHT are .TRUE. must be at least zero; if not, XERBLA */ 00184 /* will be called. */ 00185 /* Not modified. */ 00186 00187 /* C, S - COMPLEX */ 00188 /* Specify the Givens rotation to be applied. If LROWS is */ 00189 /* true, then the matrix ( c s ) */ 00190 /* ( _ _ ) */ 00191 /* (-s c ) is applied from the left; */ 00192 /* if false, then the transpose (not conjugated) thereof is */ 00193 /* applied from the right. Note that in contrast to the */ 00194 /* output of CROTG or to most versions of CROT, both C and S */ 00195 /* are complex. For a Givens rotation, |C|**2 + |S|**2 should */ 00196 /* be 1, but this is not checked. */ 00197 /* Not modified. */ 00198 00199 /* A - COMPLEX array. */ 00200 /* The array containing the rows/columns to be rotated. The */ 00201 /* first element of A should be the upper left element to */ 00202 /* be rotated. */ 00203 /* Read and modified. */ 00204 00205 /* LDA - INTEGER */ 00206 /* The "effective" leading dimension of A. If A contains */ 00207 /* a matrix stored in GE, HE, or SY format, then this is just */ 00208 /* the leading dimension of A as dimensioned in the calling */ 00209 /* routine. If A contains a matrix stored in band (GB, HB, or */ 00210 /* SB) format, then this should be *one less* than the leading */ 00211 /* dimension used in the calling routine. Thus, if A were */ 00212 /* dimensioned A(LDA,*) in CLAROT, then A(1,j) would be the */ 00213 /* j-th element in the first of the two rows to be rotated, */ 00214 /* and A(2,j) would be the j-th in the second, regardless of */ 00215 /* how the array may be stored in the calling routine. [A */ 00216 /* cannot, however, actually be dimensioned thus, since for */ 00217 /* band format, the row number may exceed LDA, which is not */ 00218 /* legal FORTRAN.] */ 00219 /* If LROWS=.TRUE., then LDA must be at least 1, otherwise */ 00220 /* it must be at least NL minus the number of .TRUE. values */ 00221 /* in XLEFT and XRIGHT. */ 00222 /* Not modified. */ 00223 00224 /* XLEFT - COMPLEX */ 00225 /* If LLEFT is .TRUE., then XLEFT will be used and modified */ 00226 /* instead of A(2,1) (if LROWS=.TRUE.) or A(1,2) */ 00227 /* (if LROWS=.FALSE.). */ 00228 /* Read and modified. */ 00229 00230 /* XRIGHT - COMPLEX */ 00231 /* If LRIGHT is .TRUE., then XRIGHT will be used and modified */ 00232 /* instead of A(1,NL) (if LROWS=.TRUE.) or A(NL,1) */ 00233 /* (if LROWS=.FALSE.). */ 00234 /* Read and modified. */ 00235 00236 /* ===================================================================== */ 00237 00238 /* .. Local Scalars .. */ 00239 /* .. */ 00240 /* .. Local Arrays .. */ 00241 /* .. */ 00242 /* .. External Subroutines .. */ 00243 /* .. */ 00244 /* .. Intrinsic Functions .. */ 00245 /* .. */ 00246 /* .. Executable Statements .. */ 00247 00248 /* Set up indices, arrays for ends */ 00249 00250 /* Parameter adjustments */ 00251 --a; 00252 00253 /* Function Body */ 00254 if (*lrows) { 00255 iinc = *lda; 00256 inext = 1; 00257 } else { 00258 iinc = 1; 00259 inext = *lda; 00260 } 00261 00262 if (*lleft) { 00263 nt = 1; 00264 ix = iinc + 1; 00265 iy = *lda + 2; 00266 xt[0].r = a[1].r, xt[0].i = a[1].i; 00267 yt[0].r = xleft->r, yt[0].i = xleft->i; 00268 } else { 00269 nt = 0; 00270 ix = 1; 00271 iy = inext + 1; 00272 } 00273 00274 if (*lright) { 00275 iyt = inext + 1 + (*nl - 1) * iinc; 00276 ++nt; 00277 i__1 = nt - 1; 00278 xt[i__1].r = xright->r, xt[i__1].i = xright->i; 00279 i__1 = nt - 1; 00280 i__2 = iyt; 00281 yt[i__1].r = a[i__2].r, yt[i__1].i = a[i__2].i; 00282 } 00283 00284 /* Check for errors */ 00285 00286 if (*nl < nt) { 00287 xerbla_("CLAROT", &c__4); 00288 return 0; 00289 } 00290 if (*lda <= 0 || ! (*lrows) && *lda < *nl - nt) { 00291 xerbla_("CLAROT", &c__8); 00292 return 0; 00293 } 00294 00295 /* Rotate */ 00296 00297 /* CROT( NL-NT, A(IX),IINC, A(IY),IINC, C, S ) with complex C, S */ 00298 00299 i__1 = *nl - nt - 1; 00300 for (j = 0; j <= i__1; ++j) { 00301 i__2 = ix + j * iinc; 00302 q__2.r = c__->r * a[i__2].r - c__->i * a[i__2].i, q__2.i = c__->r * a[ 00303 i__2].i + c__->i * a[i__2].r; 00304 i__3 = iy + j * iinc; 00305 q__3.r = s->r * a[i__3].r - s->i * a[i__3].i, q__3.i = s->r * a[i__3] 00306 .i + s->i * a[i__3].r; 00307 q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i; 00308 tempx.r = q__1.r, tempx.i = q__1.i; 00309 i__2 = iy + j * iinc; 00310 r_cnjg(&q__4, s); 00311 q__3.r = -q__4.r, q__3.i = -q__4.i; 00312 i__3 = ix + j * iinc; 00313 q__2.r = q__3.r * a[i__3].r - q__3.i * a[i__3].i, q__2.i = q__3.r * a[ 00314 i__3].i + q__3.i * a[i__3].r; 00315 r_cnjg(&q__6, c__); 00316 i__4 = iy + j * iinc; 00317 q__5.r = q__6.r * a[i__4].r - q__6.i * a[i__4].i, q__5.i = q__6.r * a[ 00318 i__4].i + q__6.i * a[i__4].r; 00319 q__1.r = q__2.r + q__5.r, q__1.i = q__2.i + q__5.i; 00320 a[i__2].r = q__1.r, a[i__2].i = q__1.i; 00321 i__2 = ix + j * iinc; 00322 a[i__2].r = tempx.r, a[i__2].i = tempx.i; 00323 /* L10: */ 00324 } 00325 00326 /* CROT( NT, XT,1, YT,1, C, S ) with complex C, S */ 00327 00328 i__1 = nt; 00329 for (j = 1; j <= i__1; ++j) { 00330 i__2 = j - 1; 00331 q__2.r = c__->r * xt[i__2].r - c__->i * xt[i__2].i, q__2.i = c__->r * 00332 xt[i__2].i + c__->i * xt[i__2].r; 00333 i__3 = j - 1; 00334 q__3.r = s->r * yt[i__3].r - s->i * yt[i__3].i, q__3.i = s->r * yt[ 00335 i__3].i + s->i * yt[i__3].r; 00336 q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i; 00337 tempx.r = q__1.r, tempx.i = q__1.i; 00338 i__2 = j - 1; 00339 r_cnjg(&q__4, s); 00340 q__3.r = -q__4.r, q__3.i = -q__4.i; 00341 i__3 = j - 1; 00342 q__2.r = q__3.r * xt[i__3].r - q__3.i * xt[i__3].i, q__2.i = q__3.r * 00343 xt[i__3].i + q__3.i * xt[i__3].r; 00344 r_cnjg(&q__6, c__); 00345 i__4 = j - 1; 00346 q__5.r = q__6.r * yt[i__4].r - q__6.i * yt[i__4].i, q__5.i = q__6.r * 00347 yt[i__4].i + q__6.i * yt[i__4].r; 00348 q__1.r = q__2.r + q__5.r, q__1.i = q__2.i + q__5.i; 00349 yt[i__2].r = q__1.r, yt[i__2].i = q__1.i; 00350 i__2 = j - 1; 00351 xt[i__2].r = tempx.r, xt[i__2].i = tempx.i; 00352 /* L20: */ 00353 } 00354 00355 /* Stuff values back into XLEFT, XRIGHT, etc. */ 00356 00357 if (*lleft) { 00358 a[1].r = xt[0].r, a[1].i = xt[0].i; 00359 xleft->r = yt[0].r, xleft->i = yt[0].i; 00360 } 00361 00362 if (*lright) { 00363 i__1 = nt - 1; 00364 xright->r = xt[i__1].r, xright->i = xt[i__1].i; 00365 i__1 = iyt; 00366 i__2 = nt - 1; 00367 a[i__1].r = yt[i__2].r, a[i__1].i = yt[i__2].i; 00368 } 00369 00370 return 0; 00371 00372 /* End of CLAROT */ 00373 00374 } /* clarot_ */