268f476d8f35953d5f8f73b84ad56eae65702fcb
[charm.git] / src / util / ckcomplex.h
1 #ifndef __CKCOMPLEX_H__
2 #define __CKCOMPLEX_H__
3
4 #include "pup.h"
5
6 #if USE_FFTW_DECLS
7 #include "fftw.h"
8 typedef fftw_real  RealType;
9 #else
10 typedef double     RealType;
11 #endif
12
13 #include <cmath>
14
15 struct ckcomplex {
16     RealType  re;
17     RealType  im;   
18
19  
20     inline ckcomplex(RealType _re=0., RealType _im=0.): re(_re), im(_im){}
21     //    inline ckcomplex(RealType r) {re=r; im=0;}
22     //    inline ckcomplex(RealType r,RealType i) {re=r; im=i;}
23     
24     inline ~ckcomplex() {}
25
26     inline RealType getMagSqr(void) const { 
27       return re*re+im*im; 
28     }
29
30     inline ckcomplex operator+(ckcomplex a) { 
31       return ckcomplex(re+a.re,im+a.im); 
32     }
33
34     // note: not a member
35     inline friend ckcomplex operator-(ckcomplex lhs, ckcomplex rhs) {
36       return ckcomplex(lhs.re - rhs.re, lhs.im - rhs.im);
37     }
38
39     inline ckcomplex conj(void) { 
40         return ckcomplex(re, -im); 
41     }
42
43     inline void operator+=(ckcomplex a) { 
44       re+=a.re; im+=a.im; 
45     }
46     
47     // note: not a member
48     inline friend ckcomplex operator*(RealType lhs, ckcomplex rhs) {
49       return ckcomplex(rhs.re*lhs, rhs.im*lhs);
50     }
51
52     inline friend ckcomplex operator/(ckcomplex lhs, RealType rhs) {
53         return ckcomplex(lhs.re/rhs, lhs.im/rhs);
54     }
55
56     inline ckcomplex operator*(RealType a) { 
57       return ckcomplex(re*a, im*a); 
58     } 
59
60     inline bool notzero() const { return( (0.0 != re) ? true : (0.0 != im)); }
61
62     inline void operator*=(ckcomplex a) {        
63         RealType treal, tim;
64         treal = re * a.re - im * a.im;
65         tim = re * a.im + im * a.re;
66         re = treal;
67         im = tim;
68     }
69
70     inline ckcomplex operator*(ckcomplex a) {
71       return ckcomplex( re * a.re - im * a.im, re * a.im + im * a.re); 
72     }
73
74
75     inline void operator -= (ckcomplex a) {
76       re -= a.re; im -= a.im;
77     }
78
79
80     inline ckcomplex multiplyByi () {
81         return ckcomplex(-im, re);
82     }
83         
84     inline void * operator new[] (size_t size){
85         void *buf = malloc(size);
86         return buf;
87     }
88     
89     inline void operator delete[] (void *buf){
90         free(buf);
91     }
92 };
93
94 typedef ckcomplex complex;
95
96 PUPbytes(ckcomplex)
97
98
99 /// Overload std::isfinite for complex numbers. @note: Doesn't seem to be part of the standard
100 inline int isfinite(complex aNum) { return ( std::isfinite(aNum.re) && std::isfinite(aNum.im) ); }
101
102 /// Like std::norm, return the square of the distance from (0,0) in the complex number plane
103 inline RealType norm(complex aNum) { return ( aNum.re*aNum.re + aNum.im*aNum.im ); }
104 /// Return the distance from (0,0) in the complex plane
105 inline RealType abs(complex aNum) { return ( sqrt(norm(aNum)) ); }
106
107 #endif
108