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_LAMCH_HEADER
00038 #define TEMPLATE_LAPACK_LAMCH_HEADER
00039
00040
00041 #include <stdio.h>
00042 #include <iostream>
00043 #include <limits>
00044
00045
00046
00047 template<class Treal>
00048 Treal template_lapack_d_sign(const Treal *a, const Treal *b)
00049 {
00050 Treal x;
00051 x = (*a >= 0 ? *a : - *a);
00052 return( *b >= 0 ? x : -x);
00053 }
00054
00055
00056
00057
00058 #define log10e 0.43429448190325182765
00059 template<class Treal>
00060 Treal template_blas_lg10(Treal *x)
00061 {
00062 return( log10e * template_blas_log(*x) );
00063 }
00064
00065
00066
00067
00068
00069
00070
00071
00072 template<class Treal>
00073 int template_lapack_lassq(const integer *n, const Treal *x, const integer *incx,
00074 Treal *scale, Treal *sumsq)
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 integer i__1, i__2;
00129 Treal d__1;
00130
00131 Treal absxi;
00132 integer ix;
00133
00134 --x;
00135
00136
00137 if (*n > 0) {
00138 i__1 = (*n - 1) * *incx + 1;
00139 i__2 = *incx;
00140 for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) {
00141 if (x[ix] != 0.) {
00142 absxi = (d__1 = x[ix], absMACRO(d__1));
00143 if (*scale < absxi) {
00144
00145 d__1 = *scale / absxi;
00146 *sumsq = *sumsq * (d__1 * d__1) + 1;
00147 *scale = absxi;
00148 } else {
00149
00150 d__1 = absxi / *scale;
00151 *sumsq += d__1 * d__1;
00152 }
00153 }
00154
00155 }
00156 }
00157 return 0;
00158
00159
00160
00161 }
00162
00163
00164
00165
00166
00167 template<class Treal>
00168 double template_lapack_pow_di(Treal *ap, integer *bp)
00169 {
00170 Treal pow, x;
00171 integer n;
00172 unsigned long u;
00173
00174 pow = 1;
00175 x = *ap;
00176 n = *bp;
00177
00178 if(n != 0)
00179 {
00180 if(n < 0)
00181 {
00182 n = -n;
00183 x = 1/x;
00184 }
00185 for(u = n; ; )
00186 {
00187 if(u & 01)
00188 pow *= x;
00189 if(u >>= 1)
00190 x *= x;
00191 else
00192 break;
00193 }
00194 }
00195 return(pow);
00196 }
00197
00198
00199
00200
00201 template<class Treal>
00202 Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
00203 {
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 Treal rmach, ret_val;
00249
00250
00251 rmach = 0;
00252 if (template_blas_lsame(cmach, "E")) {
00253 rmach = template_blas_get_machine_epsilon<Treal>();
00254 } else if (template_blas_lsame(cmach, "S")) {
00255 rmach = template_blas_get_num_limit_min<Treal>();
00256 } else if (template_blas_lsame(cmach, "B")) {
00257
00258 rmach = 2.0;
00259 } else if (template_blas_lsame(cmach, "P")) {
00260
00261 rmach = 2.0 * template_blas_get_machine_epsilon<Treal>();
00262 } else if (template_blas_lsame(cmach, "N")) {
00263 std::cout << "ERROR in template_lapack_lamch: case N not implemented." << std::endl;
00264 throw "ERROR in template_lapack_lamch: case N not implemented.";
00265 } else if (template_blas_lsame(cmach, "R")) {
00266 std::cout << "ERROR in template_lapack_lamch: case R not implemented." << std::endl;
00267 throw "ERROR in template_lapack_lamch: case R not implemented.";
00268 } else if (template_blas_lsame(cmach, "M")) {
00269 std::cout << "ERROR in template_lapack_lamch: case M not implemented." << std::endl;
00270 throw "ERROR in template_lapack_lamch: case M not implemented.";
00271 } else if (template_blas_lsame(cmach, "U")) {
00272 std::cout << "ERROR in template_lapack_lamch: case U not implemented." << std::endl;
00273 throw "ERROR in template_lapack_lamch: case U not implemented.";
00274 } else if (template_blas_lsame(cmach, "L")) {
00275 std::cout << "ERROR in template_lapack_lamch: case L not implemented." << std::endl;
00276 throw "ERROR in template_lapack_lamch: case L not implemented.";
00277 } else if (template_blas_lsame(cmach, "O")) {
00278 std::cout << "ERROR in template_lapack_lamch: case O not implemented." << std::endl;
00279 throw "ERROR in template_lapack_lamch: case O not implemented.";
00280 }
00281
00282 ret_val = rmach;
00283 return ret_val;
00284
00285
00286
00287 }
00288
00289
00290
00291 #endif