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