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
00038
00039
00040 #ifndef _TEUCHOS_SCALARTRAITS_HPP_
00041 #define _TEUCHOS_SCALARTRAITS_HPP_
00042
00047 #include "Teuchos_ConfigDefs.hpp"
00048
00049 #ifndef HAVE_NUMERIC_LIMITS
00050 #include "Teuchos_LAPACK.hpp"
00051 #endif
00052
00053 #ifdef HAVE_TEUCHOS_ARPREC
00054 #include "arprec/mp_real.h"
00055 #endif
00056
00057 #ifdef HAVE_TEUCHOS_GNU_MP
00058 #include "gmp.h"
00059 #include "gmpxx.h"
00060 #endif
00061
00084
00085
00086
00087 namespace Teuchos {
00088
00089 template<class T>
00090 struct UndefinedScalarTraits
00091 {
00093 static inline T notDefined() { return T::this_type_is_missing_a_specialization(); }
00094 };
00095
00096 template<class T>
00097 struct ScalarTraits
00098 {
00100 typedef T magnitudeType;
00102 static const bool isComplex = false;
00104 static const bool isComparable = false;
00106 static const bool hasMachineParameters = false;
00108 static inline magnitudeType eps() { return UndefinedScalarTraits<T>::notDefined(); }
00110 static inline magnitudeType sfmin() { return UndefinedScalarTraits<T>::notDefined(); }
00112 static inline magnitudeType base() { return UndefinedScalarTraits<T>::notDefined(); }
00114 static inline magnitudeType prec() { return UndefinedScalarTraits<T>::notDefined(); }
00116 static inline magnitudeType t() { return UndefinedScalarTraits<T>::notDefined(); }
00118 static inline magnitudeType rnd() { return UndefinedScalarTraits<T>::notDefined(); }
00120 static inline magnitudeType emin() { return UndefinedScalarTraits<T>::notDefined(); }
00122 static inline magnitudeType rmin() { return UndefinedScalarTraits<T>::notDefined(); }
00124 static inline magnitudeType emax() { return UndefinedScalarTraits<T>::notDefined(); }
00126 static inline magnitudeType rmax() { return UndefinedScalarTraits<T>::notDefined(); }
00128 static inline magnitudeType magnitude(T a) { return UndefinedScalarTraits<T>::notDefined(); }
00130 static inline T zero() { return UndefinedScalarTraits<T>::notDefined(); }
00132 static inline T one() { return UndefinedScalarTraits<T>::notDefined(); }
00134 static inline magnitudeType real(T a) { return UndefinedScalarTraits<T>::notDefined(); }
00136 static inline magnitudeType imag(T a) { return UndefinedScalarTraits<T>::notDefined(); }
00138 static inline T conjugate(T a) { return UndefinedScalarTraits<T>::notDefined(); }
00140 static inline T nan() { return UndefinedScalarTraits<T>::notDefined(); }
00142 static inline bool isnaninf(const T& x) { return UndefinedScalarTraits<T>::notDefined(); }
00144 static inline void seedrandom(unsigned int s) { int i; T t = &i; }
00146 static inline T random() { return UndefinedScalarTraits<T>::notDefined(); }
00148 static inline std::string name() { (void)UndefinedScalarTraits<T>::notDefined(); return 0; }
00150 static inline T squareroot(T x) { return UndefinedScalarTraits<T>::notDefined(); }
00152 static inline T pow(T x, T y) { return UndefinedScalarTraits<T>::notDefined(); }
00153 };
00154
00155 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00156
00157
00158 void throwScalarTraitsNanInfError( const std::string &errMsg );
00159
00160
00161 #define TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR( VALUE, MSG ) \
00162 if (isnaninf(VALUE)) { \
00163 std::ostringstream omsg; \
00164 omsg << MSG; \
00165 throwScalarTraitsNanInfError(omsg.str()); \
00166 }
00167
00168
00169 template<>
00170 struct ScalarTraits<char>
00171 {
00172 typedef char magnitudeType;
00173 static const bool isComplex = false;
00174 static const bool isComparable = true;
00175 static const bool hasMachineParameters = false;
00176
00177 static inline magnitudeType magnitude(char a) { return static_cast<char>(std::fabs(static_cast<double>(a))); }
00178 static inline char zero() { return 0; }
00179 static inline char one() { return 1; }
00180 static inline char conjugate(char x) { return x; }
00181 static inline char real(char x) { return x; }
00182 static inline char imag(char) { return 0; }
00183 static inline bool isnaninf(char ) { return false; }
00184 static inline void seedrandom(unsigned int s) {
00185 std::srand(s);
00186 #ifdef __APPLE__
00187
00188
00189 random();
00190 #endif
00191 }
00192
00193 static inline char random() { return std::rand(); }
00194 static inline std::string name() { return "char"; }
00195 static inline char squareroot(char x) { return (char) std::sqrt((double) x); }
00196 static inline char pow(char x, char y) { return (char) std::pow((double)x,(double)y); }
00197 };
00198
00199 template<>
00200 struct ScalarTraits<short int>
00201 {
00202 typedef short int magnitudeType;
00203 static const bool isComplex = false;
00204 static const bool isComparable = true;
00205 static const bool hasMachineParameters = false;
00206
00207 static inline magnitudeType magnitude(short int a) { return static_cast<short int>(std::fabs(static_cast<double>(a))); }
00208 static inline short int zero() { return 0; }
00209 static inline short int one() { return 1; }
00210 static inline short int conjugate(short int x) { return x; }
00211 static inline short int real(short int x) { return x; }
00212 static inline short int imag(short int) { return 0; }
00213 static inline void seedrandom(unsigned int s) {
00214 std::srand(s);
00215 #ifdef __APPLE__
00216
00217
00218 random();
00219 #endif
00220 }
00221
00222 static inline short int random() { return std::rand(); }
00223 static inline std::string name() { return "short int"; }
00224 static inline short int squareroot(short int x) { return (short int) std::sqrt((double) x); }
00225 static inline short int pow(short int x, short int y) { return (short int) std::pow((double)x,(double)y); }
00226 };
00227
00228 template<>
00229 struct ScalarTraits<int>
00230 {
00231 typedef int magnitudeType;
00232 static const bool isComplex = false;
00233 static const bool isComparable = true;
00234 static const bool hasMachineParameters = false;
00235
00236 static inline magnitudeType magnitude(int a) { return static_cast<int>(std::fabs(static_cast<double>(a))); }
00237 static inline int zero() { return 0; }
00238 static inline int one() { return 1; }
00239 static inline int conjugate(int x) { return x; }
00240 static inline int real(int x) { return x; }
00241 static inline int imag(int) { return 0; }
00242 static inline bool isnaninf(int) { return false; }
00243 static inline void seedrandom(unsigned int s) {
00244 std::srand(s);
00245 #ifdef __APPLE__
00246
00247
00248 random();
00249 #endif
00250 }
00251
00252 static inline int random() { return std::rand(); }
00253 static inline std::string name() { return "int"; }
00254 static inline int squareroot(int x) { return (int) std::sqrt((double) x); }
00255 static inline int pow(int x, int y) { return (int) std::pow((double)x,(double)y); }
00256 };
00257
00258 template<>
00259 struct ScalarTraits<long int>
00260 {
00261 typedef long int magnitudeType;
00262 static const bool isComplex = false;
00263 static const bool isComparable = true;
00264 static const bool hasMachineParameters = false;
00265
00266 static inline magnitudeType magnitude(long int a) { return static_cast<long int>(std::fabs(static_cast<double>(a))); }
00267 static inline long int zero() { return 0; }
00268 static inline long int one() { return 1; }
00269 static inline long int conjugate(long int x) { return x; }
00270 static inline long int real(long int x) { return x; }
00271 static inline long int imag(long int) { return 0; }
00272 static inline void seedrandom(unsigned int s) {
00273 std::srand(s);
00274 #ifdef __APPLE__
00275
00276
00277 random();
00278 #endif
00279 }
00280
00281 static inline long int random() { return std::rand(); }
00282 static inline std::string name() { return "long int"; }
00283 static inline long int squareroot(long int x) { return (long int) std::sqrt((double) x); }
00284 static inline long int pow(long int x, long int y) { return (long int) std::pow((double)x,(double)y); }
00285 };
00286
00287 #ifndef __sun
00288 extern const float flt_nan;
00289 #endif
00290
00291 template<>
00292 struct ScalarTraits<float>
00293 {
00294 typedef float magnitudeType;
00295 static const bool isComplex = false;
00296 static const bool isComparable = true;
00297 static const bool hasMachineParameters = true;
00298 static inline float eps() {
00299 #ifdef HAVE_NUMERIC_LIMITS
00300 return std::numeric_limits<float>::epsilon();
00301 #else
00302 LAPACK<int, float> lp; return lp.LAMCH('E');
00303 #endif
00304 }
00305 static inline float sfmin() {
00306 #ifdef HAVE_NUMERIC_LIMITS
00307 return std::numeric_limits<float>::min();
00308 #else
00309 LAPACK<int, float> lp; return lp.LAMCH('S');
00310 #endif
00311 }
00312 static inline float base() {
00313 #ifdef HAVE_NUMERIC_LIMITS
00314 return std::numeric_limits<float>::radix;
00315 #else
00316 LAPACK<int, float> lp; return lp.LAMCH('B');
00317 #endif
00318 }
00319 static inline float prec() {
00320 #ifdef HAVE_NUMERIC_LIMITS
00321 return eps()*base();
00322 #else
00323 LAPACK<int, float> lp; return lp.LAMCH('P');
00324 #endif
00325 }
00326 static inline float t() {
00327 #ifdef HAVE_NUMERIC_LIMITS
00328 return std::numeric_limits<float>::digits;
00329 #else
00330 LAPACK<int, float> lp; return lp.LAMCH('N');
00331 #endif
00332 }
00333 static inline float rnd() {
00334 #ifdef HAVE_NUMERIC_LIMITS
00335 return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? float(1.0) : float(0.0) );
00336 #else
00337 LAPACK<int, float> lp; return lp.LAMCH('R');
00338 #endif
00339 }
00340 static inline float emin() {
00341 #ifdef HAVE_NUMERIC_LIMITS
00342 return std::numeric_limits<float>::min_exponent;
00343 #else
00344 LAPACK<int, float> lp; return lp.LAMCH('M');
00345 #endif
00346 }
00347 static inline float rmin() {
00348 #ifdef HAVE_NUMERIC_LIMITS
00349 return std::numeric_limits<float>::min();
00350 #else
00351 LAPACK<int, float> lp; return lp.LAMCH('U');
00352 #endif
00353 }
00354 static inline float emax() {
00355 #ifdef HAVE_NUMERIC_LIMITS
00356 return std::numeric_limits<float>::max_exponent;
00357 #else
00358 LAPACK<int, float> lp; return lp.LAMCH('L');
00359 #endif
00360 }
00361 static inline float rmax() {
00362 #ifdef HAVE_NUMERIC_LIMITS
00363 return std::numeric_limits<float>::max();
00364 #else
00365 LAPACK<int, float> lp; return lp.LAMCH('O');
00366 #endif
00367 }
00368 static inline magnitudeType magnitude(float a)
00369 {
00370 #ifdef TEUCHOS_DEBUG
00371 TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00372 a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00373 #endif
00374 return std::fabs(a);
00375 }
00376 static inline float zero() { return(0.0); }
00377 static inline float one() { return(1.0); }
00378 static inline float conjugate(float x) { return(x); }
00379 static inline float real(float x) { return x; }
00380 static inline float imag(float) { return 0; }
00381 static inline float nan() {
00382 #ifdef __sun
00383 return 0.0/std::sin(0.0);
00384 #else
00385 return flt_nan;
00386 #endif
00387 }
00388 static inline bool isnaninf(float x) {
00389 const float tol = 1e-6;
00390 if( !(x <= tol) && !(x > tol) ) return true;
00391 float z=0.0*x; if( !(z <= tol) && !(z > tol) ) return true;
00392 return false;
00393 }
00394 static inline void seedrandom(unsigned int s) {
00395 std::srand(s);
00396 #ifdef __APPLE__
00397
00398
00399 random();
00400 #endif
00401 }
00402 static inline float random() { float rnd = (float) std::rand() / RAND_MAX; return (float)(-1.0 + 2.0 * rnd); }
00403 static inline std::string name() { return "float"; }
00404 static inline float squareroot(float x)
00405 {
00406 #ifdef TEUCHOS_DEBUG
00407 TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00408 x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00409 #endif
00410 errno = 0;
00411 const float rtn = std::sqrt(x);
00412 if (errno)
00413 return nan();
00414 return rtn;
00415 }
00416 static inline float pow(float x, float y) { return std::pow(x,y); }
00417 };
00418
00419 #ifndef __sun
00420 extern const double dbl_nan;
00421 #endif
00422
00423 template<>
00424 struct ScalarTraits<double>
00425 {
00426 typedef double magnitudeType;
00427 static const bool isComplex = false;
00428 static const bool isComparable = true;
00429 static const bool hasMachineParameters = true;
00430 static inline double eps() {
00431 #ifdef HAVE_NUMERIC_LIMITS
00432 return std::numeric_limits<double>::epsilon();
00433 #else
00434 LAPACK<int, double> lp; return lp.LAMCH('E');
00435 #endif
00436 }
00437 static inline double sfmin() {
00438 #ifdef HAVE_NUMERIC_LIMITS
00439 return std::numeric_limits<double>::min();
00440 #else
00441 LAPACK<int, double> lp; return lp.LAMCH('S');
00442 #endif
00443 }
00444 static inline double base() {
00445 #ifdef HAVE_NUMERIC_LIMITS
00446 return std::numeric_limits<double>::radix;
00447 #else
00448 LAPACK<int, double> lp; return lp.LAMCH('B');
00449 #endif
00450 }
00451 static inline double prec() {
00452 #ifdef HAVE_NUMERIC_LIMITS
00453 return eps()*base();
00454 #else
00455 LAPACK<int, double> lp; return lp.LAMCH('P');
00456 #endif
00457 }
00458 static inline double t() {
00459 #ifdef HAVE_NUMERIC_LIMITS
00460 return std::numeric_limits<double>::digits;
00461 #else
00462 LAPACK<int, double> lp; return lp.LAMCH('N');
00463 #endif
00464 }
00465 static inline double rnd() {
00466 #ifdef HAVE_NUMERIC_LIMITS
00467 return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
00468 #else
00469 LAPACK<int, double> lp; return lp.LAMCH('R');
00470 #endif
00471 }
00472 static inline double emin() {
00473 #ifdef HAVE_NUMERIC_LIMITS
00474 return std::numeric_limits<double>::min_exponent;
00475 #else
00476 LAPACK<int, double> lp; return lp.LAMCH('M');
00477 #endif
00478 }
00479 static inline double rmin() {
00480 #ifdef HAVE_NUMERIC_LIMITS
00481 return std::numeric_limits<double>::min();
00482 #else
00483 LAPACK<int, double> lp; return lp.LAMCH('U');
00484 #endif
00485 }
00486 static inline double emax() {
00487 #ifdef HAVE_NUMERIC_LIMITS
00488 return std::numeric_limits<double>::max_exponent;
00489 #else
00490 LAPACK<int, double> lp; return lp.LAMCH('L');
00491 #endif
00492 }
00493 static inline double rmax() {
00494 #ifdef HAVE_NUMERIC_LIMITS
00495 return std::numeric_limits<double>::max();
00496 #else
00497 LAPACK<int, double> lp; return lp.LAMCH('O');
00498 #endif
00499 }
00500 static inline magnitudeType magnitude(double a)
00501 {
00502 #ifdef TEUCHOS_DEBUG
00503 TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00504 a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00505 #endif
00506 return std::fabs(a);
00507 }
00508 static inline double zero() { return 0.0; }
00509 static inline double one() { return 1.0; }
00510 static inline double conjugate(double x) { return(x); }
00511 static inline double real(double x) { return(x); }
00512 static inline double imag(double) { return(0); }
00513 static inline double nan() {
00514 #ifdef __sun
00515 return 0.0/std::sin(0.0);
00516 #else
00517 return dbl_nan;
00518 #endif
00519 }
00520 static inline bool isnaninf(double x) {
00521 const double tol = 1e-6;
00522 if( !(x <= tol) && !(x > tol) ) return true;
00523 double z=0.0*x; if( !(z <= tol) && !(z > tol) ) return true;
00524 return false;
00525 }
00526 static inline void seedrandom(unsigned int s) {
00527 std::srand(s);
00528 #ifdef __APPLE__
00529
00530
00531 random();
00532 #endif
00533 }
00534 static inline double random() { double rnd = (double) std::rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); }
00535 static inline std::string name() { return "double"; }
00536 static inline double squareroot(double x)
00537 {
00538 #ifdef TEUCHOS_DEBUG
00539 TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00540 x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00541 #endif
00542 errno = 0;
00543 const double rtn = std::sqrt(x);
00544 if (errno)
00545 return nan();
00546 return rtn;
00547 }
00548 static inline double pow(double x, double y) { return std::pow(x,y); }
00549 };
00550
00551 #ifdef HAVE_TEUCHOS_GNU_MP
00552
00553 extern gmp_randclass gmp_rng;
00554
00555 template<>
00556 struct ScalarTraits<mpf_class>
00557 {
00558 typedef mpf_class magnitudeType;
00559 static const bool isComplex = false;
00560 static const bool isComparable = true;
00561 static const bool hasMachineParameters = false;
00562
00563 static magnitudeType magnitude(mpf_class a) { return abs(a); }
00564 static inline mpf_class zero() { mpf_class zero = 0.0; return zero; }
00565 static inline mpf_class one() { mpf_class one = 1.0; return one; }
00566 static inline mpf_class conjugate(mpf_class x) { return x; }
00567 static inline mpf_class real(mpf_class x) { return(x); }
00568 static inline mpf_class imag(mpf_class x) { return(0); }
00569 static inline bool isnaninf(mpf_class x) { return false; }
00570 static inline void seedrandom(unsigned int s) {
00571 unsigned long int seedVal = static_cast<unsigned long int>(s);
00572 gmp_rng.seed( seedVal );
00573 }
00574 static inline mpf_class random() {
00575 return gmp_rng.get_f();
00576 }
00577 static inline std::string name() { return "mpf_class"; }
00578 static inline mpf_class squareroot(mpf_class x) { return sqrt(x); }
00579 static inline mpf_class pow(mpf_class x, mpf_class y) { return pow(x,y); }
00580
00581 };
00582
00583 #endif
00584
00585 #ifdef HAVE_TEUCHOS_ARPREC
00586
00587 template<>
00588 struct ScalarTraits<mp_real>
00589 {
00590 typedef mp_real magnitudeType;
00591 static const bool isComplex = false;
00592 static const bool isComparable = true;
00593 static const bool hasMachineParameters = false;
00594
00595 static magnitudeType magnitude(mp_real a) { return abs(a); }
00596 static inline mp_real zero() { mp_real zero = 0.0; return zero; }
00597 static inline mp_real one() { mp_real one = 1.0; return one; }
00598 static inline mp_real conjugate(mp_real x) { return x; }
00599 static inline mp_real real(mp_real x) { return(x); }
00600 static inline mp_real imag(mp_real x) { return(0.); }
00601 static inline bool isnaninf(mp_real x) { return false; }
00602 static inline void seedrandom(unsigned int s) {
00603 long int seedVal = static_cast<long int>(s);
00604 srand48(seedVal);
00605 }
00606 static inline mp_real random() { return mp_rand(); }
00607 static inline std::string name() { return "mp_real"; }
00608 static inline mp_real squareroot(mp_real x) { return sqrt(x); }
00609 static inline mp_real pow(mp_real x, mp_real y) { return pow(x,y); }
00610
00611 };
00612
00613 #endif // HAVE_TEUCHOS_ARPREC
00614
00615 #if ( defined(HAVE_COMPLEX) || defined(HAVE_COMPLEX_H) ) && defined(HAVE_TEUCHOS_COMPLEX)
00616
00617
00618 template<class T>
00619 struct ScalarTraits<
00620 #if defined(HAVE_COMPLEX)
00621 std::complex<T>
00622 #elif defined(HAVE_COMPLEX_H)
00623 std::complex<T>
00624 #endif
00625 >
00626 {
00627 #if defined(HAVE_COMPLEX)
00628 typedef std::complex<T> ComplexT;
00629 #elif defined(HAVE_COMPLEX_H)
00630 typedef std::complex<T> ComplexT;
00631 #endif
00632 typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
00633 static const bool isComplex = true;
00634 static const bool isComparable = false;
00635 static const bool hasMachineParameters = true;
00636 static inline magnitudeType eps() { return ScalarTraits<magnitudeType>::eps(); }
00637 static inline magnitudeType sfmin() { return ScalarTraits<magnitudeType>::sfmin(); }
00638 static inline magnitudeType base() { return ScalarTraits<magnitudeType>::base(); }
00639 static inline magnitudeType prec() { return ScalarTraits<magnitudeType>::prec(); }
00640 static inline magnitudeType t() { return ScalarTraits<magnitudeType>::t(); }
00641 static inline magnitudeType rnd() { return ScalarTraits<magnitudeType>::rnd(); }
00642 static inline magnitudeType emin() { return ScalarTraits<magnitudeType>::emin(); }
00643 static inline magnitudeType rmin() { return ScalarTraits<magnitudeType>::rmin(); }
00644 static inline magnitudeType emax() { return ScalarTraits<magnitudeType>::emax(); }
00645 static inline magnitudeType rmax() { return ScalarTraits<magnitudeType>::rmax(); }
00646 static magnitudeType magnitude(ComplexT a)
00647 {
00648 #ifdef TEUCHOS_DEBUG
00649 TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00650 a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00651 #endif
00652 return std::abs(a);
00653 }
00654 static inline ComplexT zero() { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
00655 static inline ComplexT one() { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
00656 static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
00657 static inline magnitudeType real(ComplexT a) { return a.real(); }
00658 static inline magnitudeType imag(ComplexT a) { return a.imag(); }
00659 static inline ComplexT nan() { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
00660 static inline bool isnaninf(ComplexT x) { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
00661 static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
00662 static inline ComplexT random()
00663 {
00664 const T rnd1 = ScalarTraits<magnitudeType>::random();
00665 const T rnd2 = ScalarTraits<magnitudeType>::random();
00666 return ComplexT(rnd1,rnd2);
00667 }
00668 static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
00669
00670 static inline ComplexT squareroot(ComplexT x)
00671 {
00672 #ifdef TEUCHOS_DEBUG
00673 TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00674 x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00675 #endif
00676 typedef ScalarTraits<magnitudeType> STMT;
00677 const T r = x.real(), i = x.imag();
00678 const T a = STMT::squareroot((r*r)+(i*i));
00679 const T nr = STMT::squareroot((a+r)/2);
00680 const T ni = STMT::squareroot((a-r)/2);
00681 return ComplexT(nr,ni);
00682 }
00683 static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(x,y); }
00684 };
00685
00686 #endif // HAVE_COMPLEX || HAVE_COMPLEX_H
00687
00688 #endif // DOXYGEN_SHOULD_SKIP_THIS
00689
00690 }
00691
00692 #endif // _TEUCHOS_SCALARTRAITS_HPP_