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