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