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_LAR1V_HEADER
00038 #define TEMPLATE_LAPACK_LAR1V_HEADER
00039
00040 template<class Treal>
00041 int template_lapack_lar1v(integer *n, integer *b1, integer *bn, Treal
00042 *lambda, Treal *d__, Treal *l, Treal *ld, Treal *
00043 lld, Treal *pivmin, Treal *gaptol, Treal *z__, logical
00044 *wantnc, integer *negcnt, Treal *ztz, Treal *mingma,
00045 integer *r__, integer *isuppz, Treal *nrminv, Treal *resid,
00046 Treal *rqcorr, Treal *work)
00047 {
00048
00049 integer i__1;
00050 Treal d__1, d__2, d__3;
00051
00052
00053
00054 integer i__;
00055 Treal s;
00056 integer r1, r2;
00057 Treal eps, tmp;
00058 integer neg1, neg2, indp, inds;
00059 Treal dplus;
00060 integer indlpl, indumn;
00061 Treal dminus;
00062 logical sawnan1, sawnan2;
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
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
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 --work;
00200 --isuppz;
00201 --z__;
00202 --lld;
00203 --ld;
00204 --l;
00205 --d__;
00206
00207
00208 eps = template_lapack_lamch("Precision", (Treal)0);
00209 if (*r__ == 0) {
00210 r1 = *b1;
00211 r2 = *bn;
00212 } else {
00213 r1 = *r__;
00214 r2 = *r__;
00215 }
00216
00217 indlpl = 0;
00218
00219 indumn = *n;
00220 inds = (*n << 1) + 1;
00221 indp = *n * 3 + 1;
00222 if (*b1 == 1) {
00223 work[inds] = 0.;
00224 } else {
00225 work[inds + *b1 - 1] = lld[*b1 - 1];
00226 }
00227
00228
00229
00230
00231 sawnan1 = FALSE_;
00232 neg1 = 0;
00233 s = work[inds + *b1 - 1] - *lambda;
00234 i__1 = r1 - 1;
00235 for (i__ = *b1; i__ <= i__1; ++i__) {
00236 dplus = d__[i__] + s;
00237 work[indlpl + i__] = ld[i__] / dplus;
00238 if (dplus < 0.) {
00239 ++neg1;
00240 }
00241 work[inds + i__] = s * work[indlpl + i__] * l[i__];
00242 s = work[inds + i__] - *lambda;
00243
00244 }
00245 sawnan1 = template_lapack_isnan(&s);
00246 if (sawnan1) {
00247 goto L60;
00248 }
00249 i__1 = r2 - 1;
00250 for (i__ = r1; i__ <= i__1; ++i__) {
00251 dplus = d__[i__] + s;
00252 work[indlpl + i__] = ld[i__] / dplus;
00253 work[inds + i__] = s * work[indlpl + i__] * l[i__];
00254 s = work[inds + i__] - *lambda;
00255
00256 }
00257 sawnan1 = template_lapack_isnan(&s);
00258
00259 L60:
00260 if (sawnan1) {
00261
00262 neg1 = 0;
00263 s = work[inds + *b1 - 1] - *lambda;
00264 i__1 = r1 - 1;
00265 for (i__ = *b1; i__ <= i__1; ++i__) {
00266 dplus = d__[i__] + s;
00267 if (absMACRO(dplus) < *pivmin) {
00268 dplus = -(*pivmin);
00269 }
00270 work[indlpl + i__] = ld[i__] / dplus;
00271 if (dplus < 0.) {
00272 ++neg1;
00273 }
00274 work[inds + i__] = s * work[indlpl + i__] * l[i__];
00275 if (work[indlpl + i__] == 0.) {
00276 work[inds + i__] = lld[i__];
00277 }
00278 s = work[inds + i__] - *lambda;
00279
00280 }
00281 i__1 = r2 - 1;
00282 for (i__ = r1; i__ <= i__1; ++i__) {
00283 dplus = d__[i__] + s;
00284 if (absMACRO(dplus) < *pivmin) {
00285 dplus = -(*pivmin);
00286 }
00287 work[indlpl + i__] = ld[i__] / dplus;
00288 work[inds + i__] = s * work[indlpl + i__] * l[i__];
00289 if (work[indlpl + i__] == 0.) {
00290 work[inds + i__] = lld[i__];
00291 }
00292 s = work[inds + i__] - *lambda;
00293
00294 }
00295 }
00296
00297
00298
00299
00300 sawnan2 = FALSE_;
00301 neg2 = 0;
00302 work[indp + *bn - 1] = d__[*bn] - *lambda;
00303 i__1 = r1;
00304 for (i__ = *bn - 1; i__ >= i__1; --i__) {
00305 dminus = lld[i__] + work[indp + i__];
00306 tmp = d__[i__] / dminus;
00307 if (dminus < 0.) {
00308 ++neg2;
00309 }
00310 work[indumn + i__] = l[i__] * tmp;
00311 work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
00312
00313 }
00314 tmp = work[indp + r1 - 1];
00315 sawnan2 = template_lapack_isnan(&tmp);
00316 if (sawnan2) {
00317
00318 neg2 = 0;
00319 i__1 = r1;
00320 for (i__ = *bn - 1; i__ >= i__1; --i__) {
00321 dminus = lld[i__] + work[indp + i__];
00322 if (absMACRO(dminus) < *pivmin) {
00323 dminus = -(*pivmin);
00324 }
00325 tmp = d__[i__] / dminus;
00326 if (dminus < 0.) {
00327 ++neg2;
00328 }
00329 work[indumn + i__] = l[i__] * tmp;
00330 work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
00331 if (tmp == 0.) {
00332 work[indp + i__ - 1] = d__[i__] - *lambda;
00333 }
00334
00335 }
00336 }
00337
00338
00339
00340
00341 *mingma = work[inds + r1 - 1] + work[indp + r1 - 1];
00342 if (*mingma < 0.) {
00343 ++neg1;
00344 }
00345 if (*wantnc) {
00346 *negcnt = neg1 + neg2;
00347 } else {
00348 *negcnt = -1;
00349 }
00350 if (absMACRO(*mingma) == 0.) {
00351 *mingma = eps * work[inds + r1 - 1];
00352 }
00353 *r__ = r1;
00354 i__1 = r2 - 1;
00355 for (i__ = r1; i__ <= i__1; ++i__) {
00356 tmp = work[inds + i__] + work[indp + i__];
00357 if (tmp == 0.) {
00358 tmp = eps * work[inds + i__];
00359 }
00360 if (absMACRO(tmp) <= absMACRO(*mingma)) {
00361 *mingma = tmp;
00362 *r__ = i__ + 1;
00363 }
00364
00365 }
00366
00367
00368
00369 isuppz[1] = *b1;
00370 isuppz[2] = *bn;
00371 z__[*r__] = 1.;
00372 *ztz = 1.;
00373
00374
00375
00376 if (! sawnan1 && ! sawnan2) {
00377 i__1 = *b1;
00378 for (i__ = *r__ - 1; i__ >= i__1; --i__) {
00379 z__[i__] = -(work[indlpl + i__] * z__[i__ + 1]);
00380 if (((d__1 = z__[i__], absMACRO(d__1)) + (d__2 = z__[i__ + 1], absMACRO(
00381 d__2))) * (d__3 = ld[i__], absMACRO(d__3)) < *gaptol) {
00382 z__[i__] = 0.;
00383 isuppz[1] = i__ + 1;
00384 goto L220;
00385 }
00386 *ztz += z__[i__] * z__[i__];
00387
00388 }
00389 L220:
00390 ;
00391 } else {
00392
00393 i__1 = *b1;
00394 for (i__ = *r__ - 1; i__ >= i__1; --i__) {
00395 if (z__[i__ + 1] == 0.) {
00396 z__[i__] = -(ld[i__ + 1] / ld[i__]) * z__[i__ + 2];
00397 } else {
00398 z__[i__] = -(work[indlpl + i__] * z__[i__ + 1]);
00399 }
00400 if (((d__1 = z__[i__], absMACRO(d__1)) + (d__2 = z__[i__ + 1], absMACRO(
00401 d__2))) * (d__3 = ld[i__], absMACRO(d__3)) < *gaptol) {
00402 z__[i__] = 0.;
00403 isuppz[1] = i__ + 1;
00404 goto L240;
00405 }
00406 *ztz += z__[i__] * z__[i__];
00407
00408 }
00409 L240:
00410 ;
00411 }
00412
00413 if (! sawnan1 && ! sawnan2) {
00414 i__1 = *bn - 1;
00415 for (i__ = *r__; i__ <= i__1; ++i__) {
00416 z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
00417 if (((d__1 = z__[i__], absMACRO(d__1)) + (d__2 = z__[i__ + 1], absMACRO(
00418 d__2))) * (d__3 = ld[i__], absMACRO(d__3)) < *gaptol) {
00419 z__[i__ + 1] = 0.;
00420 isuppz[2] = i__;
00421 goto L260;
00422 }
00423 *ztz += z__[i__ + 1] * z__[i__ + 1];
00424
00425 }
00426 L260:
00427 ;
00428 } else {
00429
00430 i__1 = *bn - 1;
00431 for (i__ = *r__; i__ <= i__1; ++i__) {
00432 if (z__[i__] == 0.) {
00433 z__[i__ + 1] = -(ld[i__ - 1] / ld[i__]) * z__[i__ - 1];
00434 } else {
00435 z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
00436 }
00437 if (((d__1 = z__[i__], absMACRO(d__1)) + (d__2 = z__[i__ + 1], absMACRO(
00438 d__2))) * (d__3 = ld[i__], absMACRO(d__3)) < *gaptol) {
00439 z__[i__ + 1] = 0.;
00440 isuppz[2] = i__;
00441 goto L280;
00442 }
00443 *ztz += z__[i__ + 1] * z__[i__ + 1];
00444
00445 }
00446 L280:
00447 ;
00448 }
00449
00450
00451
00452 tmp = 1. / *ztz;
00453 *nrminv = template_blas_sqrt(tmp);
00454 *resid = absMACRO(*mingma) * *nrminv;
00455 *rqcorr = *mingma * tmp;
00456
00457
00458 return 0;
00459
00460
00461
00462 }
00463
00464 #endif