00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #ifndef TEMPLATE_LAPACK_LARRF_HEADER
00038 #define TEMPLATE_LAPACK_LARRF_HEADER
00039
00040
00041 template<class Treal>
00042 int template_lapack_larrf(integer *n, Treal *d__, Treal *l,
00043 Treal *ld, integer *clstrt, integer *clend, Treal *w,
00044 Treal *wgap, Treal *werr, Treal *spdiam, Treal *
00045 clgapl, Treal *clgapr, Treal *pivmin, Treal *sigma,
00046 Treal *dplus, Treal *lplus, Treal *work, integer *info)
00047 {
00048
00049 integer i__1;
00050 Treal d__1, d__2, d__3;
00051
00052
00053 integer i__;
00054 Treal s, bestshift, smlgrowth, eps, tmp, max1, max2, rrr1, rrr2,
00055 znm2, growthbound, fail, fact, oldp;
00056 integer indx;
00057 Treal prod;
00058 integer ktry;
00059 Treal fail2, avgap, ldmax, rdmax;
00060 integer shift;
00061 logical dorrr1;
00062 Treal ldelta;
00063 logical nofail;
00064 Treal mingap, lsigma, rdelta;
00065 logical forcer;
00066 Treal rsigma, clwdth;
00067 logical sawnan1, sawnan2, tryrrr1;
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
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 integer c__1 = 1;
00172
00173
00174
00175
00176 --work;
00177 --lplus;
00178 --dplus;
00179 --werr;
00180 --wgap;
00181 --w;
00182 --ld;
00183 --l;
00184 --d__;
00185
00186
00187 indx = 0;
00188
00189 *info = 0;
00190 fact = 2.;
00191 eps = template_lapack_lamch("Precision", (Treal)0);
00192 shift = 0;
00193 forcer = FALSE_;
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 nofail = TRUE_;
00206
00207
00208 clwdth = (d__1 = w[*clend] - w[*clstrt], absMACRO(d__1)) + werr[*clend] + werr[
00209 *clstrt];
00210 avgap = clwdth / (Treal) (*clend - *clstrt);
00211 mingap = minMACRO(*clgapl,*clgapr);
00212
00213
00214 d__1 = w[*clstrt], d__2 = w[*clend];
00215 lsigma = minMACRO(d__1,d__2) - werr[*clstrt];
00216
00217 d__1 = w[*clstrt], d__2 = w[*clend];
00218 rsigma = maxMACRO(d__1,d__2) + werr[*clend];
00219
00220 lsigma -= absMACRO(lsigma) * 4. * eps;
00221 rsigma += absMACRO(rsigma) * 4. * eps;
00222
00223 ldmax = mingap * .25 + *pivmin * 2.;
00224 rdmax = mingap * .25 + *pivmin * 2.;
00225
00226 d__1 = avgap, d__2 = wgap[*clstrt];
00227 ldelta = maxMACRO(d__1,d__2) / fact;
00228
00229 d__1 = avgap, d__2 = wgap[*clend - 1];
00230 rdelta = maxMACRO(d__1,d__2) / fact;
00231
00232
00233
00234 s = template_lapack_lamch("S", (Treal)0);
00235 smlgrowth = 1. / s;
00236 fail = (Treal) (*n - 1) * mingap / (*spdiam * eps);
00237 fail2 = (Treal) (*n - 1) * mingap / (*spdiam * template_blas_sqrt(eps));
00238 bestshift = lsigma;
00239
00240
00241 ktry = 0;
00242 growthbound = *spdiam * 8.;
00243 L5:
00244 sawnan1 = FALSE_;
00245 sawnan2 = FALSE_;
00246
00247 ldelta = minMACRO(ldmax,ldelta);
00248 rdelta = minMACRO(rdmax,rdelta);
00249
00250
00251
00252 s = -lsigma;
00253 dplus[1] = d__[1] + s;
00254 if (absMACRO(dplus[1]) < *pivmin) {
00255 dplus[1] = -(*pivmin);
00256
00257
00258 sawnan1 = TRUE_;
00259 }
00260 max1 = absMACRO(dplus[1]);
00261 i__1 = *n - 1;
00262 for (i__ = 1; i__ <= i__1; ++i__) {
00263 lplus[i__] = ld[i__] / dplus[i__];
00264 s = s * lplus[i__] * l[i__] - lsigma;
00265 dplus[i__ + 1] = d__[i__ + 1] + s;
00266 if ((d__1 = dplus[i__ + 1], absMACRO(d__1)) < *pivmin) {
00267 dplus[i__ + 1] = -(*pivmin);
00268
00269
00270 sawnan1 = TRUE_;
00271 }
00272
00273 d__2 = max1, d__3 = (d__1 = dplus[i__ + 1], absMACRO(d__1));
00274 max1 = maxMACRO(d__2,d__3);
00275
00276 }
00277 sawnan1 = sawnan1 || template_lapack_isnan(&max1);
00278 if (forcer || ( max1 <= growthbound && ! sawnan1 ) ) {
00279 *sigma = lsigma;
00280 shift = 1;
00281 goto L100;
00282 }
00283
00284 s = -rsigma;
00285 work[1] = d__[1] + s;
00286 if (absMACRO(work[1]) < *pivmin) {
00287 work[1] = -(*pivmin);
00288
00289
00290 sawnan2 = TRUE_;
00291 }
00292 max2 = absMACRO(work[1]);
00293 i__1 = *n - 1;
00294 for (i__ = 1; i__ <= i__1; ++i__) {
00295 work[*n + i__] = ld[i__] / work[i__];
00296 s = s * work[*n + i__] * l[i__] - rsigma;
00297 work[i__ + 1] = d__[i__ + 1] + s;
00298 if ((d__1 = work[i__ + 1], absMACRO(d__1)) < *pivmin) {
00299 work[i__ + 1] = -(*pivmin);
00300
00301
00302 sawnan2 = TRUE_;
00303 }
00304
00305 d__2 = max2, d__3 = (d__1 = work[i__ + 1], absMACRO(d__1));
00306 max2 = maxMACRO(d__2,d__3);
00307
00308 }
00309 sawnan2 = sawnan2 || template_lapack_isnan(&max2);
00310 if (forcer || ( max2 <= growthbound && ! sawnan2 ) ) {
00311 *sigma = rsigma;
00312 shift = 2;
00313 goto L100;
00314 }
00315
00316
00317 if (sawnan1 && sawnan2) {
00318
00319 goto L50;
00320 } else {
00321 if (! sawnan1) {
00322 indx = 1;
00323 if (max1 <= smlgrowth) {
00324 smlgrowth = max1;
00325 bestshift = lsigma;
00326 }
00327 }
00328 if (! sawnan2) {
00329 if (sawnan1 || max2 <= max1) {
00330 indx = 2;
00331 }
00332 if (max2 <= smlgrowth) {
00333 smlgrowth = max2;
00334 bestshift = rsigma;
00335 }
00336 }
00337 }
00338
00339
00340
00341
00342
00343 if (clwdth < mingap / 128. && minMACRO(max1,max2) < fail2 && ! sawnan1 && !
00344 sawnan2) {
00345 dorrr1 = TRUE_;
00346 } else {
00347 dorrr1 = FALSE_;
00348 }
00349 tryrrr1 = TRUE_;
00350 if (tryrrr1 && dorrr1) {
00351 if (indx == 1) {
00352 tmp = (d__1 = dplus[*n], absMACRO(d__1));
00353 znm2 = 1.;
00354 prod = 1.;
00355 oldp = 1.;
00356 for (i__ = *n - 1; i__ >= 1; --i__) {
00357 if (prod <= eps) {
00358 prod = dplus[i__ + 1] * work[*n + i__ + 1] / (dplus[i__] *
00359 work[*n + i__]) * oldp;
00360 } else {
00361 prod *= (d__1 = work[*n + i__], absMACRO(d__1));
00362 }
00363 oldp = prod;
00364
00365 d__1 = prod;
00366 znm2 += d__1 * d__1;
00367
00368 d__2 = tmp, d__3 = (d__1 = dplus[i__] * prod, absMACRO(d__1));
00369 tmp = maxMACRO(d__2,d__3);
00370
00371 }
00372 rrr1 = tmp / (*spdiam * template_blas_sqrt(znm2));
00373 if (rrr1 <= 8.) {
00374 *sigma = lsigma;
00375 shift = 1;
00376 goto L100;
00377 }
00378 } else if (indx == 2) {
00379 tmp = (d__1 = work[*n], absMACRO(d__1));
00380 znm2 = 1.;
00381 prod = 1.;
00382 oldp = 1.;
00383 for (i__ = *n - 1; i__ >= 1; --i__) {
00384 if (prod <= eps) {
00385 prod = work[i__ + 1] * lplus[i__ + 1] / (work[i__] *
00386 lplus[i__]) * oldp;
00387 } else {
00388 prod *= (d__1 = lplus[i__], absMACRO(d__1));
00389 }
00390 oldp = prod;
00391
00392 d__1 = prod;
00393 znm2 += d__1 * d__1;
00394
00395 d__2 = tmp, d__3 = (d__1 = work[i__] * prod, absMACRO(d__1));
00396 tmp = maxMACRO(d__2,d__3);
00397
00398 }
00399 rrr2 = tmp / (*spdiam * template_blas_sqrt(znm2));
00400 if (rrr2 <= 8.) {
00401 *sigma = rsigma;
00402 shift = 2;
00403 goto L100;
00404 }
00405 }
00406 }
00407 L50:
00408 if (ktry < 1) {
00409
00410
00411
00412 d__1 = lsigma - ldelta, d__2 = lsigma - ldmax;
00413 lsigma = maxMACRO(d__1,d__2);
00414
00415 d__1 = rsigma + rdelta, d__2 = rsigma + rdmax;
00416 rsigma = minMACRO(d__1,d__2);
00417 ldelta *= 2.;
00418 rdelta *= 2.;
00419 ++ktry;
00420 goto L5;
00421 } else {
00422
00423
00424 if (smlgrowth < fail || nofail) {
00425 lsigma = bestshift;
00426 rsigma = bestshift;
00427 forcer = TRUE_;
00428 goto L5;
00429 } else {
00430 *info = 1;
00431 return 0;
00432 }
00433 }
00434 L100:
00435 if (shift == 1) {
00436 } else if (shift == 2) {
00437
00438 template_blas_copy(n, &work[1], &c__1, &dplus[1], &c__1);
00439 i__1 = *n - 1;
00440 template_blas_copy(&i__1, &work[*n + 1], &c__1, &lplus[1], &c__1);
00441 }
00442 return 0;
00443
00444
00445
00446 }
00447
00448 #endif