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 zlatsp_(char *uplo, integer *n, doublecomplex *x,
00022 integer *iseed)
00023 {
00024
00025 integer i__1, i__2, i__3;
00026 doublecomplex z__1, z__2, z__3;
00027
00028
00029 double sqrt(doublereal), z_abs(doublecomplex *);
00030
00031
00032 doublecomplex a, b, c__;
00033 integer j;
00034 doublecomplex r__;
00035 integer n5, jj;
00036 doublereal beta, alpha, alpha3;
00037 extern VOID zlarnd_(doublecomplex *, integer *,
00038 integer *);
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 --iseed;
00100 --x;
00101
00102
00103 alpha = (sqrt(17.) + 1.) / 8.;
00104 beta = alpha - .001;
00105 alpha3 = alpha * alpha * alpha;
00106
00107
00108
00109 i__1 = *n * (*n + 1) / 2;
00110 for (j = 1; j <= i__1; ++j) {
00111 i__2 = j;
00112 x[i__2].r = 0., x[i__2].i = 0.;
00113
00114 }
00115
00116
00117
00118 if (*(unsigned char *)uplo == 'U') {
00119 n5 = *n / 5;
00120 n5 = *n - n5 * 5 + 1;
00121
00122 jj = *n * (*n + 1) / 2;
00123 i__1 = n5;
00124 for (j = *n; j >= i__1; j += -5) {
00125 zlarnd_(&z__2, &c__5, &iseed[1]);
00126 z__1.r = alpha3 * z__2.r, z__1.i = alpha3 * z__2.i;
00127 a.r = z__1.r, a.i = z__1.i;
00128 zlarnd_(&z__2, &c__5, &iseed[1]);
00129 z__1.r = z__2.r / alpha, z__1.i = z__2.i / alpha;
00130 b.r = z__1.r, b.i = z__1.i;
00131 z__3.r = b.r * 2., z__3.i = b.i * 2.;
00132 z__2.r = z__3.r * 0. - z__3.i * 1., z__2.i = z__3.r * 1. + z__3.i
00133 * 0.;
00134 z__1.r = a.r - z__2.r, z__1.i = a.i - z__2.i;
00135 c__.r = z__1.r, c__.i = z__1.i;
00136 z__1.r = c__.r / beta, z__1.i = c__.i / beta;
00137 r__.r = z__1.r, r__.i = z__1.i;
00138 i__2 = jj;
00139 x[i__2].r = a.r, x[i__2].i = a.i;
00140 i__2 = jj - 2;
00141 x[i__2].r = b.r, x[i__2].i = b.i;
00142 jj -= j;
00143 i__2 = jj;
00144 zlarnd_(&z__1, &c__2, &iseed[1]);
00145 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00146 i__2 = jj - 1;
00147 x[i__2].r = r__.r, x[i__2].i = r__.i;
00148 jj -= j - 1;
00149 i__2 = jj;
00150 x[i__2].r = c__.r, x[i__2].i = c__.i;
00151 jj -= j - 2;
00152 i__2 = jj;
00153 zlarnd_(&z__1, &c__2, &iseed[1]);
00154 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00155 jj -= j - 3;
00156 i__2 = jj;
00157 zlarnd_(&z__1, &c__2, &iseed[1]);
00158 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00159 if (z_abs(&x[jj + (j - 3)]) > z_abs(&x[jj])) {
00160 i__2 = jj + (j - 4);
00161 i__3 = jj + (j - 3);
00162 z__1.r = x[i__3].r * 2., z__1.i = x[i__3].i * 2.;
00163 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00164 } else {
00165 i__2 = jj + (j - 4);
00166 i__3 = jj;
00167 z__1.r = x[i__3].r * 2., z__1.i = x[i__3].i * 2.;
00168 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00169 }
00170 jj -= j - 4;
00171
00172 }
00173
00174
00175
00176 j = n5 - 1;
00177 if (j > 2) {
00178 zlarnd_(&z__2, &c__5, &iseed[1]);
00179 z__1.r = alpha3 * z__2.r, z__1.i = alpha3 * z__2.i;
00180 a.r = z__1.r, a.i = z__1.i;
00181 zlarnd_(&z__2, &c__5, &iseed[1]);
00182 z__1.r = z__2.r / alpha, z__1.i = z__2.i / alpha;
00183 b.r = z__1.r, b.i = z__1.i;
00184 z__3.r = b.r * 2., z__3.i = b.i * 2.;
00185 z__2.r = z__3.r * 0. - z__3.i * 1., z__2.i = z__3.r * 1. + z__3.i
00186 * 0.;
00187 z__1.r = a.r - z__2.r, z__1.i = a.i - z__2.i;
00188 c__.r = z__1.r, c__.i = z__1.i;
00189 z__1.r = c__.r / beta, z__1.i = c__.i / beta;
00190 r__.r = z__1.r, r__.i = z__1.i;
00191 i__1 = jj;
00192 x[i__1].r = a.r, x[i__1].i = a.i;
00193 i__1 = jj - 2;
00194 x[i__1].r = b.r, x[i__1].i = b.i;
00195 jj -= j;
00196 i__1 = jj;
00197 zlarnd_(&z__1, &c__2, &iseed[1]);
00198 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00199 i__1 = jj - 1;
00200 x[i__1].r = r__.r, x[i__1].i = r__.i;
00201 jj -= j - 1;
00202 i__1 = jj;
00203 x[i__1].r = c__.r, x[i__1].i = c__.i;
00204 jj -= j - 2;
00205 j += -3;
00206 }
00207 if (j > 1) {
00208 i__1 = jj;
00209 zlarnd_(&z__1, &c__2, &iseed[1]);
00210 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00211 i__1 = jj - j;
00212 zlarnd_(&z__1, &c__2, &iseed[1]);
00213 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00214 if (z_abs(&x[jj]) > z_abs(&x[jj - j])) {
00215 i__1 = jj - 1;
00216 i__2 = jj;
00217 z__1.r = x[i__2].r * 2., z__1.i = x[i__2].i * 2.;
00218 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00219 } else {
00220 i__1 = jj - 1;
00221 i__2 = jj - j;
00222 z__1.r = x[i__2].r * 2., z__1.i = x[i__2].i * 2.;
00223 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00224 }
00225 jj = jj - j - (j - 1);
00226 j += -2;
00227 } else if (j == 1) {
00228 i__1 = jj;
00229 zlarnd_(&z__1, &c__2, &iseed[1]);
00230 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00231 --j;
00232 }
00233
00234
00235
00236 } else {
00237 n5 = *n / 5;
00238 n5 *= 5;
00239
00240 jj = 1;
00241 i__1 = n5;
00242 for (j = 1; j <= i__1; j += 5) {
00243 zlarnd_(&z__2, &c__5, &iseed[1]);
00244 z__1.r = alpha3 * z__2.r, z__1.i = alpha3 * z__2.i;
00245 a.r = z__1.r, a.i = z__1.i;
00246 zlarnd_(&z__2, &c__5, &iseed[1]);
00247 z__1.r = z__2.r / alpha, z__1.i = z__2.i / alpha;
00248 b.r = z__1.r, b.i = z__1.i;
00249 z__3.r = b.r * 2., z__3.i = b.i * 2.;
00250 z__2.r = z__3.r * 0. - z__3.i * 1., z__2.i = z__3.r * 1. + z__3.i
00251 * 0.;
00252 z__1.r = a.r - z__2.r, z__1.i = a.i - z__2.i;
00253 c__.r = z__1.r, c__.i = z__1.i;
00254 z__1.r = c__.r / beta, z__1.i = c__.i / beta;
00255 r__.r = z__1.r, r__.i = z__1.i;
00256 i__2 = jj;
00257 x[i__2].r = a.r, x[i__2].i = a.i;
00258 i__2 = jj + 2;
00259 x[i__2].r = b.r, x[i__2].i = b.i;
00260 jj += *n - j + 1;
00261 i__2 = jj;
00262 zlarnd_(&z__1, &c__2, &iseed[1]);
00263 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00264 i__2 = jj + 1;
00265 x[i__2].r = r__.r, x[i__2].i = r__.i;
00266 jj += *n - j;
00267 i__2 = jj;
00268 x[i__2].r = c__.r, x[i__2].i = c__.i;
00269 jj += *n - j - 1;
00270 i__2 = jj;
00271 zlarnd_(&z__1, &c__2, &iseed[1]);
00272 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00273 jj += *n - j - 2;
00274 i__2 = jj;
00275 zlarnd_(&z__1, &c__2, &iseed[1]);
00276 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00277 if (z_abs(&x[jj - (*n - j - 2)]) > z_abs(&x[jj])) {
00278 i__2 = jj - (*n - j - 2) + 1;
00279 i__3 = jj - (*n - j - 2);
00280 z__1.r = x[i__3].r * 2., z__1.i = x[i__3].i * 2.;
00281 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00282 } else {
00283 i__2 = jj - (*n - j - 2) + 1;
00284 i__3 = jj;
00285 z__1.r = x[i__3].r * 2., z__1.i = x[i__3].i * 2.;
00286 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00287 }
00288 jj += *n - j - 3;
00289
00290 }
00291
00292
00293
00294 j = n5 + 1;
00295 if (j < *n - 1) {
00296 zlarnd_(&z__2, &c__5, &iseed[1]);
00297 z__1.r = alpha3 * z__2.r, z__1.i = alpha3 * z__2.i;
00298 a.r = z__1.r, a.i = z__1.i;
00299 zlarnd_(&z__2, &c__5, &iseed[1]);
00300 z__1.r = z__2.r / alpha, z__1.i = z__2.i / alpha;
00301 b.r = z__1.r, b.i = z__1.i;
00302 z__3.r = b.r * 2., z__3.i = b.i * 2.;
00303 z__2.r = z__3.r * 0. - z__3.i * 1., z__2.i = z__3.r * 1. + z__3.i
00304 * 0.;
00305 z__1.r = a.r - z__2.r, z__1.i = a.i - z__2.i;
00306 c__.r = z__1.r, c__.i = z__1.i;
00307 z__1.r = c__.r / beta, z__1.i = c__.i / beta;
00308 r__.r = z__1.r, r__.i = z__1.i;
00309 i__1 = jj;
00310 x[i__1].r = a.r, x[i__1].i = a.i;
00311 i__1 = jj + 2;
00312 x[i__1].r = b.r, x[i__1].i = b.i;
00313 jj += *n - j + 1;
00314 i__1 = jj;
00315 zlarnd_(&z__1, &c__2, &iseed[1]);
00316 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00317 i__1 = jj + 1;
00318 x[i__1].r = r__.r, x[i__1].i = r__.i;
00319 jj += *n - j;
00320 i__1 = jj;
00321 x[i__1].r = c__.r, x[i__1].i = c__.i;
00322 jj += *n - j - 1;
00323 j += 3;
00324 }
00325 if (j < *n) {
00326 i__1 = jj;
00327 zlarnd_(&z__1, &c__2, &iseed[1]);
00328 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00329 i__1 = jj + (*n - j + 1);
00330 zlarnd_(&z__1, &c__2, &iseed[1]);
00331 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00332 if (z_abs(&x[jj]) > z_abs(&x[jj + (*n - j + 1)])) {
00333 i__1 = jj + 1;
00334 i__2 = jj;
00335 z__1.r = x[i__2].r * 2., z__1.i = x[i__2].i * 2.;
00336 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00337 } else {
00338 i__1 = jj + 1;
00339 i__2 = jj + (*n - j + 1);
00340 z__1.r = x[i__2].r * 2., z__1.i = x[i__2].i * 2.;
00341 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00342 }
00343 jj = jj + (*n - j + 1) + (*n - j);
00344 j += 2;
00345 } else if (j == *n) {
00346 i__1 = jj;
00347 zlarnd_(&z__1, &c__2, &iseed[1]);
00348 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00349 jj += *n - j + 1;
00350 ++j;
00351 }
00352 }
00353
00354 return 0;
00355
00356
00357
00358 }