NAMD
Vector.h
Go to the documentation of this file.
1 
7 #ifndef VECTOR_H
8 #define VECTOR_H
9 
10 #include <math.h>
11 #include <stdio.h>
12 #include "common.h"
13 
14 class Vector;
15 
16 class FloatVector {
17  public:
18  float x,y,z;
19  inline FloatVector(void) { ; }
20  inline FloatVector(const Vector &v);
21 };
22 
23 #ifdef ARCH_POWERPC
24 
25 #include "builtins.h"
26 
27 inline double namd_rsqrt(double x)
28 {
29  double r0, r1, r2;
30 
31  /*------------------------------------------*/
32  /* use reciprocal sqrt estimate instruction */
33  /*------------------------------------------*/
34  r0 = __frsqrte(x);
35 
36  /*----------------------*/
37  /* 1st newton iteration */
38  /*----------------------*/
39  r1 = r0 + 0.5*r0*(1.0 - x*r0*r0);
40 
41  /*----------------------*/
42  /* 2nd newton iteration */
43  /*----------------------*/
44  r2 = r1 + 0.5*r1*(1.0 - x*r1*r1);
45 
46  return r2;
47 }
48 
49 inline double namd_reciprocal (double x) {
50  register double rx;
51 
52  rx = __fres(x); // 0th estimate (13 bits)
53  rx = rx + rx*(1.0 - x*rx); // 1st Newton iteration (26 bits)
54  rx = rx + rx*(1.0 - x*rx); // 2nd Newton iteration (52 bits = full mantissa for a double)
55 
56  return rx;
57 }
58 
59 #else
60 #define namd_rsqrt(x) (1.0 / sqrt (x))
61 #define namd_reciprocal(x) (1.0 / x)
62 #endif
63 
64 class Vector {
65  public:
67 
68  #ifndef NAMD_AVXTILES
69  inline Vector(void) : x(-99999), y(-99999), z(-99999) { ; }
70  #else
71  inline Vector(void) { ; }
72  #endif
73 
74  inline Vector( BigReal newx, BigReal newy, BigReal newz)
75  : x(newx), y(newy), z(newz) { ; }
76 
77  inline Vector( BigReal newv ) // allow Vector v = 0; etc.
78  : x(newv), y(newv), z(newv) { ; }
79 
80  inline Vector(const FloatVector &v) : x(v.x), y(v.y), z(v.z) { ; }
81 
82  inline BigReal &operator[](int i) {
83  return i==0 ? x
84  :i==1 ? y
85  :i==2 ? z
86  :(NAMD_die("Vector reference out of bounds."), x);
87 
88  }
89 
90  // v1 = const;
91  inline Vector& operator=(const BigReal &v2) {
92  x = v2;
93  y = v2;
94  z = v2;
95  return *this;
96  }
97 
98  // v1 += v2;
99  inline void operator+=(const Vector &v2) {
100  x += v2.x;
101  y += v2.y;
102  z += v2.z;
103  }
104 
105  // v1 -= v2;
106  inline void operator-=(const Vector &v2) {
107  x -= v2.x;
108  y -= v2.y;
109  z -= v2.z;
110  }
111 
112  // v1 *= const
113  inline void operator*=(const BigReal &v2) {
114  x *= v2;
115  y *= v2;
116  z *= v2;
117  }
118 
119  // v1 /= const
120  inline void operator/=(const BigReal& v2) {
121  BigReal v2_recip = namd_reciprocal(v2);
122  x *= v2_recip;
123  y *= v2_recip;
124  z *= v2_recip;
125  }
126 
127  inline friend int operator == (const Vector& v1, const Vector& v2) {
128  return v1.x == v2.x && v1.y == v2.y && v1.z == v2.z;
129  }
130  inline friend int operator != (const Vector& v1, const Vector& v2) {
131  // return !(v1.x == v2.x && v1.y == v2.y && v1.z == v2.z);
132  return v1.x != v2.x || v1.y != v2.y || v1.z != v2.z;
133  }
134 
135  // addition of two vectors
136  inline friend Vector operator+(const Vector& v1, const Vector& v2) {
137  return Vector( v1.x+v2.x, v1.y+v2.y, v1.z+v2.z);
138  }
139 
140  // negation
141  inline friend Vector operator-(const Vector &v1) {
142  return Vector( -v1.x, -v1.y, -v1.z);
143  }
144 
145  // subtraction
146  inline friend Vector operator-(const Vector &v1, const Vector &v2) {
147  return Vector( v1.x-v2.x, v1.y-v2.y, v1.z-v2.z);
148  }
149  // inner ("dot") product
150  inline friend BigReal operator*(const Vector &v1, const Vector &v2) {
151  return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
152  }
153  // scalar product
154  inline friend Vector operator*(const BigReal &f, const Vector &v1) {
155  return Vector(f*v1.x, f*v1.y, f*v1.z);
156  }
157  // scalar product
158  inline friend Vector operator*(const Vector &v1, const BigReal &f) {
159  return Vector(f*v1.x, f*v1.y, f*v1.z);
160  }
161  // division by a scalar
162  inline friend Vector operator/(const Vector &v1, const BigReal &f) {
163 // if (!f)
164 // NAMD_die("Division by 0 on a vector operation.");
165  return Vector(v1.x/f, v1.y/f, v1.z/f);
166  }
167 
168  // return the norm
169  inline BigReal length(void) const {
170  return sqrt(x*x+y*y+z*z);
171  }
172 
173  inline BigReal length2(void) const {
174  return (x*x + y*y + z*z);
175  }
176 
177  inline BigReal rlength (void) {
178  return namd_rsqrt (x*x + y*y + z*z);
179  }
180 
181  // return the unit vector in the same direction
182  inline Vector unit(void) const {
183  return Vector(x, y, z)/length();
184  }
185 
186 
187  // one cross product v3 = cross(v1, v2)
188  inline friend Vector cross(const Vector &v1, const Vector &v2) {
189  return Vector( v1.y*v2.z-v2.y*v1.z,
190  // -v1.x*v2.z+v2.x*v1.z,
191  v2.x*v1.z-v1.x*v2.z,
192  v1.x*v2.y-v2.x*v1.y );
193  }
194 
195  // multiplying a cross product by a scalar is very common
196  // one cross product v3 = k*cross(v1, v2)
197  inline friend Vector cross(const Real &k, const Vector &v1, const Vector &v2) {
198  return Vector( k*(v1.y*v2.z-v2.y*v1.z),
199  // k*(-v1.x*v2.z+v2.x*v1.z),
200  k*(v2.x*v1.z-v1.x*v2.z),
201  k*(v1.x*v2.y-v2.x*v1.y) );
202  }
203 
204  inline friend Vector cross(const BigReal &k, const Vector &v1, const Vector &v2) {
205  return Vector( k*(v1.y*v2.z-v2.y*v1.z),
206  // k*(-v1.x*v2.z+v2.x*v1.z),
207  k*(v2.x*v1.z-v1.x*v2.z),
208  k*(v1.x*v2.y-v2.x*v1.y) );
209  }
210 
211  // A = A x B -- why do you want this function, anyway?
212  void cross(const Vector &v2) {
213  BigReal xx = y*v2.z-v2.y*z;
214  // BigReal yy = -x*v2.z+v2.x*z;
215  BigReal yy = v2.x*z-x*v2.z;
216  z = x*v2.y-v2.x*y;
217  y=yy;
218  x=xx;
219  }
220 
221  // returns (*this) * V2
222  BigReal dot(const Vector &v2) {
223  return x*v2.x + y*v2.y + z*v2.z;
224  }
225 
226  // set the vector based on a string. If bad, return FALSE
227  // the string can be in the form "x y z" or "x, y, z"
228  Bool set(const char *s) {
229  double a[3]; // read into doubles, since I don't know what
230  char tmp[100]; // a "BigReal" is in real life
231  // cheap way to get commas, etc. a poor regex
232  int i=sscanf(s, "%lf%99[ \t,]%lf%99[ \t,]%lf%99s",
233  a, tmp, a+1, tmp, a+2, tmp);
234  if (i != 5) return FALSE;
235  const char *t = s; // now count commas (for "1,,,,2, , 3")
236  int flg = 0; // and check for "1 2,,3"
237  i = 0;
238  for (;*t;t++) {
239  if (*t == ',') {
240  if (flg == 0) { // expecting non-whitespace
241  return FALSE; // so error
242  }
243  flg = 0; // now expect non-whitespace
244  i++; // and increment comma counter
245  }
246  else if (*t != ' ' && *t != '\t') { // got non-whitespace
247  flg = 1; // so the next can be whitespace or commas
248  }
249  }
250  if (i == 0 || i == 2) { // allow "1 2 3" or "1, 2,3" forms
251  x = a[0]; y = a[1]; z = a[2];
252  return TRUE;
253  }
254  return FALSE;
255  }
256 };
257 
258 class zVector : public Vector {
259  public:
260  inline zVector(void) : Vector(0,0,0) { ; }
261  inline zVector(const Vector &v) : Vector(v) { ; }
262 };
263 
264 
265 class AlignVector : public Vector {
266  public:
268  inline AlignVector(void) : Vector(0,0,0) { pad = 0.0; }
269  inline AlignVector(const Vector &v) : Vector(v) { pad = 0.0; }
270 
271  inline AlignVector( BigReal newx, BigReal newy, BigReal newz)
272  : Vector (newx, newy, newz) { pad = 0.0; }
273 
274  inline AlignVector( BigReal newv ) // allow Vector v = 0; etc.
275  : Vector (newv) { pad = 0.0; }
276 
277  inline AlignVector(const FloatVector &v) : Vector (v) { pad = 0.0; }
278 };
279 
280 
281 inline FloatVector::FloatVector(const Vector &v) : x(v.x), y(v.y), z(v.z) { ; }
282 
283 //#define TEST_VECTOR_CLASS
284 #ifdef TEST_VECTOR_CLASS
285 main()
286 {
287  Vector v1(1.1,2.2, 3.3);
288  Vector v2(-1, 55, 32.1);
289  Vector v3(v1+2*v2);
290  Vector v4;
291  std::cout << v1 << " " << v2 << " " << v3 << " " << v4 << '\n';
292  std::cout << v1*v2 << " " << v3-v1-2*v2 <<" "<< v2 * v3 <<" "<< v3*v2 <<'\n';
293  v4 = v3*5 - v2/4;
294  std::cout << v4 << " " << v3*5.0 - v2/4.0 << '\n';
295  std::cout << v4[0] << " " << v4[1] << " " << v4[2] << '\n';
296 // std::cout.flush();
297 // std::cout << v4[3];
298  std::cout << cross(v1, v2) << '\n';
299  std::cout << v1 << '\n';
300  v1 += v2;
301  std::cout << v1 << '\n';
302  v1 -= v2;
303  std::cout << v1 << '\n';
304  {
305  Vector v1(1.0, 2.0, 3.0); // some more examples, but I was too lazy to
306  Vector v2 = v1.unit(); // fix the names
307  std::cout << v2 << '\n';
308  std::cout << v2.dot(v1) << '\n';
309  std::cout << v1.length() << '\n';
310  v1 *= -1;
311  v1 += v2*14;
312  std::cout << v1 << '\n';
313  }
314 }
315 #endif
316 
317 #endif
318 
void operator/=(const BigReal &v2)
Definition: Vector.h:120
Vector & operator=(const BigReal &v2)
Definition: Vector.h:91
friend Vector cross(const Real &k, const Vector &v1, const Vector &v2)
Definition: Vector.h:197
friend int operator==(const Vector &v1, const Vector &v2)
Definition: Vector.h:127
Definition: Vector.h:64
friend BigReal operator*(const Vector &v1, const Vector &v2)
Definition: Vector.h:150
AlignVector(const FloatVector &v)
Definition: Vector.h:277
float Real
Definition: common.h:109
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
BigReal z
Definition: Vector.h:66
float y
Definition: Vector.h:18
#define FALSE
Definition: common.h:118
friend int operator!=(const Vector &v1, const Vector &v2)
Definition: Vector.h:130
AlignVector(BigReal newx, BigReal newy, BigReal newz)
Definition: Vector.h:271
AlignVector(const Vector &v)
Definition: Vector.h:269
BigReal length(void) const
Definition: Vector.h:169
friend Vector cross(const Vector &v1, const Vector &v2)
Definition: Vector.h:188
void operator-=(const Vector &v2)
Definition: Vector.h:106
AlignVector(void)
Definition: Vector.h:268
#define namd_rsqrt(x)
Definition: Vector.h:60
friend Vector operator/(const Vector &v1, const BigReal &f)
Definition: Vector.h:162
BigReal rlength(void)
Definition: Vector.h:177
BigReal dot(const Vector &v2)
Definition: Vector.h:222
float x
Definition: Vector.h:18
float z
Definition: Vector.h:18
Vector(const FloatVector &v)
Definition: Vector.h:80
BigReal pad
Definition: Vector.h:267
friend Vector operator*(const Vector &v1, const BigReal &f)
Definition: Vector.h:158
gridSize z
int Bool
Definition: common.h:133
Vector(BigReal newv)
Definition: Vector.h:77
friend Vector operator-(const Vector &v1, const Vector &v2)
Definition: Vector.h:146
virial xx
BigReal x
Definition: Vector.h:66
void NAMD_die(const char *err_msg)
Definition: common.C:85
friend Vector operator-(const Vector &v1)
Definition: Vector.h:141
friend Vector operator*(const BigReal &f, const Vector &v1)
Definition: Vector.h:154
virial yy
zVector(const Vector &v)
Definition: Vector.h:261
Bool set(const char *s)
Definition: Vector.h:228
Vector(void)
Definition: Vector.h:69
BigReal length2(void) const
Definition: Vector.h:173
BigReal & operator[](int i)
Definition: Vector.h:82
FloatVector(void)
Definition: Vector.h:19
void cross(const Vector &v2)
Definition: Vector.h:212
BigReal y
Definition: Vector.h:66
int main(int argc, char *argv[])
Definition: diffbinpdb.c:15
friend Vector operator+(const Vector &v1, const Vector &v2)
Definition: Vector.h:136
zVector(void)
Definition: Vector.h:260
void operator+=(const Vector &v2)
Definition: Vector.h:99
gridSize y
Vector(BigReal newx, BigReal newy, BigReal newz)
Definition: Vector.h:74
friend Vector cross(const BigReal &k, const Vector &v1, const Vector &v2)
Definition: Vector.h:204
gridSize x
#define namd_reciprocal(x)
Definition: Vector.h:61
#define TRUE
Definition: common.h:119
void operator*=(const BigReal &v2)
Definition: Vector.h:113
Vector unit(void) const
Definition: Vector.h:182
AlignVector(BigReal newv)
Definition: Vector.h:274
double BigReal
Definition: common.h:114