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_LASQ4_HEADER
00038 #define TEMPLATE_LAPACK_LASQ4_HEADER
00039
00040 template<class Treal>
00041 int template_lapack_lasq4(integer *i0, integer *n0, Treal *z__,
00042 integer *pp, integer *n0in, Treal *dmin__, Treal *dmin1,
00043 Treal *dmin2, Treal *dn, Treal *dn1, Treal *dn2,
00044 Treal *tau, integer *ttype, Treal *g)
00045 {
00046
00047 integer i__1;
00048 Treal d__1, d__2;
00049
00050
00051
00052 Treal s = 0;
00053 Treal a2, b1, b2;
00054 integer i4, nn, np;
00055 Treal gam, gap1, gap2;
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
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 --z__;
00141
00142
00143 if (*dmin__ <= 0.) {
00144 *tau = -(*dmin__);
00145 *ttype = -1;
00146 return 0;
00147 }
00148
00149 nn = (*n0 << 2) + *pp;
00150 if (*n0in == *n0) {
00151
00152
00153
00154 if (*dmin__ == *dn || *dmin__ == *dn1) {
00155
00156 b1 = template_blas_sqrt(z__[nn - 3]) * template_blas_sqrt(z__[nn - 5]);
00157 b2 = template_blas_sqrt(z__[nn - 7]) * template_blas_sqrt(z__[nn - 9]);
00158 a2 = z__[nn - 7] + z__[nn - 5];
00159
00160
00161
00162 if (*dmin__ == *dn && *dmin1 == *dn1) {
00163 gap2 = *dmin2 - a2 - *dmin2 * .25;
00164 if (gap2 > 0. && gap2 > b2) {
00165 gap1 = a2 - *dn - b2 / gap2 * b2;
00166 } else {
00167 gap1 = a2 - *dn - (b1 + b2);
00168 }
00169 if (gap1 > 0. && gap1 > b1) {
00170
00171 d__1 = *dn - b1 / gap1 * b1, d__2 = *dmin__ * .5;
00172 s = maxMACRO(d__1,d__2);
00173 *ttype = -2;
00174 } else {
00175 s = 0.;
00176 if (*dn > b1) {
00177 s = *dn - b1;
00178 }
00179 if (a2 > b1 + b2) {
00180
00181 d__1 = s, d__2 = a2 - (b1 + b2);
00182 s = minMACRO(d__1,d__2);
00183 }
00184
00185 d__1 = s, d__2 = *dmin__ * .333;
00186 s = maxMACRO(d__1,d__2);
00187 *ttype = -3;
00188 }
00189 } else {
00190
00191
00192
00193 *ttype = -4;
00194 s = *dmin__ * .25;
00195 if (*dmin__ == *dn) {
00196 gam = *dn;
00197 a2 = 0.;
00198 if (z__[nn - 5] > z__[nn - 7]) {
00199 return 0;
00200 }
00201 b2 = z__[nn - 5] / z__[nn - 7];
00202 np = nn - 9;
00203 } else {
00204 np = nn - (*pp << 1);
00205 b2 = z__[np - 2];
00206 gam = *dn1;
00207 if (z__[np - 4] > z__[np - 2]) {
00208 return 0;
00209 }
00210 a2 = z__[np - 4] / z__[np - 2];
00211 if (z__[nn - 9] > z__[nn - 11]) {
00212 return 0;
00213 }
00214 b2 = z__[nn - 9] / z__[nn - 11];
00215 np = nn - 13;
00216 }
00217
00218
00219
00220 a2 += b2;
00221 i__1 = (*i0 << 2) - 1 + *pp;
00222 for (i4 = np; i4 >= i__1; i4 += -4) {
00223 if (b2 == 0.) {
00224 goto L20;
00225 }
00226 b1 = b2;
00227 if (z__[i4] > z__[i4 - 2]) {
00228 return 0;
00229 }
00230 b2 *= z__[i4] / z__[i4 - 2];
00231 a2 += b2;
00232 if (maxMACRO(b2,b1) * 100. < a2 || .563 < a2) {
00233 goto L20;
00234 }
00235
00236 }
00237 L20:
00238 a2 *= 1.05;
00239
00240
00241
00242 if (a2 < .563) {
00243 s = gam * (1. - template_blas_sqrt(a2)) / (a2 + 1.);
00244 }
00245 }
00246 } else if (*dmin__ == *dn2) {
00247
00248
00249
00250 *ttype = -5;
00251 s = *dmin__ * .25;
00252
00253
00254
00255 np = nn - (*pp << 1);
00256 b1 = z__[np - 2];
00257 b2 = z__[np - 6];
00258 gam = *dn2;
00259 if (z__[np - 8] > b2 || z__[np - 4] > b1) {
00260 return 0;
00261 }
00262 a2 = z__[np - 8] / b2 * (z__[np - 4] / b1 + 1.);
00263
00264
00265
00266 if (*n0 - *i0 > 2) {
00267 b2 = z__[nn - 13] / z__[nn - 15];
00268 a2 += b2;
00269 i__1 = (*i0 << 2) - 1 + *pp;
00270 for (i4 = nn - 17; i4 >= i__1; i4 += -4) {
00271 if (b2 == 0.) {
00272 goto L40;
00273 }
00274 b1 = b2;
00275 if (z__[i4] > z__[i4 - 2]) {
00276 return 0;
00277 }
00278 b2 *= z__[i4] / z__[i4 - 2];
00279 a2 += b2;
00280 if (maxMACRO(b2,b1) * 100. < a2 || .563 < a2) {
00281 goto L40;
00282 }
00283
00284 }
00285 L40:
00286 a2 *= 1.05;
00287 }
00288
00289 if (a2 < .563) {
00290 s = gam * (1. - template_blas_sqrt(a2)) / (a2 + 1.);
00291 }
00292 } else {
00293
00294
00295
00296 if (*ttype == -6) {
00297 *g += (1. - *g) * .333;
00298 } else if (*ttype == -18) {
00299 *g = .083250000000000005;
00300 } else {
00301 *g = .25;
00302 }
00303 s = *g * *dmin__;
00304 *ttype = -6;
00305 }
00306
00307 } else if (*n0in == *n0 + 1) {
00308
00309
00310
00311 if (*dmin1 == *dn1 && *dmin2 == *dn2) {
00312
00313
00314
00315 *ttype = -7;
00316 s = *dmin1 * .333;
00317 if (z__[nn - 5] > z__[nn - 7]) {
00318 return 0;
00319 }
00320 b1 = z__[nn - 5] / z__[nn - 7];
00321 b2 = b1;
00322 if (b2 == 0.) {
00323 goto L60;
00324 }
00325 i__1 = (*i0 << 2) - 1 + *pp;
00326 for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
00327 a2 = b1;
00328 if (z__[i4] > z__[i4 - 2]) {
00329 return 0;
00330 }
00331 b1 *= z__[i4] / z__[i4 - 2];
00332 b2 += b1;
00333 if (maxMACRO(b1,a2) * 100. < b2) {
00334 goto L60;
00335 }
00336
00337 }
00338 L60:
00339 b2 = template_blas_sqrt(b2 * 1.05);
00340
00341 d__1 = b2;
00342 a2 = *dmin1 / (d__1 * d__1 + 1.);
00343 gap2 = *dmin2 * .5 - a2;
00344 if (gap2 > 0. && gap2 > b2 * a2) {
00345
00346 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
00347 s = maxMACRO(d__1,d__2);
00348 } else {
00349
00350 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
00351 s = maxMACRO(d__1,d__2);
00352 *ttype = -8;
00353 }
00354 } else {
00355
00356
00357
00358 s = *dmin1 * .25;
00359 if (*dmin1 == *dn1) {
00360 s = *dmin1 * .5;
00361 }
00362 *ttype = -9;
00363 }
00364
00365 } else if (*n0in == *n0 + 2) {
00366
00367
00368
00369
00370
00371 if (*dmin2 == *dn2 && z__[nn - 5] * 2. < z__[nn - 7]) {
00372 *ttype = -10;
00373 s = *dmin2 * .333;
00374 if (z__[nn - 5] > z__[nn - 7]) {
00375 return 0;
00376 }
00377 b1 = z__[nn - 5] / z__[nn - 7];
00378 b2 = b1;
00379 if (b2 == 0.) {
00380 goto L80;
00381 }
00382 i__1 = (*i0 << 2) - 1 + *pp;
00383 for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
00384 if (z__[i4] > z__[i4 - 2]) {
00385 return 0;
00386 }
00387 b1 *= z__[i4] / z__[i4 - 2];
00388 b2 += b1;
00389 if (b1 * 100. < b2) {
00390 goto L80;
00391 }
00392
00393 }
00394 L80:
00395 b2 = template_blas_sqrt(b2 * 1.05);
00396
00397 d__1 = b2;
00398 a2 = *dmin2 / (d__1 * d__1 + 1.);
00399 gap2 = z__[nn - 7] + z__[nn - 9] - template_blas_sqrt(z__[nn - 11]) * template_blas_sqrt(z__[
00400 nn - 9]) - a2;
00401 if (gap2 > 0. && gap2 > b2 * a2) {
00402
00403 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
00404 s = maxMACRO(d__1,d__2);
00405 } else {
00406
00407 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
00408 s = maxMACRO(d__1,d__2);
00409 }
00410 } else {
00411 s = *dmin2 * .25;
00412 *ttype = -11;
00413 }
00414 } else if (*n0in > *n0 + 2) {
00415
00416
00417
00418 s = 0.;
00419 *ttype = -12;
00420 }
00421
00422 *tau = s;
00423 return 0;
00424
00425
00426
00427 }
00428
00429 #endif