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 #ifndef TEMPLATE_LAPACK_LARRF_HEADER
00036 #define TEMPLATE_LAPACK_LARRF_HEADER
00037
00038
00039 template<class Treal>
00040 int template_lapack_larrf(integer *n, Treal *d__, Treal *l,
00041 Treal *ld, integer *clstrt, integer *clend, Treal *w,
00042 Treal *wgap, Treal *werr, Treal *spdiam, Treal *
00043 clgapl, Treal *clgapr, Treal *pivmin, Treal *sigma,
00044 Treal *dplus, Treal *lplus, Treal *work, integer *info)
00045 {
00046
00047 integer i__1;
00048 Treal d__1, d__2, d__3;
00049
00050
00051 integer i__;
00052 Treal s, bestshift, smlgrowth, eps, tmp, max1, max2, rrr1, rrr2,
00053 znm2, growthbound, fail, fact, oldp;
00054 integer indx;
00055 Treal prod;
00056 integer ktry;
00057 Treal fail2, avgap, ldmax, rdmax;
00058 integer shift;
00059 logical dorrr1;
00060 Treal ldelta;
00061 logical nofail;
00062 Treal mingap, lsigma, rdelta;
00063 logical forcer;
00064 Treal rsigma, clwdth;
00065 logical sawnan1, sawnan2, tryrrr1;
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
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 integer c__1 = 1;
00170
00171
00172
00173
00174 --work;
00175 --lplus;
00176 --dplus;
00177 --werr;
00178 --wgap;
00179 --w;
00180 --ld;
00181 --l;
00182 --d__;
00183
00184
00185 indx = 0;
00186
00187 *info = 0;
00188 fact = 2.;
00189 eps = template_lapack_lamch("Precision", (Treal)0);
00190 shift = 0;
00191 forcer = FALSE_;
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203 nofail = TRUE_;
00204
00205
00206 clwdth = (d__1 = w[*clend] - w[*clstrt], absMACRO(d__1)) + werr[*clend] + werr[
00207 *clstrt];
00208 avgap = clwdth / (Treal) (*clend - *clstrt);
00209 mingap = minMACRO(*clgapl,*clgapr);
00210
00211
00212 d__1 = w[*clstrt], d__2 = w[*clend];
00213 lsigma = minMACRO(d__1,d__2) - werr[*clstrt];
00214
00215 d__1 = w[*clstrt], d__2 = w[*clend];
00216 rsigma = maxMACRO(d__1,d__2) + werr[*clend];
00217
00218 lsigma -= absMACRO(lsigma) * 4. * eps;
00219 rsigma += absMACRO(rsigma) * 4. * eps;
00220
00221 ldmax = mingap * .25 + *pivmin * 2.;
00222 rdmax = mingap * .25 + *pivmin * 2.;
00223
00224 d__1 = avgap, d__2 = wgap[*clstrt];
00225 ldelta = maxMACRO(d__1,d__2) / fact;
00226
00227 d__1 = avgap, d__2 = wgap[*clend - 1];
00228 rdelta = maxMACRO(d__1,d__2) / fact;
00229
00230
00231
00232 s = template_lapack_lamch("S", (Treal)0);
00233 smlgrowth = 1. / s;
00234 fail = (Treal) (*n - 1) * mingap / (*spdiam * eps);
00235 fail2 = (Treal) (*n - 1) * mingap / (*spdiam * template_blas_sqrt(eps));
00236 bestshift = lsigma;
00237
00238
00239 ktry = 0;
00240 growthbound = *spdiam * 8.;
00241 L5:
00242 sawnan1 = FALSE_;
00243 sawnan2 = FALSE_;
00244
00245 ldelta = minMACRO(ldmax,ldelta);
00246 rdelta = minMACRO(rdmax,rdelta);
00247
00248
00249
00250 s = -lsigma;
00251 dplus[1] = d__[1] + s;
00252 if (absMACRO(dplus[1]) < *pivmin) {
00253 dplus[1] = -(*pivmin);
00254
00255
00256 sawnan1 = TRUE_;
00257 }
00258 max1 = absMACRO(dplus[1]);
00259 i__1 = *n - 1;
00260 for (i__ = 1; i__ <= i__1; ++i__) {
00261 lplus[i__] = ld[i__] / dplus[i__];
00262 s = s * lplus[i__] * l[i__] - lsigma;
00263 dplus[i__ + 1] = d__[i__ + 1] + s;
00264 if ((d__1 = dplus[i__ + 1], absMACRO(d__1)) < *pivmin) {
00265 dplus[i__ + 1] = -(*pivmin);
00266
00267
00268 sawnan1 = TRUE_;
00269 }
00270
00271 d__2 = max1, d__3 = (d__1 = dplus[i__ + 1], absMACRO(d__1));
00272 max1 = maxMACRO(d__2,d__3);
00273
00274 }
00275 sawnan1 = sawnan1 || template_lapack_isnan(&max1);
00276 if (forcer || ( max1 <= growthbound && ! sawnan1 ) ) {
00277 *sigma = lsigma;
00278 shift = 1;
00279 goto L100;
00280 }
00281
00282 s = -rsigma;
00283 work[1] = d__[1] + s;
00284 if (absMACRO(work[1]) < *pivmin) {
00285 work[1] = -(*pivmin);
00286
00287
00288 sawnan2 = TRUE_;
00289 }
00290 max2 = absMACRO(work[1]);
00291 i__1 = *n - 1;
00292 for (i__ = 1; i__ <= i__1; ++i__) {
00293 work[*n + i__] = ld[i__] / work[i__];
00294 s = s * work[*n + i__] * l[i__] - rsigma;
00295 work[i__ + 1] = d__[i__ + 1] + s;
00296 if ((d__1 = work[i__ + 1], absMACRO(d__1)) < *pivmin) {
00297 work[i__ + 1] = -(*pivmin);
00298
00299
00300 sawnan2 = TRUE_;
00301 }
00302
00303 d__2 = max2, d__3 = (d__1 = work[i__ + 1], absMACRO(d__1));
00304 max2 = maxMACRO(d__2,d__3);
00305
00306 }
00307 sawnan2 = sawnan2 || template_lapack_isnan(&max2);
00308 if (forcer || ( max2 <= growthbound && ! sawnan2 ) ) {
00309 *sigma = rsigma;
00310 shift = 2;
00311 goto L100;
00312 }
00313
00314
00315 if (sawnan1 && sawnan2) {
00316
00317 goto L50;
00318 } else {
00319 if (! sawnan1) {
00320 indx = 1;
00321 if (max1 <= smlgrowth) {
00322 smlgrowth = max1;
00323 bestshift = lsigma;
00324 }
00325 }
00326 if (! sawnan2) {
00327 if (sawnan1 || max2 <= max1) {
00328 indx = 2;
00329 }
00330 if (max2 <= smlgrowth) {
00331 smlgrowth = max2;
00332 bestshift = rsigma;
00333 }
00334 }
00335 }
00336
00337
00338
00339
00340
00341 if (clwdth < mingap / 128. && minMACRO(max1,max2) < fail2 && ! sawnan1 && !
00342 sawnan2) {
00343 dorrr1 = TRUE_;
00344 } else {
00345 dorrr1 = FALSE_;
00346 }
00347 tryrrr1 = TRUE_;
00348 if (tryrrr1 && dorrr1) {
00349 if (indx == 1) {
00350 tmp = (d__1 = dplus[*n], absMACRO(d__1));
00351 znm2 = 1.;
00352 prod = 1.;
00353 oldp = 1.;
00354 for (i__ = *n - 1; i__ >= 1; --i__) {
00355 if (prod <= eps) {
00356 prod = dplus[i__ + 1] * work[*n + i__ + 1] / (dplus[i__] *
00357 work[*n + i__]) * oldp;
00358 } else {
00359 prod *= (d__1 = work[*n + i__], absMACRO(d__1));
00360 }
00361 oldp = prod;
00362
00363 d__1 = prod;
00364 znm2 += d__1 * d__1;
00365
00366 d__2 = tmp, d__3 = (d__1 = dplus[i__] * prod, absMACRO(d__1));
00367 tmp = maxMACRO(d__2,d__3);
00368
00369 }
00370 rrr1 = tmp / (*spdiam * template_blas_sqrt(znm2));
00371 if (rrr1 <= 8.) {
00372 *sigma = lsigma;
00373 shift = 1;
00374 goto L100;
00375 }
00376 } else if (indx == 2) {
00377 tmp = (d__1 = work[*n], absMACRO(d__1));
00378 znm2 = 1.;
00379 prod = 1.;
00380 oldp = 1.;
00381 for (i__ = *n - 1; i__ >= 1; --i__) {
00382 if (prod <= eps) {
00383 prod = work[i__ + 1] * lplus[i__ + 1] / (work[i__] *
00384 lplus[i__]) * oldp;
00385 } else {
00386 prod *= (d__1 = lplus[i__], absMACRO(d__1));
00387 }
00388 oldp = prod;
00389
00390 d__1 = prod;
00391 znm2 += d__1 * d__1;
00392
00393 d__2 = tmp, d__3 = (d__1 = work[i__] * prod, absMACRO(d__1));
00394 tmp = maxMACRO(d__2,d__3);
00395
00396 }
00397 rrr2 = tmp / (*spdiam * template_blas_sqrt(znm2));
00398 if (rrr2 <= 8.) {
00399 *sigma = rsigma;
00400 shift = 2;
00401 goto L100;
00402 }
00403 }
00404 }
00405 L50:
00406 if (ktry < 1) {
00407
00408
00409
00410 d__1 = lsigma - ldelta, d__2 = lsigma - ldmax;
00411 lsigma = maxMACRO(d__1,d__2);
00412
00413 d__1 = rsigma + rdelta, d__2 = rsigma + rdmax;
00414 rsigma = minMACRO(d__1,d__2);
00415 ldelta *= 2.;
00416 rdelta *= 2.;
00417 ++ktry;
00418 goto L5;
00419 } else {
00420
00421
00422 if (smlgrowth < fail || nofail) {
00423 lsigma = bestshift;
00424 rsigma = bestshift;
00425 forcer = TRUE_;
00426 goto L5;
00427 } else {
00428 *info = 1;
00429 return 0;
00430 }
00431 }
00432 L100:
00433 if (shift == 1) {
00434 } else if (shift == 2) {
00435
00436 template_blas_copy(n, &work[1], &c__1, &dplus[1], &c__1);
00437 i__1 = *n - 1;
00438 template_blas_copy(&i__1, &work[*n + 1], &c__1, &lplus[1], &c__1);
00439 }
00440 return 0;
00441
00442
00443
00444 }
00445
00446 #endif