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_LAGTS_HEADER
00038 #define TEMPLATE_LAPACK_LAGTS_HEADER
00039
00040
00041 template<class Treal>
00042 int template_lapack_lagts(const integer *job, const integer *n, const Treal *a,
00043 const Treal *b, const Treal *c__, const Treal *d__, const integer *in,
00044 Treal *y, Treal *tol, integer *info)
00045 {
00046
00047
00048
00049
00050
00051
00052
00053
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
00139
00140 integer i__1;
00141 Treal d__1, d__2, d__3, d__4, d__5;
00142
00143 Treal temp, pert;
00144 integer k;
00145 Treal absak, sfmin, ak;
00146 Treal bignum, eps;
00147
00148 --y;
00149 --in;
00150 --d__;
00151 --c__;
00152 --b;
00153 --a;
00154
00155
00156 *info = 0;
00157 if (absMACRO(*job) > 2 || *job == 0) {
00158 *info = -1;
00159 } else if (*n < 0) {
00160 *info = -2;
00161 }
00162 if (*info != 0) {
00163 i__1 = -(*info);
00164 template_blas_erbla("LAGTS ", &i__1);
00165 return 0;
00166 }
00167
00168 if (*n == 0) {
00169 return 0;
00170 }
00171
00172 eps = template_lapack_lamch("Epsilon", (Treal)0);
00173 sfmin = template_lapack_lamch("Safe minimum", (Treal)0);
00174 bignum = 1. / sfmin;
00175
00176 if (*job < 0) {
00177 if (*tol <= 0.) {
00178 *tol = absMACRO(a[1]);
00179 if (*n > 1) {
00180
00181 d__1 = *tol, d__2 = absMACRO(a[2]), d__1 = maxMACRO(d__1,d__2), d__2 =
00182 absMACRO(b[1]);
00183 *tol = maxMACRO(d__1,d__2);
00184 }
00185 i__1 = *n;
00186 for (k = 3; k <= i__1; ++k) {
00187
00188 d__4 = *tol, d__5 = (d__1 = a[k], absMACRO(d__1)), d__4 = maxMACRO(d__4,
00189 d__5), d__5 = (d__2 = b[k - 1], absMACRO(d__2)), d__4 =
00190 maxMACRO(d__4,d__5), d__5 = (d__3 = d__[k - 2], absMACRO(d__3));
00191 *tol = maxMACRO(d__4,d__5);
00192
00193 }
00194 *tol *= eps;
00195 if (*tol == 0.) {
00196 *tol = eps;
00197 }
00198 }
00199 }
00200
00201 if (absMACRO(*job) == 1) {
00202 i__1 = *n;
00203 for (k = 2; k <= i__1; ++k) {
00204 if (in[k - 1] == 0) {
00205 y[k] -= c__[k - 1] * y[k - 1];
00206 } else {
00207 temp = y[k - 1];
00208 y[k - 1] = y[k];
00209 y[k] = temp - c__[k - 1] * y[k];
00210 }
00211
00212 }
00213 if (*job == 1) {
00214 for (k = *n; k >= 1; --k) {
00215 if (k <= *n - 2) {
00216 temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
00217 } else if (k == *n - 1) {
00218 temp = y[k] - b[k] * y[k + 1];
00219 } else {
00220 temp = y[k];
00221 }
00222 ak = a[k];
00223 absak = absMACRO(ak);
00224 if (absak < 1.) {
00225 if (absak < sfmin) {
00226 if (absak == 0. || absMACRO(temp) * sfmin > absak) {
00227 *info = k;
00228 return 0;
00229 } else {
00230 temp *= bignum;
00231 ak *= bignum;
00232 }
00233 } else if (absMACRO(temp) > absak * bignum) {
00234 *info = k;
00235 return 0;
00236 }
00237 }
00238 y[k] = temp / ak;
00239
00240 }
00241 } else {
00242 for (k = *n; k >= 1; --k) {
00243 if (k <= *n - 2) {
00244 temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
00245 } else if (k == *n - 1) {
00246 temp = y[k] - b[k] * y[k + 1];
00247 } else {
00248 temp = y[k];
00249 }
00250 ak = a[k];
00251 pert = template_lapack_d_sign(tol, &ak);
00252 L40:
00253 absak = absMACRO(ak);
00254 if (absak < 1.) {
00255 if (absak < sfmin) {
00256 if (absak == 0. || absMACRO(temp) * sfmin > absak) {
00257 ak += pert;
00258 pert *= 2;
00259 goto L40;
00260 } else {
00261 temp *= bignum;
00262 ak *= bignum;
00263 }
00264 } else if (absMACRO(temp) > absak * bignum) {
00265 ak += pert;
00266 pert *= 2;
00267 goto L40;
00268 }
00269 }
00270 y[k] = temp / ak;
00271
00272 }
00273 }
00274 } else {
00275
00276
00277
00278 if (*job == 2) {
00279 i__1 = *n;
00280 for (k = 1; k <= i__1; ++k) {
00281 if (k >= 3) {
00282 temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
00283 } else if (k == 2) {
00284 temp = y[k] - b[k - 1] * y[k - 1];
00285 } else {
00286 temp = y[k];
00287 }
00288 ak = a[k];
00289 absak = absMACRO(ak);
00290 if (absak < 1.) {
00291 if (absak < sfmin) {
00292 if (absak == 0. || absMACRO(temp) * sfmin > absak) {
00293 *info = k;
00294 return 0;
00295 } else {
00296 temp *= bignum;
00297 ak *= bignum;
00298 }
00299 } else if (absMACRO(temp) > absak * bignum) {
00300 *info = k;
00301 return 0;
00302 }
00303 }
00304 y[k] = temp / ak;
00305
00306 }
00307 } else {
00308 i__1 = *n;
00309 for (k = 1; k <= i__1; ++k) {
00310 if (k >= 3) {
00311 temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
00312 } else if (k == 2) {
00313 temp = y[k] - b[k - 1] * y[k - 1];
00314 } else {
00315 temp = y[k];
00316 }
00317 ak = a[k];
00318 pert = template_lapack_d_sign(tol, &ak);
00319 L70:
00320 absak = absMACRO(ak);
00321 if (absak < 1.) {
00322 if (absak < sfmin) {
00323 if (absak == 0. || absMACRO(temp) * sfmin > absak) {
00324 ak += pert;
00325 pert *= 2;
00326 goto L70;
00327 } else {
00328 temp *= bignum;
00329 ak *= bignum;
00330 }
00331 } else if (absMACRO(temp) > absak * bignum) {
00332 ak += pert;
00333 pert *= 2;
00334 goto L70;
00335 }
00336 }
00337 y[k] = temp / ak;
00338
00339 }
00340 }
00341
00342 for (k = *n; k >= 2; --k) {
00343 if (in[k - 1] == 0) {
00344 y[k - 1] -= c__[k - 1] * y[k];
00345 } else {
00346 temp = y[k - 1];
00347 y[k - 1] = y[k];
00348 y[k] = temp - c__[k - 1] * y[k];
00349 }
00350
00351 }
00352 }
00353
00354
00355
00356 return 0;
00357 }
00358
00359 #endif