00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
00017
00018 static integer c__5 = 5;
00019 static integer c__2 = 2;
00020
00021 int clatsy_(char *uplo, integer *n, complex *x, integer *ldx,
00022 integer *iseed)
00023 {
00024
00025 integer x_dim1, x_offset, i__1, i__2, i__3;
00026 complex q__1, q__2, q__3;
00027
00028
00029 double sqrt(doublereal), c_abs(complex *);
00030
00031
00032 complex a, b, c__;
00033 integer i__, j;
00034 complex r__;
00035 integer n5;
00036 real beta, alpha, alpha3;
00037 extern VOID clarnd_(complex *, integer *, integer *);
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100 x_dim1 = *ldx;
00101 x_offset = 1 + x_dim1;
00102 x -= x_offset;
00103 --iseed;
00104
00105
00106 alpha = (sqrt(17.f) + 1.f) / 8.f;
00107 beta = alpha - .001f;
00108 alpha3 = alpha * alpha * alpha;
00109
00110
00111
00112 if (*(unsigned char *)uplo == 'U') {
00113
00114
00115
00116 i__1 = *n;
00117 for (j = 1; j <= i__1; ++j) {
00118 i__2 = j;
00119 for (i__ = 1; i__ <= i__2; ++i__) {
00120 i__3 = i__ + j * x_dim1;
00121 x[i__3].r = 0.f, x[i__3].i = 0.f;
00122
00123 }
00124
00125 }
00126 n5 = *n / 5;
00127 n5 = *n - n5 * 5 + 1;
00128
00129 i__1 = n5;
00130 for (i__ = *n; i__ >= i__1; i__ += -5) {
00131 clarnd_(&q__2, &c__5, &iseed[1]);
00132 q__1.r = alpha3 * q__2.r, q__1.i = alpha3 * q__2.i;
00133 a.r = q__1.r, a.i = q__1.i;
00134 clarnd_(&q__2, &c__5, &iseed[1]);
00135 q__1.r = q__2.r / alpha, q__1.i = q__2.i / alpha;
00136 b.r = q__1.r, b.i = q__1.i;
00137 q__3.r = b.r * 2.f, q__3.i = b.i * 2.f;
00138 q__2.r = q__3.r * 0.f - q__3.i * 1.f, q__2.i = q__3.r * 1.f +
00139 q__3.i * 0.f;
00140 q__1.r = a.r - q__2.r, q__1.i = a.i - q__2.i;
00141 c__.r = q__1.r, c__.i = q__1.i;
00142 q__1.r = c__.r / beta, q__1.i = c__.i / beta;
00143 r__.r = q__1.r, r__.i = q__1.i;
00144 i__2 = i__ + i__ * x_dim1;
00145 x[i__2].r = a.r, x[i__2].i = a.i;
00146 i__2 = i__ - 2 + i__ * x_dim1;
00147 x[i__2].r = b.r, x[i__2].i = b.i;
00148 i__2 = i__ - 2 + (i__ - 1) * x_dim1;
00149 x[i__2].r = r__.r, x[i__2].i = r__.i;
00150 i__2 = i__ - 2 + (i__ - 2) * x_dim1;
00151 x[i__2].r = c__.r, x[i__2].i = c__.i;
00152 i__2 = i__ - 1 + (i__ - 1) * x_dim1;
00153 clarnd_(&q__1, &c__2, &iseed[1]);
00154 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00155 i__2 = i__ - 3 + (i__ - 3) * x_dim1;
00156 clarnd_(&q__1, &c__2, &iseed[1]);
00157 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00158 i__2 = i__ - 4 + (i__ - 4) * x_dim1;
00159 clarnd_(&q__1, &c__2, &iseed[1]);
00160 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00161 if (c_abs(&x[i__ - 3 + (i__ - 3) * x_dim1]) > c_abs(&x[i__ - 4 + (
00162 i__ - 4) * x_dim1])) {
00163 i__2 = i__ - 4 + (i__ - 3) * x_dim1;
00164 i__3 = i__ - 3 + (i__ - 3) * x_dim1;
00165 q__1.r = x[i__3].r * 2.f, q__1.i = x[i__3].i * 2.f;
00166 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00167 } else {
00168 i__2 = i__ - 4 + (i__ - 3) * x_dim1;
00169 i__3 = i__ - 4 + (i__ - 4) * x_dim1;
00170 q__1.r = x[i__3].r * 2.f, q__1.i = x[i__3].i * 2.f;
00171 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00172 }
00173
00174 }
00175
00176
00177
00178 i__ = n5 - 1;
00179 if (i__ > 2) {
00180 clarnd_(&q__2, &c__5, &iseed[1]);
00181 q__1.r = alpha3 * q__2.r, q__1.i = alpha3 * q__2.i;
00182 a.r = q__1.r, a.i = q__1.i;
00183 clarnd_(&q__2, &c__5, &iseed[1]);
00184 q__1.r = q__2.r / alpha, q__1.i = q__2.i / alpha;
00185 b.r = q__1.r, b.i = q__1.i;
00186 q__3.r = b.r * 2.f, q__3.i = b.i * 2.f;
00187 q__2.r = q__3.r * 0.f - q__3.i * 1.f, q__2.i = q__3.r * 1.f +
00188 q__3.i * 0.f;
00189 q__1.r = a.r - q__2.r, q__1.i = a.i - q__2.i;
00190 c__.r = q__1.r, c__.i = q__1.i;
00191 q__1.r = c__.r / beta, q__1.i = c__.i / beta;
00192 r__.r = q__1.r, r__.i = q__1.i;
00193 i__1 = i__ + i__ * x_dim1;
00194 x[i__1].r = a.r, x[i__1].i = a.i;
00195 i__1 = i__ - 2 + i__ * x_dim1;
00196 x[i__1].r = b.r, x[i__1].i = b.i;
00197 i__1 = i__ - 2 + (i__ - 1) * x_dim1;
00198 x[i__1].r = r__.r, x[i__1].i = r__.i;
00199 i__1 = i__ - 2 + (i__ - 2) * x_dim1;
00200 x[i__1].r = c__.r, x[i__1].i = c__.i;
00201 i__1 = i__ - 1 + (i__ - 1) * x_dim1;
00202 clarnd_(&q__1, &c__2, &iseed[1]);
00203 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00204 i__ += -3;
00205 }
00206 if (i__ > 1) {
00207 i__1 = i__ + i__ * x_dim1;
00208 clarnd_(&q__1, &c__2, &iseed[1]);
00209 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00210 i__1 = i__ - 1 + (i__ - 1) * x_dim1;
00211 clarnd_(&q__1, &c__2, &iseed[1]);
00212 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00213 if (c_abs(&x[i__ + i__ * x_dim1]) > c_abs(&x[i__ - 1 + (i__ - 1) *
00214 x_dim1])) {
00215 i__1 = i__ - 1 + i__ * x_dim1;
00216 i__2 = i__ + i__ * x_dim1;
00217 q__1.r = x[i__2].r * 2.f, q__1.i = x[i__2].i * 2.f;
00218 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00219 } else {
00220 i__1 = i__ - 1 + i__ * x_dim1;
00221 i__2 = i__ - 1 + (i__ - 1) * x_dim1;
00222 q__1.r = x[i__2].r * 2.f, q__1.i = x[i__2].i * 2.f;
00223 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00224 }
00225 i__ += -2;
00226 } else if (i__ == 1) {
00227 i__1 = i__ + i__ * x_dim1;
00228 clarnd_(&q__1, &c__2, &iseed[1]);
00229 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00230 --i__;
00231 }
00232
00233
00234
00235 } else {
00236
00237
00238
00239 i__1 = *n;
00240 for (j = 1; j <= i__1; ++j) {
00241 i__2 = *n;
00242 for (i__ = j; i__ <= i__2; ++i__) {
00243 i__3 = i__ + j * x_dim1;
00244 x[i__3].r = 0.f, x[i__3].i = 0.f;
00245
00246 }
00247
00248 }
00249 n5 = *n / 5;
00250 n5 *= 5;
00251
00252 i__1 = n5;
00253 for (i__ = 1; i__ <= i__1; i__ += 5) {
00254 clarnd_(&q__2, &c__5, &iseed[1]);
00255 q__1.r = alpha3 * q__2.r, q__1.i = alpha3 * q__2.i;
00256 a.r = q__1.r, a.i = q__1.i;
00257 clarnd_(&q__2, &c__5, &iseed[1]);
00258 q__1.r = q__2.r / alpha, q__1.i = q__2.i / alpha;
00259 b.r = q__1.r, b.i = q__1.i;
00260 q__3.r = b.r * 2.f, q__3.i = b.i * 2.f;
00261 q__2.r = q__3.r * 0.f - q__3.i * 1.f, q__2.i = q__3.r * 1.f +
00262 q__3.i * 0.f;
00263 q__1.r = a.r - q__2.r, q__1.i = a.i - q__2.i;
00264 c__.r = q__1.r, c__.i = q__1.i;
00265 q__1.r = c__.r / beta, q__1.i = c__.i / beta;
00266 r__.r = q__1.r, r__.i = q__1.i;
00267 i__2 = i__ + i__ * x_dim1;
00268 x[i__2].r = a.r, x[i__2].i = a.i;
00269 i__2 = i__ + 2 + i__ * x_dim1;
00270 x[i__2].r = b.r, x[i__2].i = b.i;
00271 i__2 = i__ + 2 + (i__ + 1) * x_dim1;
00272 x[i__2].r = r__.r, x[i__2].i = r__.i;
00273 i__2 = i__ + 2 + (i__ + 2) * x_dim1;
00274 x[i__2].r = c__.r, x[i__2].i = c__.i;
00275 i__2 = i__ + 1 + (i__ + 1) * x_dim1;
00276 clarnd_(&q__1, &c__2, &iseed[1]);
00277 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00278 i__2 = i__ + 3 + (i__ + 3) * x_dim1;
00279 clarnd_(&q__1, &c__2, &iseed[1]);
00280 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00281 i__2 = i__ + 4 + (i__ + 4) * x_dim1;
00282 clarnd_(&q__1, &c__2, &iseed[1]);
00283 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00284 if (c_abs(&x[i__ + 3 + (i__ + 3) * x_dim1]) > c_abs(&x[i__ + 4 + (
00285 i__ + 4) * x_dim1])) {
00286 i__2 = i__ + 4 + (i__ + 3) * x_dim1;
00287 i__3 = i__ + 3 + (i__ + 3) * x_dim1;
00288 q__1.r = x[i__3].r * 2.f, q__1.i = x[i__3].i * 2.f;
00289 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00290 } else {
00291 i__2 = i__ + 4 + (i__ + 3) * x_dim1;
00292 i__3 = i__ + 4 + (i__ + 4) * x_dim1;
00293 q__1.r = x[i__3].r * 2.f, q__1.i = x[i__3].i * 2.f;
00294 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00295 }
00296
00297 }
00298
00299
00300
00301 i__ = n5 + 1;
00302 if (i__ < *n - 1) {
00303 clarnd_(&q__2, &c__5, &iseed[1]);
00304 q__1.r = alpha3 * q__2.r, q__1.i = alpha3 * q__2.i;
00305 a.r = q__1.r, a.i = q__1.i;
00306 clarnd_(&q__2, &c__5, &iseed[1]);
00307 q__1.r = q__2.r / alpha, q__1.i = q__2.i / alpha;
00308 b.r = q__1.r, b.i = q__1.i;
00309 q__3.r = b.r * 2.f, q__3.i = b.i * 2.f;
00310 q__2.r = q__3.r * 0.f - q__3.i * 1.f, q__2.i = q__3.r * 1.f +
00311 q__3.i * 0.f;
00312 q__1.r = a.r - q__2.r, q__1.i = a.i - q__2.i;
00313 c__.r = q__1.r, c__.i = q__1.i;
00314 q__1.r = c__.r / beta, q__1.i = c__.i / beta;
00315 r__.r = q__1.r, r__.i = q__1.i;
00316 i__1 = i__ + i__ * x_dim1;
00317 x[i__1].r = a.r, x[i__1].i = a.i;
00318 i__1 = i__ + 2 + i__ * x_dim1;
00319 x[i__1].r = b.r, x[i__1].i = b.i;
00320 i__1 = i__ + 2 + (i__ + 1) * x_dim1;
00321 x[i__1].r = r__.r, x[i__1].i = r__.i;
00322 i__1 = i__ + 2 + (i__ + 2) * x_dim1;
00323 x[i__1].r = c__.r, x[i__1].i = c__.i;
00324 i__1 = i__ + 1 + (i__ + 1) * x_dim1;
00325 clarnd_(&q__1, &c__2, &iseed[1]);
00326 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00327 i__ += 3;
00328 }
00329 if (i__ < *n) {
00330 i__1 = i__ + i__ * x_dim1;
00331 clarnd_(&q__1, &c__2, &iseed[1]);
00332 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00333 i__1 = i__ + 1 + (i__ + 1) * x_dim1;
00334 clarnd_(&q__1, &c__2, &iseed[1]);
00335 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00336 if (c_abs(&x[i__ + i__ * x_dim1]) > c_abs(&x[i__ + 1 + (i__ + 1) *
00337 x_dim1])) {
00338 i__1 = i__ + 1 + i__ * x_dim1;
00339 i__2 = i__ + i__ * x_dim1;
00340 q__1.r = x[i__2].r * 2.f, q__1.i = x[i__2].i * 2.f;
00341 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00342 } else {
00343 i__1 = i__ + 1 + i__ * x_dim1;
00344 i__2 = i__ + 1 + (i__ + 1) * x_dim1;
00345 q__1.r = x[i__2].r * 2.f, q__1.i = x[i__2].i * 2.f;
00346 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00347 }
00348 i__ += 2;
00349 } else if (i__ == *n) {
00350 i__1 = i__ + i__ * x_dim1;
00351 clarnd_(&q__1, &c__2, &iseed[1]);
00352 x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00353 ++i__;
00354 }
00355 }
00356
00357 return 0;
00358
00359
00360
00361 }