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