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_LASQ3_HEADER
00038 #define TEMPLATE_LAPACK_LASQ3_HEADER
00039
00040
00041 #include "template_lapack_lasq4.h"
00042 #include "template_lapack_lasq5.h"
00043 #include "template_lapack_lasq6.h"
00044
00045
00046 template<class Treal>
00047 int template_lapack_lasq3(integer *i0, integer *n0, Treal *z__,
00048 integer *pp, Treal *dmin__, Treal *sigma, Treal *desig,
00049 Treal *qmax, integer *nfail, integer *iter, integer *ndiv,
00050 logical *ieee, integer *ttype, Treal *dmin1, Treal *dmin2,
00051 Treal *dn, Treal *dn1, Treal *dn2, Treal *g,
00052 Treal *tau)
00053 {
00054
00055 integer i__1;
00056 Treal d__1, d__2;
00057
00058
00059
00060 Treal s, t;
00061 integer j4, nn;
00062 Treal eps, tol;
00063 integer n0in, ipn4;
00064 Treal tol2, temp;
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 --z__;
00154
00155
00156 n0in = *n0;
00157 eps = template_lapack_lamch("Precision", (Treal)0);
00158 tol = eps * 100.;
00159
00160 d__1 = tol;
00161 tol2 = d__1 * d__1;
00162
00163
00164
00165 L10:
00166
00167 if (*n0 < *i0) {
00168 return 0;
00169 }
00170 if (*n0 == *i0) {
00171 goto L20;
00172 }
00173 nn = (*n0 << 2) + *pp;
00174 if (*n0 == *i0 + 1) {
00175 goto L40;
00176 }
00177
00178
00179
00180 if (z__[nn - 5] > tol2 * (*sigma + z__[nn - 3]) && z__[nn - (*pp << 1) -
00181 4] > tol2 * z__[nn - 7]) {
00182 goto L30;
00183 }
00184
00185 L20:
00186
00187 z__[(*n0 << 2) - 3] = z__[(*n0 << 2) + *pp - 3] + *sigma;
00188 --(*n0);
00189 goto L10;
00190
00191
00192
00193 L30:
00194
00195 if (z__[nn - 9] > tol2 * *sigma && z__[nn - (*pp << 1) - 8] > tol2 * z__[
00196 nn - 11]) {
00197 goto L50;
00198 }
00199
00200 L40:
00201
00202 if (z__[nn - 3] > z__[nn - 7]) {
00203 s = z__[nn - 3];
00204 z__[nn - 3] = z__[nn - 7];
00205 z__[nn - 7] = s;
00206 }
00207 if (z__[nn - 5] > z__[nn - 3] * tol2) {
00208 t = (z__[nn - 7] - z__[nn - 3] + z__[nn - 5]) * .5;
00209 s = z__[nn - 3] * (z__[nn - 5] / t);
00210 if (s <= t) {
00211 s = z__[nn - 3] * (z__[nn - 5] / (t * (template_blas_sqrt(s / t + 1.) + 1.)));
00212 } else {
00213 s = z__[nn - 3] * (z__[nn - 5] / (t + template_blas_sqrt(t) * template_blas_sqrt(t + s)));
00214 }
00215 t = z__[nn - 7] + (s + z__[nn - 5]);
00216 z__[nn - 3] *= z__[nn - 7] / t;
00217 z__[nn - 7] = t;
00218 }
00219 z__[(*n0 << 2) - 7] = z__[nn - 7] + *sigma;
00220 z__[(*n0 << 2) - 3] = z__[nn - 3] + *sigma;
00221 *n0 += -2;
00222 goto L10;
00223
00224 L50:
00225 if (*pp == 2) {
00226 *pp = 0;
00227 }
00228
00229
00230
00231 if (*dmin__ <= 0. || *n0 < n0in) {
00232 if (z__[(*i0 << 2) + *pp - 3] * 1.5 < z__[(*n0 << 2) + *pp - 3]) {
00233 ipn4 = ( *i0 + *n0 ) << 2;
00234 i__1 = ( *i0 + *n0 - 1 ) << 1;
00235 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
00236 temp = z__[j4 - 3];
00237 z__[j4 - 3] = z__[ipn4 - j4 - 3];
00238 z__[ipn4 - j4 - 3] = temp;
00239 temp = z__[j4 - 2];
00240 z__[j4 - 2] = z__[ipn4 - j4 - 2];
00241 z__[ipn4 - j4 - 2] = temp;
00242 temp = z__[j4 - 1];
00243 z__[j4 - 1] = z__[ipn4 - j4 - 5];
00244 z__[ipn4 - j4 - 5] = temp;
00245 temp = z__[j4];
00246 z__[j4] = z__[ipn4 - j4 - 4];
00247 z__[ipn4 - j4 - 4] = temp;
00248
00249 }
00250 if (*n0 - *i0 <= 4) {
00251 z__[(*n0 << 2) + *pp - 1] = z__[(*i0 << 2) + *pp - 1];
00252 z__[(*n0 << 2) - *pp] = z__[(*i0 << 2) - *pp];
00253 }
00254
00255 d__1 = *dmin2, d__2 = z__[(*n0 << 2) + *pp - 1];
00256 *dmin2 = minMACRO(d__1,d__2);
00257
00258 d__1 = z__[(*n0 << 2) + *pp - 1], d__2 = z__[(*i0 << 2) + *pp - 1]
00259 , d__1 = minMACRO(d__1,d__2), d__2 = z__[(*i0 << 2) + *pp + 3];
00260 z__[(*n0 << 2) + *pp - 1] = minMACRO(d__1,d__2);
00261
00262 d__1 = z__[(*n0 << 2) - *pp], d__2 = z__[(*i0 << 2) - *pp], d__1 =
00263 minMACRO(d__1,d__2), d__2 = z__[(*i0 << 2) - *pp + 4];
00264 z__[(*n0 << 2) - *pp] = minMACRO(d__1,d__2);
00265
00266 d__1 = *qmax, d__2 = z__[(*i0 << 2) + *pp - 3], d__1 = maxMACRO(d__1,
00267 d__2), d__2 = z__[(*i0 << 2) + *pp + 1];
00268 *qmax = maxMACRO(d__1,d__2);
00269 *dmin__ = -0.;
00270 }
00271 }
00272
00273
00274
00275 template_lapack_lasq4(i0, n0, &z__[1], pp, &n0in, dmin__, dmin1, dmin2, dn, dn1, dn2,
00276 tau, ttype, g);
00277
00278
00279
00280 L70:
00281
00282 template_lapack_lasq5(i0, n0, &z__[1], pp, tau, dmin__, dmin1, dmin2, dn, dn1, dn2,
00283 ieee);
00284
00285 *ndiv += *n0 - *i0 + 2;
00286 ++(*iter);
00287
00288
00289
00290 if (*dmin__ >= 0. && *dmin1 > 0.) {
00291
00292
00293
00294 goto L90;
00295
00296 } else if (*dmin__ < 0. && *dmin1 > 0. && z__[( ( *n0 - 1 ) << 2) - *pp] < tol
00297 * (*sigma + *dn1) && absMACRO(*dn) < tol * *sigma) {
00298
00299
00300
00301 z__[( ( *n0 - 1 ) << 2) - *pp + 2] = 0.;
00302 *dmin__ = 0.;
00303 goto L90;
00304 } else if (*dmin__ < 0.) {
00305
00306
00307
00308 ++(*nfail);
00309 if (*ttype < -22) {
00310
00311
00312
00313 *tau = 0.;
00314 } else if (*dmin1 > 0.) {
00315
00316
00317
00318 *tau = (*tau + *dmin__) * (1. - eps * 2.);
00319 *ttype += -11;
00320 } else {
00321
00322
00323
00324 *tau *= .25;
00325 *ttype += -12;
00326 }
00327 goto L70;
00328 } else if (template_lapack_isnan(dmin__)) {
00329
00330
00331
00332 if (*tau == 0.) {
00333 goto L80;
00334 } else {
00335 *tau = 0.;
00336 goto L70;
00337 }
00338 } else {
00339
00340
00341
00342 goto L80;
00343 }
00344
00345
00346
00347 L80:
00348 template_lapack_lasq6(i0, n0, &z__[1], pp, dmin__, dmin1, dmin2, dn, dn1, dn2);
00349 *ndiv += *n0 - *i0 + 2;
00350 ++(*iter);
00351 *tau = 0.;
00352
00353 L90:
00354 if (*tau < *sigma) {
00355 *desig += *tau;
00356 t = *sigma + *desig;
00357 *desig -= t - *sigma;
00358 } else {
00359 t = *sigma + *tau;
00360 *desig = *sigma - (t - *tau) + *desig;
00361 }
00362 *sigma = t;
00363
00364 return 0;
00365
00366
00367
00368 }
00369
00370 #endif