1 /////////////////////////////////////////////////////////////////////////////// 2 // // 3 // The Template Matrix/Vector Library for C++ was created by Mike Jarvis // 4 // Copyright (C) 1998 - 2016 // 5 // All rights reserved // 6 // // 7 // The project is hosted at https://code.google.com/p/tmv-cpp/ // 8 // where you can find the current version and current documention. // 9 // // 10 // For concerns or problems with the software, Mike may be contacted at // 11 // mike_jarvis17 [at] gmail. // 12 // // 13 // This software is licensed under a FreeBSD license. The file // 14 // TMV_LICENSE should have bee included with this distribution. // 15 // It not, you can get a copy from https://code.google.com/p/tmv-cpp/. // 16 // // 17 // Essentially, you can use this software however you want provided that // 18 // you include the TMV_LICENSE file in any distribution that uses it. // 19 // // 20 /////////////////////////////////////////////////////////////////////////////// 21 22 23 //----------------------------------------------------------------------------- 24 // 25 // This file defines some basic helper functions used by the TMV Matrix 26 // and Vector classes. 27 // 28 // Probably the only thing someone might care about here is the 29 // class tmv::tmv_exception which gets thrown on run-time errors, 30 // so you could catch it if you don't want the program to crash. 31 // 32 // Also, it is a good idea to compile any code that uses TMV with 33 // TMV_DEBUG defined. (defined below in this file) 34 // 35 // This will turn on some basic debugging in the code and will help 36 // you catch programming mistakes. 37 // For example, if you have two vectors: 38 // v1 of length 5 and v2 of lengths 10, 39 // then v1 = v2 will only produce a runtime error if TMV_DEBUG is 40 // turned on. Otherwise some random data will be overwritten and 41 // who knows what will happen then. 42 // 43 // Once you know your code is working properly, you can recompile without 44 // the TMV_DEBUG flag to speed up the code. 45 // (And actually, the slow down isn't much, so it is probably 46 // a good idea to just leave it on to help diagnose any bug that 47 // might turn up.) 48 // 49 // A comment here about the letters used to identify each type of Matrix. 50 // The writeCompact routines for sparse Matrix types have an identifying 51 // letter(s) before the output. The same letter (or combination) is 52 // used in the names of the composite arithmetic types like SumDD, etc. 53 // Here is the list of letters used for reference: 54 // 55 // (* indicates that the type is not yet implemented.) 56 // 57 // X Scalar 58 // V Vector 59 // M Matrix 60 // D DiagMatrix 61 // U UpperTriMatrix 62 // L LowerTriMatrix 63 // B BandMatrix 64 // S SymMatrix 65 // H HermMatrix 66 // v SmallVector 67 // m SmallMatrix 68 // sB SymBandMatrix 69 // hB HermBandMatrix 70 // *kD BlockDiagMatrix 71 // *skD SymBlockDiagMatrix 72 // *hkD HermBlockDiagMatrix 73 // *W SparseVector 74 // *Q SparseMatrix 75 // *sQ SymSparseMatrix 76 // *hQ HermSparseMatrix 77 // *kQ BlockSparseMatrix 78 // *skQ SymBlockSparseMatrix 79 // *hkQ HermBlockSparseMatrix 80 // 81 // TODO BlockDiagMatrix and Sym varieties 82 // TODO SparseMatrix and Sym, Block varieties 83 84 85 #ifndef TMV_Base_H 86 #define TMV_Base_H 87 88 #if ((!defined(NDEBUG) && !defined(TMV_NDEBUG)) || defined(TMV_EXTRA_DEBUG)) 89 #define TMV_DEBUG 90 #endif 91 92 #ifndef TMV_NO_WARN 93 #define TMV_WARN 94 #endif 95 96 #include <iosfwd> 97 #include <limits> 98 #include <cmath> 99 #include <complex> 100 #include <memory> 101 #include <stdexcept> 102 #include <cstddef> 103 #include <string> 104 #include <algorithm> 105 106 #if defined(__SSE2__) || defined(__SSE__) 107 #include "xmmintrin.h" 108 #endif 109 110 #ifdef TMV_DEBUG 111 #include <typeinfo> 112 #include <iostream> 113 #endif 114 115 #include <typeinfo> 116 117 #ifdef TMV_MEM_DEBUG 118 #include "tmv/mmgr.h" 119 #endif 120 121 #ifndef NO_INST_DOUBLE 122 #define INST_DOUBLE 123 #endif 124 125 #ifndef NO_INST_FLOAT 126 #define INST_FLOAT 127 #endif 128 129 #ifndef NO_INST_COMPLEX 130 #define INST_COMPLEX 131 #endif 132 133 namespace tmv { 134 TMV_Version()135 inline std::string TMV_Version() { return "0.75"; } 136 #define TMV_MAJOR_VERSION 0 137 #define TMV_MINOR_VERSION 75 138 #define TMV_VERSION_AT_LEAST(major,minor) \ 139 ( (TMV_MAJOR_VERSION > major) || \ 140 (TMV_MAJOR_VERSION == major && TMV_MINOR_VERSION >= minor) ) 141 142 // Attributes are stored as a binary field, so they can be |'ed 143 // together to calculate the full set of attributes. 144 // 145 // For each object we have a default set of attributes that is 146 // implied by A=0. This way we can write, for example: 147 // UpperTriMatrix<T,UnitDiag> U; 148 // and this will imply (UnitDiag | ColMajor | CStyle). 149 // 150 // Attributes that are always the default (whenever they are allowed) 151 // can be defined as 0, since the absence of the converse is sufficient. 152 // e.g. CStyle, ColMajor, Lower 153 154 enum IndexStyle { CStyle = 0, FortranStyle = 0x1 }; 155 enum StorageType { 156 ColMajor = 0, RowMajor = 0x2, DiagMajor = 0x4, 157 AllStorageType = 0x6 }; 158 enum DiagType { NonUnitDiag = 0, UnitDiag = 0x8 }; 159 enum UpLoType { Lower = 0, Upper = 0x10 }; 160 161 template <int A> 162 struct Attrib 163 { 164 enum { onestoragetype = !( !!(A & RowMajor) && !!(A & DiagMajor) ) }; 165 enum { viewok = A <= 0x1 }; 166 enum { vectorok = A <= 0x1 }; 167 enum { matrixok = A <= 0x3 }; 168 enum { diagmatrixok = A <= 0x1 }; 169 enum { trimatrixok = ( A <= 0xb && !(A & DiagMajor) ) }; 170 enum { bandmatrixok = A <= 0x5 }; 171 enum { symmatrixok = ( 172 A <= 0x13 && !(A & DiagMajor) && !(A & UnitDiag) ) }; 173 enum { symbandmatrixok = ( 174 A <= 0x15 && onestoragetype && !(A & UnitDiag) ) }; 175 176 enum { fort = !!(A & FortranStyle) }; 177 enum { rowmajor = !!(A & RowMajor) }; 178 enum { diagmajor = !!(A & DiagMajor) }; 179 enum { colmajor = !(rowmajor || diagmajor) }; 180 enum { unitdiag = !!(A & UnitDiag) }; 181 enum { nonunitdiag = !unitdiag }; 182 enum { upper = !!(A & Upper) }; 183 enum { lower = !upper }; 184 185 enum { stor = rowmajor ? RowMajor : diagmajor ? DiagMajor : ColMajor }; 186 textAttrib187 static std::string text() 188 { 189 std::string ret = 190 ( colmajor ? "ColMajor" : 191 rowmajor ? "RowMajor" : "DiagMajor" ); 192 ret += fort ? "|FortranStyle" : ""; 193 ret += unitdiag ? "|UnitDiag" : ""; 194 ret += upper ? "|Upper" : ""; 195 return ret; 196 } vtextAttrib197 static std::string vtext() 198 { return fort ? "FortranStyle" : "CStyle"; } 199 }; 200 201 enum ADType { Ascend, Descend }; 202 enum CompType { RealComp, AbsComp, ImagComp, ArgComp }; 203 enum StepType { Unit, Step }; 204 enum ConjType { NonConj, Conj }; 205 enum SymType { Sym, Herm }; 206 207 #define TMV_UTransOf(U) (U==int(Upper) ? Lower : Upper) 208 209 enum DivType { 210 XX=0, LU=1, CH=2, QR=4, QRP=8, SV=16, 211 // We store the divtype in a binary field integer. 212 // In addition to the above, we also use the same object to 213 // store the following other flags related to division. 214 // So these values must not clash with the above DivType values. 215 // These aren't technically DivType's but since they are 216 // stored together, I think this adds to type-safety. 217 DivInPlaceFlag = 32, 218 SaveDivFlag = 64, 219 // And finally shorthand for "one of the real DivType's": 220 DivTypeFlags = 31 221 }; 222 inline DivType operator|(DivType a, DivType b) 223 { 224 return static_cast<DivType>( 225 static_cast<int>(a) | static_cast<int>(b)); 226 } 227 inline DivType operator&(DivType a, DivType b) 228 { 229 return static_cast<DivType>( 230 static_cast<int>(a) & static_cast<int>(b)); 231 } 232 inline DivType& operator|=(DivType& a, DivType b) 233 { a = (a|b); return a; } 234 inline DivType& operator&=(DivType& a, DivType b) 235 { a = (a&b); return a; } 236 inline DivType operator~(DivType a) 237 { return static_cast<DivType>(~static_cast<int>(a)); } 238 239 #ifdef NOALIASCHECK 240 template <class M1, class M2> SameStorage(const M1 &,const M2 &)241 inline bool SameStorage(const M1& , const M2& ) 242 { return false; } 243 #else 244 template <class M1, class M2> SameStorage(const M1 & m1,const M2 & m2)245 inline bool SameStorage(const M1& m1, const M2& m2) 246 { 247 return 248 static_cast<const void*>(m1.realPart().cptr()) == 249 static_cast<const void*>(m2.realPart().cptr()); 250 } 251 #endif 252 253 template <typename T> TMV_SQR(T x)254 inline T TMV_SQR(T x) 255 { return x*x; } 256 257 template <typename T> TMV_SQRT(T x)258 inline T TMV_SQRT(T x) 259 { return T(std::sqrt(x)); } 260 261 template <typename T> TMV_EXP(T x)262 inline T TMV_EXP(T x) 263 { return T(std::exp(x)); } 264 265 template <typename T> TMV_LOG(T x)266 inline T TMV_LOG(T x) 267 { return T(std::log(x)); } 268 269 template <typename T> TMV_NORM(T x)270 inline T TMV_NORM(T x) 271 { return x*x; } 272 273 template <typename T> TMV_NORM(std::complex<T> x)274 inline T TMV_NORM(std::complex<T> x) 275 { return std::norm(x); } 276 277 template <typename T> TMV_CONJ(T x)278 inline T TMV_CONJ(T x) 279 { return x; } 280 281 template <typename T> TMV_CONJ(std::complex<T> x)282 inline std::complex<T> TMV_CONJ(std::complex<T> x) 283 { return std::conj(x); } 284 285 template <typename T> TMV_REAL(T x)286 inline T TMV_REAL(T x) 287 { return x; } 288 289 template <typename T> TMV_REAL(std::complex<T> x)290 inline T TMV_REAL(std::complex<T> x) 291 { return std::real(x); } 292 293 template <typename T> TMV_IMAG(T)294 inline T TMV_IMAG(T ) 295 { return T(0); } 296 297 template <typename T> TMV_IMAG(std::complex<T> x)298 inline T TMV_IMAG(std::complex<T> x) 299 { return std::imag(x); } 300 301 template <typename T> TMV_ARG(T x)302 inline T TMV_ARG(T x) 303 { return x >= T(0) ? T(1) : T(-1); } 304 305 template <typename T> TMV_ARG(std::complex<T> x)306 inline T TMV_ARG(std::complex<T> x) 307 { return arg(x); } 308 309 template <typename T> TMV_ABS(T x)310 inline T TMV_ABS(T x) 311 { return std::abs(x); } 312 313 template <typename T> TMV_ABS(std::complex<T> x)314 inline T TMV_ABS(std::complex<T> x) 315 { 316 // This is the same as the usual std::abs algorithm. 317 // However, I have come across implementations that don't 318 // protext against overflow, so I just duplicate it here. 319 320 T xr = x.real(); 321 T xi = x.imag(); 322 const T s = std::max(std::abs(xr),std::abs(xi)); 323 if (s == T(0)) return s; // Check for s == 0 324 xr /= s; 325 xi /= s; 326 return s * std::sqrt(xr*xr + xi*xi); 327 } 328 329 template <typename T> TMV_ABS2(T x)330 inline T TMV_ABS2(T x) 331 { return std::abs(x); } 332 333 template <typename T> TMV_ABS2(std::complex<T> x)334 inline T TMV_ABS2(std::complex<T> x) 335 { return std::abs(std::real(x)) + std::abs(std::imag(x)); } 336 337 template <typename T> TMV_MIN(T x,T y)338 inline T TMV_MIN(T x, T y) 339 { return x > y ? y : x; } 340 341 template <typename T> TMV_MAX(T x,T y)342 inline T TMV_MAX(T x, T y) 343 { return x > y ? x : y; } 344 345 template <typename T> TMV_SWAP(T & x,T & y)346 inline void TMV_SWAP(T& x, T& y) 347 { T z = x; x = y; y = z; } 348 349 template <typename T> 350 struct Traits 351 { 352 enum { isreal = true }; 353 enum { iscomplex = false }; 354 enum { isinteger = std::numeric_limits<T>::is_integer }; 355 enum { twoifcomplex = 1 }; 356 357 typedef T real_type; 358 typedef std::complex<T> complex_type; 359 }; 360 361 template <typename T> 362 struct Traits<std::complex<T> > 363 { 364 enum { isreal = false }; 365 enum { iscomplex = true }; 366 enum { isinteger = Traits<T>::isinteger }; 367 enum { twoifcomplex = 2 }; 368 369 typedef T real_type; 370 typedef std::complex<T> complex_type; 371 }; 372 373 template <typename T> 374 struct Traits<T&> : public Traits<T> {}; 375 template <typename T> 376 struct Traits<const T> : public Traits<T> {}; 377 template <typename T> 378 struct Traits<const T&> : public Traits<T> {}; 379 380 template <typename T1, typename T2> 381 struct Traits2 382 { 383 enum { sametype = false }; 384 enum { samebase = false }; 385 typedef T1 type; 386 }; 387 template <typename T1, typename T2> 388 struct Traits2<T1,std::complex<T2> > 389 { 390 enum { sametype = false }; 391 enum { samebase = false }; 392 typedef std::complex<typename Traits2<T1,T2>::type> type; 393 }; 394 template <typename T1, typename T2> 395 struct Traits2<std::complex<T1>,T2> 396 { 397 enum { sametype = false }; 398 enum { samebase = false }; 399 typedef std::complex<typename Traits2<T1,T2>::type> type; 400 }; 401 template <typename T1, typename T2> 402 struct Traits2<std::complex<T1>,std::complex<T2> > 403 { 404 enum { sametype = false }; 405 enum { samebase = false }; 406 typedef std::complex<typename Traits2<T1,T2>::type> type; 407 }; 408 template <typename T> 409 struct Traits2<T,T> 410 { 411 enum { sametype = true }; 412 enum { samebase = true }; 413 typedef T type; 414 }; 415 template <typename T> 416 struct Traits2<std::complex<T>,T> 417 { 418 enum { sametype = false }; 419 enum { samebase = true }; 420 typedef std::complex<T> type; 421 }; 422 template <typename T> 423 struct Traits2<T,std::complex<T> > 424 { 425 enum { sametype = false }; 426 enum { samebase = true }; 427 typedef std::complex<T> type; 428 }; 429 template <typename T> 430 struct Traits2<std::complex<T>,std::complex<T> > 431 { 432 enum { sametype = true }; 433 enum { samebase = true }; 434 typedef std::complex<T> type; 435 }; 436 // Specialize the real pairs for which the second value is the type to 437 // use for sums and products rather than the first. 438 template <> 439 struct Traits2<int,float> 440 { 441 enum { sametype = false }; 442 enum { samebase = false }; 443 typedef float type; 444 }; 445 template <> 446 struct Traits2<int,double> 447 { 448 enum { sametype = false }; 449 enum { samebase = false }; 450 typedef double type; 451 }; 452 template <> 453 struct Traits2<int,long double> 454 { 455 enum { sametype = false }; 456 enum { samebase = false }; 457 typedef long double type; 458 }; 459 template <> 460 struct Traits2<float,double> 461 { 462 enum { sametype = false }; 463 enum { samebase = false }; 464 typedef double type; 465 }; 466 template <> 467 struct Traits2<float,long double> 468 { 469 enum { sametype = false }; 470 enum { samebase = false }; 471 typedef long double type; 472 }; 473 template <> 474 struct Traits2<double,long double> 475 { 476 enum { sametype = false }; 477 enum { samebase = false }; 478 typedef long double type; 479 }; 480 481 482 #define TMV_RealType(T) typename tmv::Traits<T>::real_type 483 #define TMV_ComplexType(T) typename tmv::Traits<T>::complex_type 484 485 template <typename T> 486 inline bool isReal(T) 487 { return true; } 488 template <typename T> 489 inline bool isReal(std::complex<T>) 490 { return false; } 491 template <typename T> 492 inline bool isComplex(T x) 493 { return !isReal(x); } 494 template <typename T1, typename T2> 495 inline bool isSameType(T1,T2) 496 { return false; } 497 template <typename T> 498 inline bool isSameType(T,T) 499 { return true; } 500 501 #ifndef NOTHROW 502 class Error : public std::runtime_error 503 { 504 public : 505 inline Error(std::string _s1) throw() : 506 std::runtime_error("TMV Error: "+_s1) {} 507 virtual inline ~Error() throw() {} 508 // Typically, this is overridden by subclasses with more detail. 509 virtual inline void write(std::ostream& os) const throw() 510 { os << this->what() << std::endl; } 511 }; 512 513 inline std::ostream& operator<<(std::ostream& os, const Error& e) throw() 514 { e.write(os); return os; } 515 516 class FailedAssert : 517 public Error 518 { 519 public : 520 std::string failed_assert; 521 unsigned long line; 522 std::string file; 523 524 inline FailedAssert(std::string s, unsigned long l, 525 std::string f) throw() : 526 Error("Failed Assert statement "+s), 527 failed_assert(s), line(l), file(f) {} 528 virtual inline ~FailedAssert() throw() {} 529 virtual inline void write(std::ostream& os) const throw() 530 { 531 os<<"TMV Failed Assert: "<<failed_assert<<std::endl<< 532 "on line "<<line<<" in file "<<file<<std::endl; 533 } 534 }; 535 #endif 536 537 #ifndef TMV_DEBUG 538 #define TMVAssert(x) 539 #else 540 #ifdef NOTHROW 541 #define TMVAssert(x) do { if(!(x)) { \ 542 std::cerr<<"Failed Assert: "<<#x<<" on line "<<__LINE__<<" in file "<<__FILE__<<std::endl; exit(1); } } while (false) 543 #else 544 #define TMVAssert(x) do { if(!(x)) { \ 545 throw tmv::FailedAssert(#x,__LINE__,__FILE__); } } while(false) 546 #endif 547 #endif 548 #ifdef NOTHROW 549 #define TMVAssert2(x) do { if(!(x)) { \ 550 std::cerr<<"Failed Assert: "<<#x<<" on line "<<__LINE__<<" in file "<<__FILE__<<std::endl; exit(1); } } while (false) 551 #else 552 #define TMVAssert2(x) do { if(!(x)) { \ 553 throw tmv::FailedAssert(#x,__LINE__,__FILE__); } } while(false) 554 #endif 555 556 #ifdef __GNUC__ 557 #define TMV_DEPRECATED(func) func __attribute__ ((deprecated)) 558 #elif defined(_MSC_VER) 559 #define TMV_DEPRECATED(func) __declspec(deprecated) func 560 #elif defined(__PGI) 561 #define TMV_DEPRECATED(func) func __attribute__ ((deprecated)) 562 #else 563 #define TMV_DEPRECATED(func) func 564 #endif 565 566 // Use DEBUGPARAM(x) for parameters that are only used in TMVAssert 567 // statements. So then they don't give warnings when compiled with 568 // -DNDEBUG 569 #ifdef TMV_DEBUG 570 #define TMV_DEBUG_PARAM(x) x 571 #else 572 #define TMV_DEBUG_PARAM(x) 573 #endif 574 575 #ifndef NOTHROW 576 class ReadError : 577 public Error 578 { 579 public : 580 inline ReadError() throw() : 581 Error("Invalid istream input encountered.") {} 582 inline ReadError(std::string s) throw() : 583 Error("Invalid istream input encountered while reading "+s) {} 584 virtual inline ~ReadError() throw() {} 585 }; 586 587 // Each Vector, Matrix type has it's own derived ReadError class 588 // which outputs whatever has been read so far, and possibly 589 // what was expected instead of what was actually read in. 590 591 class Singular : 592 public Error 593 { 594 public : 595 inline Singular() throw() : 596 Error("Encountered singular matrix.") {} 597 inline Singular(std::string s) throw() : 598 Error("Encountered singular "+s) {} 599 virtual inline ~Singular() throw() {} 600 }; 601 602 class NonPosDef : 603 public Error 604 { 605 public: 606 inline NonPosDef() throw() : 607 Error("Invalid non-positive-definite matrix found.") {} 608 inline NonPosDef(std::string s) throw() : 609 Error("Non-positive-definite matrix found in "+s) {} 610 virtual inline ~NonPosDef() throw() {} 611 }; 612 #endif 613 614 #ifdef TMV_WARN 615 class TMV_WarnSingleton 616 { 617 // Note: This is not thread safe. 618 // If multiple threads write to the warning stream at the 619 // same time, they can clobber each other. 620 public: 621 static inline std::ostream*& inst() { 622 static std::ostream* warn = 0; 623 return warn; 624 } 625 private: 626 TMV_WarnSingleton(); 627 }; 628 629 inline void TMV_Warning(std::string s) 630 { 631 if (TMV_WarnSingleton::inst()) { 632 *TMV_WarnSingleton::inst() << "Warning:\n" << s << std::endl; 633 } 634 } 635 636 inline std::ostream* WriteWarningsTo(std::ostream* newos) 637 { 638 std::ostream* temp = TMV_WarnSingleton::inst(); 639 TMV_WarnSingleton::inst() = newos; 640 return temp; 641 } 642 643 inline void NoWarnings() 644 { TMV_WarnSingleton::inst() = 0; } 645 #else 646 inline void TMV_Warning(std::string ) {} 647 inline std::ostream* WriteWarningsTo(std::ostream* os) { return 0; } 648 inline void NoWarnings() {} 649 #endif 650 651 // Unknown is the value of _size, _step, etc. whenever it 652 // is not known at compile time. 653 // We use for this value the maximally negative int. 654 // In binary, this is a 1 followed by all zeros. 655 // Note: can't use numeric_limits<int>::min, since it is a function, 656 // not a compile-time constant. 657 // And we need Unknown as a compile-time constant. 658 const int Unknown = 1<<(sizeof(int)*8-1); 659 660 // A helper structure that acts like an int, 661 // but only bothers to make the integer if S == Unknown. 662 // It also checks the constructor if S != Unknown 663 template <int S> 664 struct CheckedInt 665 { 666 CheckedInt(ptrdiff_t TMV_DEBUG_PARAM(s)) { 667 #ifdef TMV_DEBUG 668 if (s != S) { 669 std::cerr<<"Mismatched CheckInt:\n"; 670 std::cerr<<"template parameter S = "<<S<<std::endl; 671 std::cerr<<"argument s = "<<s<<std::endl; 672 } 673 #endif 674 TMVAssert(s == S); 675 } 676 operator ptrdiff_t() const { return S; } 677 }; 678 template <> 679 struct CheckedInt<Unknown> 680 { 681 ptrdiff_t step; 682 CheckedInt(ptrdiff_t s) : step(s) {} 683 operator ptrdiff_t() const { return step; } 684 ~CheckedInt() {} 685 }; 686 687 template <typename T> 688 inline TMV_RealType(T) TMV_Epsilon() 689 { return std::numeric_limits<typename Traits<T>::real_type>::epsilon(); } 690 691 template <typename T> 692 inline bool TMV_Underflow(T x) 693 { 694 return tmv::Traits<T>::isinteger ? false : 695 TMV_ABS2(x) < std::numeric_limits<T>::min(); 696 } 697 698 template <typename T> 699 inline bool TMV_Underflow(std::complex<T> x) 700 { 701 return tmv::Traits<T>::isinteger ? false : 702 TMV_ABS2(x) < T(2) * std::numeric_limits<T>::min(); 703 } 704 705 template <typename T> 706 inline T TMV_Divide(T x, T y) 707 { return x / y; } 708 709 template <typename T> 710 inline std::complex<T> TMV_Divide(std::complex<T> x, T y) 711 { 712 // return x / y; 713 // Amazingly, some implementations get even this one wrong! 714 // So I need to do each component explicitly. 715 return std::complex<T>(x.real()/y,x.imag()/y); 716 } 717 718 template <typename T> 719 inline T TMV_InverseOf(T x) 720 { return T(1) / x; } 721 722 template <typename T> 723 inline std::complex<T> TMV_InverseOf(std::complex<T> x) 724 { 725 // The standard library's implemenation of complex division 726 // does not guard against overflow/underflow. So we have 727 // our own version here. 728 729 T xr = x.real(); 730 T xi = x.imag(); 731 if (std::abs(xr) > std::abs(xi)) { 732 xi /= xr; 733 T denom = xr*(T(1)+xi*xi); 734 return std::complex<T>(T(1)/denom,-xi/denom); 735 } else if (std::abs(xi) > T(0)) { 736 xr /= xi; 737 T denom = xi*(T(1)+xr*xr); 738 return std::complex<T>(xr/denom,T(-1)/denom); 739 } else { 740 // (xr,xi) = (0,0) 741 // division by zero. Just go ahead and do it. 742 return T(1) / xi; 743 } 744 } 745 746 747 template <typename T> 748 inline std::complex<T> TMV_Divide(T x, std::complex<T> y) 749 { return x * TMV_InverseOf(y); } 750 751 template <typename T> 752 inline std::complex<T> TMV_Divide(std::complex<T> x, std::complex<T> y) 753 { return x * TMV_InverseOf(y); } 754 755 template <typename T> 756 inline T TMV_SIGN(T x, T ) 757 { return x > 0 ? T(1) : T(-1); } 758 759 template <typename T> 760 inline std::complex<T> TMV_SIGN(std::complex<T> x, T absx) 761 { return absx > 0 ? TMV_Divide(x,absx) : std::complex<T>(1); } 762 763 764 extern bool TMV_FALSE; 765 // = false (in TMV_Vector.cpp), but without the unreachable returns 766 767 inline int TMV_ABS(const std::complex<int>& x) 768 { return int(TMV_ABS(std::complex<double>(std::real(x),std::imag(x)))); } 769 770 inline int TMV_ARG(const std::complex<int>& x) 771 { return int(TMV_ARG(std::complex<double>(std::real(x),std::imag(x)))); } 772 773 inline int TMV_SQRT(int x) 774 { return int(TMV_SQRT(double(x))); } 775 776 inline std::complex<int> TMV_SQRT(const std::complex<int>& x) 777 { 778 std::complex<double> temp = 779 TMV_SQRT(std::complex<double>(std::real(x),std::imag(x))); 780 return std::complex<int>(int(real(temp)),int(imag(temp))); 781 } 782 783 inline int TMV_EXP(int x) 784 { return int(TMV_EXP(double(x))); } 785 786 inline std::complex<int> TMV_EXP(const std::complex<int>& x) 787 { 788 std::complex<double> temp = 789 TMV_EXP(std::complex<double>(std::real(x),std::imag(x))); 790 return std::complex<int>(int(real(temp)),int(imag(temp))); 791 } 792 793 inline int TMV_LOG(int x) 794 { return int(TMV_LOG(double(x))); } 795 796 inline std::complex<int> TMV_LOG(const std::complex<int>& x) 797 { 798 std::complex<double> temp = 799 TMV_LOG(std::complex<double>(std::real(x),std::imag(x))); 800 return std::complex<int>(int(real(temp)),int(imag(temp))); 801 } 802 803 inline bool TMV_Underflow(int ) 804 { return false; } 805 806 template <typename T> 807 inline std::string TMV_Text(const T&) 808 { return std::string("Unknown (") + typeid(T).name() + ")"; } 809 810 inline std::string TMV_Text(const double&) 811 { return "double"; } 812 813 inline std::string TMV_Text(const float&) 814 { return "float"; } 815 816 inline std::string TMV_Text(const int&) 817 { return "int"; } 818 819 inline std::string TMV_Text(const long double&) 820 { return "long double"; } 821 822 template <typename T> 823 inline std::string TMV_Text(std::complex<T>) 824 { return std::string("complex<") + TMV_Text(T()) + ">"; } 825 826 inline std::string TMV_Text(ConjType c) 827 { return c == Conj ? "Conj" : c == NonConj ? "NonConj" : "VarConj"; } 828 829 inline std::string TMV_Text(IndexStyle i) 830 { return i == CStyle ? "CStyle" : "FortranStyle"; } 831 832 inline std::string TMV_Text(StepType s) 833 { return s == Unit ? "Unit" : "Step"; } 834 835 inline std::string TMV_Text(StorageType s) 836 { 837 return 838 s==ColMajor ? "ColMajor" : 839 s==RowMajor ? "RowMajor" : "DiagMajor"; 840 } 841 842 inline std::string TMV_Text(DiagType u) 843 { return u == UnitDiag ? "UnitDiag" : "NonUnitDiag"; } 844 845 inline std::string TMV_Text(UpLoType u) 846 { return u == Upper ? "Upper" : "Lower"; } 847 848 inline std::string TMV_Text(SymType s) 849 { return s == Sym ? "Sym" : "Herm"; } 850 851 inline std::string TMV_Text(DivType d) 852 { 853 return 854 d==XX ? "XX" : 855 d==LU ? "LU" : 856 d==CH ? "CH" : 857 d==QR ? "QR" : 858 d==QRP ? "QRP" : 859 d==SV ? "SV" : 860 "unkown DivType"; 861 } 862 863 //#define TMVFLDEBUG 864 865 #ifdef TMVFLDEBUG 866 #define TMV_FIRSTLAST ,_first,_last 867 #define TMV_PARAMFIRSTLAST(T) , const T* _first, const T* _last 868 #define TMV_DEFFIRSTLAST(f,l) ,_first(f),_last(l) 869 #define TMV_DEFFIRSTLAST2(f,l) :_first(f),_last(l) 870 #define TMV_FIRSTLAST1(f,l) ,f,l 871 #define TMV_SETFIRSTLAST(f,l) _first=(f); _last=(l); 872 #else 873 #define TMV_FIRSTLAST 874 #define TMV_PARAMFIRSTLAST(T) 875 #define TMV_DEFFIRSTLAST(f,l) 876 #define TMV_DEFFIRSTLAST2(f,l) 877 #define TMV_FIRSTLAST1(f,l) 878 #define TMV_SETFIRSTLAST(f,l) 879 #endif 880 881 #define TMV_ConjOf(T,C) ((isReal(T()) || C==Conj) ? NonConj : Conj) 882 883 // Use DEBUGPARAM(x) for parameters that are only used in TMVAssert 884 // statements. So then they don't give warnings when compiled with 885 // -DNDEBUG 886 #ifdef TMV_DEBUG 887 #define TMV_DEBUGPARAM(x) x 888 #else 889 #define TMV_DEBUGPARAM(x) 890 #endif 891 892 // Workaround for migration from auto_ptr -> unique_ptr. Use auto_ptr name, but all the 893 // functionality should work with either class. 894 #if __cplusplus >= 201103L 895 template <typename T> 896 using auto_ptr = std::unique_ptr<T>; 897 #else 898 using std::auto_ptr; 899 #endif 900 901 } // namespace tmv 902 903 #endif 904