API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Dec 22nd, 2014
169
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 8.32 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <stdlib.h>
  3. #include <stdio.h>
  4. #include <sys/time.h>
  5. #include <algorithm>
  6. #ifdef POWER
  7. #include <builtins.h>
  8. #endif
  9.  
  10. //----------------------------------------------------------------------
  11. const int N = 20000;
  12. const double L = 10.0;
  13. const int D = 3;
  14. const int X = 0;
  15. const int Y = 1;
  16. const int Z = 2;
  17. double q[N][D],p[N][D];
  18. const double dt = 0.001;
  19. const double C2 = 0.1;
  20. const double CUTOFF2 = L*L*0.9*0.9;
  21. //----------------------------------------------------------------------
  22. double
  23. myrand(void){
  24. return static_cast<double>(rand())/(static_cast<double>(RAND_MAX));
  25. }
  26. //----------------------------------------------------------------------
  27. double
  28. myclock(void){
  29. struct timeval t;
  30. gettimeofday(&t, NULL);
  31. return t.tv_sec + t.tv_usec*1e-6;
  32. }
  33. //----------------------------------------------------------------------
  34. void
  35. calcforce(void){
  36. for(int i=0;i<N-1;i++){
  37. for(int j=i+1;j<N;j++){
  38. double dx = q[j][X] - q[i][X];
  39. double dy = q[j][Y] - q[i][Y];
  40. double dz = q[j][Z] - q[i][Z];
  41. double r2 = (dx*dx + dy*dy + dz*dz);
  42. double r6 = r2*r2*r2;
  43. double df = ((24.0*r6-48.0)/(r6*r6*r2)+C2*8.0)*dt;
  44. p[i][X] += df*dx;
  45. p[i][Y] += df*dy;
  46. p[i][Z] += df*dz;
  47. p[j][X] -= df*dx;
  48. p[j][Y] -= df*dy;
  49. p[j][Z] -= df*dz;
  50. }
  51. }
  52. }
  53. //----------------------------------------------------------------------
  54. void
  55. calcforce_hand(void){
  56. for(int i=0;i<N-1;i++){
  57. const double qx = q[i][X];
  58. const double qy = q[i][Y];
  59. const double qz = q[i][Z];
  60. double px = p[i][X];
  61. double py = p[i][Y];
  62. double pz = p[i][Z];
  63. for(int j=i+1;j<N;j++){
  64. double dx = q[j][X] - qx;
  65. double dy = q[j][Y] - qy;
  66. double dz = q[j][Z] - qz;
  67. double r2 = (dx*dx + dy*dy + dz*dz);
  68. double r6 = r2*r2*r2;
  69. double df = ((24.0*r6-48.0)/(r6*r6*r2)+C2*8.0)*dt;
  70. px += df*dx;
  71. py += df*dy;
  72. pz += df*dz;
  73. p[j][X] -= df*dx;
  74. p[j][Y] -= df*dy;
  75. p[j][Z] -= df*dz;
  76. }
  77. p[i][X] = px;
  78. p[i][Y] = py;
  79. p[i][Z] = pz;
  80. }
  81. }
  82. //----------------------------------------------------------------------
  83. void
  84. calcforce_hand_if(void){
  85. for(int i=0;i<N-1;i++){
  86. const double qx = q[i][X];
  87. const double qy = q[i][Y];
  88. const double qz = q[i][Z];
  89. double px = p[i][X];
  90. double py = p[i][Y];
  91. double pz = p[i][Z];
  92. for(int j=i+1;j<N;j++){
  93. double dx = q[j][X] - qx;
  94. double dy = q[j][Y] - qy;
  95. double dz = q[j][Z] - qz;
  96. double r2 = (dx*dx + dy*dy + dz*dz);
  97. if(r2 > CUTOFF2) continue;
  98. double r6 = r2*r2*r2;
  99. double df = ((24.0*r6-48.0)/(r6*r6*r2)+C2*8.0)*dt;
  100. px += df*dx;
  101. py += df*dy;
  102. pz += df*dz;
  103. p[j][X] -= df*dx;
  104. p[j][Y] -= df*dy;
  105. p[j][Z] -= df*dz;
  106. }
  107. p[i][X] = px;
  108. p[i][Y] = py;
  109. p[i][Z] = pz;
  110. }
  111. }
  112. //----------------------------------------------------------------------
  113. double r2_flag[N];
  114. //----------------------------------------------------------------------
  115. void
  116. calcforce_hand_if_pre(void){
  117. for(int i=0;i<N-1;i++){
  118. const double qx = q[i][X];
  119. const double qy = q[i][Y];
  120. const double qz = q[i][Z];
  121. double px = p[i][X];
  122. double py = p[i][Y];
  123. double pz = p[i][Z];
  124. for(int j=i+1;j<N;j++){
  125. double dx = q[j][X] - qx;
  126. double dy = q[j][Y] - qy;
  127. double dz = q[j][Z] - qz;
  128. double r2 = (dx*dx + dy*dy + dz*dz);
  129. r2_flag[j] = (r2 > CUTOFF2)? 0.0: 1.0;
  130. }
  131. for(int j=i+1;j<N;j++){
  132. double dx = q[j][X] - qx;
  133. double dy = q[j][Y] - qy;
  134. double dz = q[j][Z] - qz;
  135. double r2 = (dx*dx + dy*dy + dz*dz);
  136. double r6 = r2*r2*r2;
  137. double df = ((24.0*r6-48.0)/(r6*r6*r2)+C2*8.0)*dt*r2_flag[j];
  138. px += df*dx;
  139. py += df*dy;
  140. pz += df*dz;
  141. p[j][X] -= df*dx;
  142. p[j][Y] -= df*dy;
  143. p[j][Z] -= df*dz;
  144. }
  145. p[i][X] = px;
  146. p[i][Y] = py;
  147. p[i][Z] = pz;
  148. }
  149. }
  150. //----------------------------------------------------------------------
  151. void
  152. calcforce_hand_if_min(void){
  153. for(int i=0;i<N-1;i++){
  154. const double qx = q[i][X];
  155. const double qy = q[i][Y];
  156. const double qz = q[i][Z];
  157. double px = p[i][X];
  158. double py = p[i][Y];
  159. double pz = p[i][Z];
  160. for(int j=i+1;j<N;j++){
  161. double dx = q[j][X] - qx;
  162. double dy = q[j][Y] - qy;
  163. double dz = q[j][Z] - qz;
  164. double r2 = (dx*dx + dy*dy + dz*dz);
  165. r2 = std::min(r2,CUTOFF2);
  166. double r6 = r2*r2*r2;
  167. double df = ((24.0*r6-48.0)/(r6*r6*r2)+C2*8.0)*dt;
  168. px += df*dx;
  169. py += df*dy;
  170. pz += df*dz;
  171. p[j][X] -= df*dx;
  172. p[j][Y] -= df*dy;
  173. p[j][Z] -= df*dz;
  174. }
  175. p[i][X] = px;
  176. p[i][Y] = py;
  177. p[i][Z] = pz;
  178. }
  179. }
  180.  
  181. //----------------------------------------------------------------------
  182. void
  183. calcforce_hand_if_condop(void){
  184. for(int i=0;i<N-1;i++){
  185. const double qx = q[i][X];
  186. const double qy = q[i][Y];
  187. const double qz = q[i][Z];
  188. double px = p[i][X];
  189. double py = p[i][Y];
  190. double pz = p[i][Z];
  191. for(int j=i+1;j<N;j++){
  192. double dx = q[j][X] - qx;
  193. double dy = q[j][Y] - qy;
  194. double dz = q[j][Z] - qz;
  195. double r2 = (dx*dx + dy*dy + dz*dz);
  196. r2 = (r2<CUTOFF2)? r2:CUTOFF2;
  197. double r6 = r2*r2*r2;
  198. double df = ((24.0*r6-48.0)/(r6*r6*r2)+C2*8.0)*dt;
  199. px += df*dx;
  200. py += df*dy;
  201. pz += df*dz;
  202. p[j][X] -= df*dx;
  203. p[j][Y] -= df*dy;
  204. p[j][Z] -= df*dz;
  205. }
  206. p[i][X] = px;
  207. p[i][Y] = py;
  208. p[i][Z] = pz;
  209. }
  210. }
  211.  
  212. //----------------------------------------------------------------------
  213. void
  214. calcforce_hand_if2(void){
  215. for(int i=0;i<N-1;i++){
  216. const double qx = q[i][X];
  217. const double qy = q[i][Y];
  218. const double qz = q[i][Z];
  219. double px = p[i][X];
  220. double py = p[i][Y];
  221. double pz = p[i][Z];
  222. for(int j=i+1;j<N;j++){
  223. double dx = q[j][X] - qx;
  224. double dy = q[j][Y] - qy;
  225. double dz = q[j][Z] - qz;
  226. double r2 = (dx*dx + dy*dy + dz*dz);
  227. double r6 = r2*r2*r2;
  228. double df = ((24.0*r6-48.0)/(r6*r6*r2)+C2*8.0)*dt;
  229. #ifdef POWER
  230. df = __fsel(r2 - CUTOFF2,df,0.0);
  231. #else
  232. df = (r2 > CUTOFF2)? 0.0: df;
  233. #endif
  234. px += df*dx;
  235. py += df*dy;
  236. pz += df*dz;
  237. p[j][X] -= df*dx;
  238. p[j][Y] -= df*dy;
  239. p[j][Z] -= df*dz;
  240. }
  241. p[i][X] = px;
  242. p[i][Y] = py;
  243. p[i][Z] = pz;
  244. }
  245. }
  246. //----------------------------------------------------------------------
  247. void
  248. calcforce_if(void){
  249. for(int i=0;i<N-1;i++){
  250. for(int j=i+1;j<N;j++){
  251. double dx = q[j][X] - q[i][X];
  252. double dy = q[j][Y] - q[i][Y];
  253. double dz = q[j][Z] - q[i][Z];
  254. double r2 = (dx*dx + dy*dy + dz*dz);
  255. if(r2 > CUTOFF2) continue;
  256. double r6 = r2*r2*r2;
  257. double df = ((24.0*r6-48.0)/(r6*r6*r2)+C2*8.0)*dt;
  258. p[i][X] += df*dx;
  259. p[i][Y] += df*dy;
  260. p[i][Z] += df*dz;
  261. p[j][X] -= df*dx;
  262. p[j][Y] -= df*dy;
  263. p[j][Z] -= df*dz;
  264. }
  265. }
  266. }
  267. //----------------------------------------------------------------------
  268. void
  269. calcforce_if2(void){
  270. for(int i=0;i<N-1;i++){
  271. for(int j=i+1;j<N;j++){
  272. double dx = q[j][X] - q[i][X];
  273. double dy = q[j][Y] - q[i][Y];
  274. double dz = q[j][Z] - q[i][Z];
  275. double r2 = (dx*dx + dy*dy + dz*dz);
  276. double r6 = r2*r2*r2;
  277. double df = ((24.0*r6-48.0)/(r6*r6*r2)+C2*8.0)*dt;
  278. #ifdef POWER
  279. df = __fsel(r2 - CUTOFF2,df,0.0);
  280. #else
  281. df = (r2 > CUTOFF2)? 0.0: df;
  282. #endif
  283. p[i][X] += df*dx;
  284. p[i][Y] += df*dy;
  285. p[i][Z] += df*dz;
  286. p[j][X] -= df*dx;
  287. p[j][Y] -= df*dy;
  288. p[j][Z] -= df*dz;
  289. }
  290. }
  291. }
  292.  
  293. //----------------------------------------------------------------------
  294. void
  295. measure(void(*pfunc)(),const char *name){
  296. double st = myclock();
  297. pfunc();
  298. double t = myclock()-st;
  299. printf("N=%d, %s %f [sec]\n",N,name,t);
  300. }
  301. //----------------------------------------------------------------------
  302. int
  303. main(void){
  304. for(int i=0;i<N;i++){
  305. p[i][X] = 0.0;
  306. p[i][Y] = 0.0;
  307. p[i][Z] = 0.0;
  308. q[i][X] = L*myrand();
  309. q[i][Y] = L*myrand();
  310. q[i][Z] = L*myrand();
  311. }
  312. measure(&calcforce,"calcforce");
  313. measure(&calcforce_hand,"calcforce_hand");
  314. measure(&calcforce_if,"calcforce_if");
  315. measure(&calcforce_hand_if,"calcforce_hand_if");
  316. measure(&calcforce_hand_if_min,"calcforce_hand_if_min");
  317. measure(&calcforce_hand_if_min,"calcforce_hand_if_condop");
  318. }
  319. //----------------------------------------------------------------------
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!