API tools faq
paste
Advertisement
Guest User

Untitled

a guest
May 24th, 2019
81
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 19.35 KB | None | 0 0
raw download clone embed print report
  1. program Helmholtz
  2. ! PARA
  3.  
  4. use mpi
  5. ! PARA
  6.  
  7. implicit none
  8.  
  9. ! PARA
  10.  
  11. integer, parameter :: master_id = 0
  12.  
  13. integer :: id, tag, numprocs
  14. integer :: mpi_error_code
  15. integer, dimension(MPI_STATUS_SIZE) :: status
  16.  
  17.  
  18. ! PARA
  19.  
  20. real(kind=8), parameter :: ystar(2) = (0,0), kappa=2, eta=kappa, &
  21. eulerconst = 0.5772156649, pi = 4.D0*DATAN(1.D0)
  22. integer :: n = 1, k, j
  23. real(kind=8) :: testpoint(2) = (3,3)
  24. real(kind=8) , allocatable :: matrix(:,:), rightpart (:), density(:)
  25.  
  26. ! PARA
  27.  
  28. call MPI_init(mpi_error_code)
  29. call MPI_comm_rank(MPI_COMM_WORLD, id, mpi_error_code)
  30. call MPI_comm_size(MPI_COMM_WORLD, numprocs, mpi_error_code)
  31.  
  32.  
  33. !if (id == 0) then
  34. ! PARA
  35.  
  36. do while (n<=20)
  37. ALLOCATE (matrix(4*n,4*n), rightpart(4*n), density(4*n))
  38.  
  39. if (id /= master_id) then
  40. do k = id - 1, 2*n-1, numprocs - 1 !myid, 2*n-1, numprocs
  41. do j = 0, 2*n-1
  42. Matrix(k + 1, j+1) = Lre(k, j);
  43. Matrix(k + 1, j + 2 * n+1) = -Lim(k, j);
  44. Matrix(k + 2 * n+1, j+1) = Lim(k, j);
  45. Matrix(k + 2 * n+1, j + 2 * n+1) = Lre(k, j);
  46. if (k == j) then
  47. Matrix(k+1, j + 2 * n+1) = Matrix(k+1, j + 2 * n+1) - aim(t(k));
  48. Matrix(k + 2 * n+1, j+1) = Matrix(k + 2 * n+1, j+1) + aim(t(k));
  49. endif
  50. end do
  51. tag = 0
  52. call MPI_send(Matrix(k + 1, :), 4*n, MPI_DOUBLE_PRECISION, &
  53. master_id, tag, &
  54. MPI_COMM_WORLD, mpi_error_code)
  55. tag = 1
  56. call MPI_send(Matrix(k + 2 * n+1, :), 4*n, MPI_DOUBLE_PRECISION, &
  57. master_id, tag, &
  58. MPI_COMM_WORLD, mpi_error_code)
  59. end do
  60. else
  61. do k = 0, 2*n-1 !myid, 2*n-1, numprocs
  62. tag = 0
  63. call MPI_recv(Matrix(k + 1, :), 4*n, MPI_DOUBLE_PRECISION, &
  64. mod(k, numprocs - 1) + 1, tag, &
  65. MPI_COMM_WORLD, status, mpi_error_code)
  66. tag = 1
  67. call MPI_recv(Matrix(k + 2 * n+1, :), 4*n, MPI_DOUBLE_PRECISION, &
  68. mod(k, numprocs - 1) + 1, tag, &
  69. MPI_COMM_WORLD, status, mpi_error_code)
  70. rightPart(k+1) = fre(t(k));
  71. rightPart(k + 2 * n+1) = fim(t(k));
  72. end do
  73. end if
  74. density = solve(Matrix,rightPart)
  75.  
  76.  
  77. if (id == master_id) then
  78. write(*, *) Matrix
  79. write(*, *) rightPart
  80. write(*, *) sum(Matrix(:, 1) * rightPart)
  81. write(*, *) density
  82. write(*,*) "n = ", n, " ", &
  83. Abs(resultRe(testPoint, density) - exactRe(testPoint)), " ", &
  84. Abs(resultIm(testPoint, density) - exactIm(testPoint))
  85. end if
  86. DEALLOCATE (matrix, rightpart, density)
  87. n = n*2
  88. end do
  89.  
  90. ! PARA
  91. !endif
  92. call MPI_FINALIZE(mpi_error_code)
  93.  
  94. ! PARA
  95.  
  96. contains
  97.  
  98. real(kind=8) function x1(t)
  99. real(kind=8) t
  100. x1 = 2.*cos(t)
  101. return
  102. end
  103.  
  104. real(kind=8) function x2(t)
  105. real(kind=8) t
  106. x2 = 2.*sin(t)
  107. return
  108. end
  109.  
  110. real(kind=8) function x1d(t)
  111. real(kind=8) t
  112. x1d = -2.*sin(t)
  113. return
  114. end
  115.  
  116. real(kind=8) function x2d(t)
  117. real(kind=8) t
  118. x2d = 2.*cos(t)
  119. return
  120. end
  121.  
  122. real(kind=8) function x1dd(t)
  123. real(kind=8) t
  124. x1dd = -2.*cos(t)
  125. return
  126. end
  127.  
  128. real(kind=8) function x2dd(t)
  129. real(kind=8) t
  130. x2dd = -2.*sin(t)
  131. return
  132. end
  133.  
  134. real(kind=8) function x1ddd(t)
  135. real(kind=8) t
  136. x1ddd = 2.*sin(t)
  137. return
  138. end
  139.  
  140. real(kind=8) function x2ddd(t)
  141. real(kind=8) t
  142. x2ddd = -2.*cos(t)
  143. return
  144. end
  145.  
  146. real(kind=8) function diffmodule(t, tau)
  147. real(kind=8) t, tau
  148. diffmodule = SQRT((x1(t) - x1(tau)) * (x1(t) - x1(tau)) + (x2(t) - x2(tau)) * (x2(t) - x2(tau)))
  149. return
  150. end
  151.  
  152. real(kind=8) function dermodule(t)
  153. real(kind=8) t
  154. dermodule = SQRT(x1d(t) * x1d(t) + x2d(t) * x2d(t))
  155. return
  156. end
  157.  
  158. real(kind=8) function doubledermodule(t)
  159. real(kind=8) t
  160. doubledermodule = SQRT(x1dd(t) * x1dd(t) + x2dd(t) * x2dd(t))
  161. return
  162. end
  163.  
  164. real(kind=8) function h(t, tau)
  165. real(kind=8) t, tau, numerator, denominator
  166. numerator = ((x1(tau) - x1(t)) * x2d(t) - (x2(tau) - x2(t)) * x1d(t))
  167. denominator = diffModule(t, tau)
  168. h = numerator / denominator
  169. return
  170. end
  171.  
  172. real(kind=8) function him(t, tau)
  173. real(kind=8) t, tau
  174. him = kappa * h(t, tau) * BESSEL_J1(kappa * diffModule(t, tau)) * derModule(tau) / 2.0
  175. return
  176. end
  177.  
  178. real(kind=8) function hre(t, tau)
  179. real(kind=8) t, tau
  180. hre = -kappa * h(t, tau) * BESSEL_Y1(kappa * diffModule(t, tau)) * derModule(tau) / 2.0;
  181. return
  182. end
  183.  
  184. real(kind=8) function h1(t, tau)
  185. real(kind=8) t, tau
  186. if (t /= tau) then
  187. h1 = -kappa * h(t, tau) * BESSEL_J1(kappa * diffModule(t, tau)) * derModule(tau) / (2.0 * pi)
  188. return
  189. else
  190. h1 = 0.0
  191. return
  192. end if
  193. end
  194.  
  195. real(kind=8) function h2re(t, tau)
  196. real(kind=8) t, tau
  197. if (t /= tau) then
  198. h2re = hre(t, tau) - h1(t, tau) * log(4.0 * sin((t - tau) / 2.0) * sin((t - tau) / 2.0))
  199. return
  200. else
  201. h2re = (x2d(t) * x1dd(t) - x1d(t) * x2dd(t)) / (2.0 * pi * derModule(t))
  202. return
  203. end if
  204. end
  205.  
  206. real(kind=8) function h2im(t, tau)
  207. real(kind=8) t, tau
  208. if (t /= tau) then
  209. h2im = him(t, tau)
  210. return
  211. else
  212. h2im = 0.0
  213. return
  214. end if
  215. end
  216.  
  217. real(kind=8) function mim(t, tau)
  218. real(kind=8) t, tau
  219. mim = BESSEL_J0(kappa * diffModule(t, tau)) / 2.0
  220. return
  221. end
  222.  
  223. real(kind=8) function mre(t, tau)
  224. real(kind=8) t, tau
  225. mre = -BESSEL_Y0(kappa * diffModule(t, tau)) / 2.0
  226. return
  227. end
  228.  
  229. real(kind=8) function m1(t, tau)
  230. real(kind=8) t, tau
  231. if (t /= tau) then
  232. m1 = -BESSEL_J0(kappa * diffModule(t, tau)) / (2.0 * pi)
  233. return
  234. else
  235. m1 = -1 / (2.0 * pi)
  236. return
  237. end if
  238. end
  239.  
  240. real(kind=8) function m2re(t, tau)
  241. real(kind=8) t, tau
  242. if (t /= tau) then
  243. m2re = mre(t, tau) - m1(t, tau) * log(4.0 * sin((t - tau) / 2.0) * sin((t - tau) / 2.0))
  244. return
  245. else
  246. m2re = -(eulerconst + log(kappa * derModule(t) / 2.0)) / pi
  247. return
  248. end if
  249. end
  250.  
  251. real(kind=8) function m2im(t, tau)
  252. real(kind=8) t, tau
  253. if (t /= tau) then
  254. m2im = mim(t, tau)
  255. return
  256. else
  257. m2im = 1.0 / 2.0
  258. return
  259. end if
  260. end
  261.  
  262. real(kind=8) function tildan(t, tau)
  263. real(kind=8) t, tau, numerator, denominator
  264. numerator = ((x1(t) - x1(tau)) * x1d(t) + (x2(t) - x2(tau)) * x2d(t)) * &
  265. ((x1(t) - x1(tau)) * x1d(tau) + (x2(t) - x2(tau)) * x2d(tau))
  266. denominator = diffModule(t, tau) * diffModule(t, tau)
  267. tildan = numerator / denominator
  268. return
  269. end
  270.  
  271. real(kind=8) function nre(t, tau)
  272. real(kind=8) t, tau
  273. nre = 1.0 / (4.0 * pi * sin((t - tau) / 2.0) * sin((t - tau) / 2.0)) + tildan(t, tau) * &
  274. (-kappa * kappa * BESSEL_Y0(kappa * diffModule(t, tau)) + &
  275. 2.0 * kappa * BESSEL_Y1(kappa * diffModule(t, tau)) / diffModule(t, tau)) / 2.0 &
  276. - kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * BESSEL_Y1(kappa * diffModule(t, tau)) / (2.0 * diffModule(t, tau))
  277. return
  278. end
  279.  
  280. real(kind=8) function nim(t, tau)
  281. real(kind=8) t, tau
  282. nim = tildan(t, tau) * (kappa * kappa * BESSEL_J0(kappa * diffModule(t, tau)) - &
  283. 2 * kappa * BESSEL_J1(kappa * diffModule(t, tau)) / diffModule(t, tau)) / 2.0 + &
  284. kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * BESSEL_J1(kappa * diffModule(t, tau)) / (2.0 * diffModule(t, tau))
  285. return
  286. end
  287.  
  288. real(kind=8) function n1(t, tau)
  289. real(kind=8) t, tau
  290. if (t /= tau) then
  291. n1 = -tildaN(t, tau) * (kappa * kappa * BESSEL_J0(kappa * diffModule(t, tau)) - 2 * kappa * BESSEL_J1(kappa * &
  292. diffModule(t, tau)) / diffModule(t, tau)) / (2.0 * pi) - &
  293. kappa * BESSEL_J1(kappa * diffModule(t, tau)) * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) / &
  294. (2.0 * pi * diffModule(t, tau))
  295. return
  296. else
  297. n1 = -kappa * kappa * derModule(t) * derModule(t) / (4.0 * pi)
  298. return
  299. end if
  300. end
  301.  
  302. real(kind=8) function n2re(t, tau)
  303. real(kind=8) t, tau
  304. if (t /= tau) then
  305. n2re = Nre(t, tau) - N1(t, tau) * log(4.0 * sin((t - tau) / 2.0) * sin((t - tau) / 2.0))
  306. return
  307. else
  308. n2re = -kappa * kappa * derModule(t) * derModule(t) / (4.0 * pi) - eulerconst * kappa * kappa * derModule(t) &
  309. * derModule(t) / (2.0 * pi) - kappa * kappa * derModule(t) * derModule(t) * &
  310. log(kappa * derModule(t) / 2.0) / (2.0 * pi) + 1 / (12.0 * pi) &
  311. + (x1d(t) * x1dd(t) + x2d(t) * x2dd(t)) * (x1d(t) * x1dd(t) + x2d(t) * x2dd(t)) / &
  312. (2.0 * pi * derModule(t)*derModule(t)*derModule(t)*derModule(t)) &
  313. - doubleDerModule(t) * doubleDerModule(t) / (4.0 * pi * derModule(t) * derModule(t)) &
  314. - (x1d(t) * x1ddd(t) + x2d(t) * x2ddd(t)) / (6.0 * pi * derModule(t) * derModule(t))
  315. return
  316. end if
  317. end
  318.  
  319. real(kind=8) function n2im(t, tau)
  320. real(kind=8) t, tau
  321. if (t /= tau) then
  322. n2im = Nim(t, tau)
  323. return
  324. else
  325. n2im = kappa * kappa * derModule(t) * derModule(t) / 4.0
  326. return
  327. end if
  328. end
  329.  
  330. real(kind=8) function k1re(t, tau)
  331. real(kind=8) t, tau
  332. k1re = kappa * kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * M1(t, tau) - N1(t, tau)
  333. return
  334. end
  335.  
  336. real(kind=8) function k1im(t, tau)
  337. real(kind=8) t, tau
  338. k1im = -eta * H1(t, tau)
  339. return
  340. end
  341.  
  342. real(kind=8) function k2re(t, tau)
  343. real(kind=8) t, tau
  344. k2re = kappa * kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * M2re(t, tau) - N2re(t, tau) + eta * H2im(t, tau)
  345. return
  346. end
  347.  
  348. real(kind=8) function k2im(t, tau)
  349. real(kind=8) t, tau
  350. k2im = kappa * kappa * (x1d(t) * x1d(tau) + x2d(t) * x2d(tau)) * M2im(t, tau) - N2im(t, tau) - eta * H2re(t, tau)
  351. return
  352. end
  353.  
  354. real(kind=8) function tcapital(j)
  355. integer j
  356. if (j == 0) then
  357. tcapital = -n / 2.0
  358. return
  359. else if (mod(j,2) == 0) then
  360. tcapital = 0.0
  361. return
  362. else
  363. tcapital = 1 / (2.0 * n * sin(t(j) / 2.0) * sin(t(j) / 2.0))
  364. return
  365. end if
  366. end
  367.  
  368. real(kind=8) function t(j)
  369. integer j
  370. t = j * pi / n
  371. end
  372.  
  373. real(kind=8) function r(j)
  374. integer j, m
  375. real(kind=8) sum
  376. sum = 0.0
  377. do m = 1, n-1
  378. sum = sum + cos(m * j * pi / n) / m
  379. end do
  380. r = -2 * pi * sum / n + (-1)**j * pi / (n * n)
  381. end
  382.  
  383. real(kind=8) function aim(t)
  384. real(kind=8) t
  385. aim = eta * derModule(t)
  386. return
  387. end
  388.  
  389. real(kind=8) function lre(k, j)
  390. integer k, j
  391. lre = tcapital(abs(k - j)) + r(abs(k - j)) * K1re(t(k), t(j)) + pi * K2re(t(k), t(j)) / n
  392. return
  393. end
  394.  
  395. real(kind=8) function lim(k, j)
  396. integer k, j
  397. lim = r(abs(k - j)) * K1im(t(k), t(j)) + pi * K2im(t(k), t(j)) / n
  398. return
  399. end
  400.  
  401. real(kind=8) function uinfre(xhat, density)
  402. real(kind=8), dimension(:) :: xhat, density
  403. real(kind=8) :: sum = 0.0
  404. integer j
  405. do j = 0, 2 * n - 1
  406. sum = sum + (kappa * (xHat(1) * x2d(t(j)) - xHat(2) * x1d(t(j))) + eta) * derModule(t(j)) * &
  407. Ere(xHat, j, density)/(derModule(t(j))* derModule(t(j)))
  408. end do
  409. uinfre = pi * sum / (n * SQRT(8 * pi * kappa))
  410. return
  411. end
  412.  
  413. real(kind=8) function ere(xhat, j, density)
  414. real(kind=8), dimension(:) :: xhat, density
  415. real(kind=8) scalarprod
  416. integer j
  417. scalarprod = (xHat(1) * x1(t(j)) + xHat(2) * x2(t(j))) / SQRT(xHat(1) * xHat(1) + xHat(2) * xHat(2))
  418. ere = cos(kappa * scalarprod + pi / 4.0) * density(j+1) + sin(kappa * scalarprod + pi / 4.0) * density(j + 2 * n+1)
  419. return
  420. end
  421.  
  422. real(kind=8) function uinfim(xhat, density)
  423. real(kind=8), dimension(:) :: xhat, density
  424. real(kind=8) :: sum = 0.0
  425. integer j
  426. do j = 0, 2 * n - 1
  427. sum = sum + (kappa * (xHat(1) * x2d(t(j)) - xHat(2) * x1d(t(j))) + eta) * (derModule(t(j))) * Eim(xHat, j, density)
  428. end do
  429. uinfim = pi * sum / (n * SQRT(8 * pi * kappa))
  430. return
  431. end
  432.  
  433. real(kind=8) function eim(xhat, j, density)
  434. real(kind=8), dimension(:) :: xhat, density
  435. real(kind=8) scalarprod
  436. integer j
  437. scalarprod = (xHat(1) * x1(t(j)) + xHat(2) * x2(t(j))) / SQRT(xHat(1) * xHat(1) + xHat(2) * xHat(2))
  438. eim = -sin(kappa * scalarprod + pi / 4.0) * density(j+1) + cos(kappa * scalarprod + pi / 4.0) * density(j + 2 * n + 1)
  439. return
  440. end
  441.  
  442. real(kind=8) function fre(t)
  443. real(kind=8) t
  444. real(kind=8), dimension(2) :: x
  445. x(1) = x1(t)
  446. x(2) = x2(t)
  447. fre = kappa * BESSEL_Y1(kappa * module(x, yStar)) * ((x1(t) - yStar(1)) * x2d(t) - (x2(t) - yStar(2)) * x1d(t)) / &
  448. (2 * module(x, yStar));
  449. end
  450.  
  451. real(kind=8) function fim(t)
  452. real(kind=8) t
  453. real(kind=8), dimension(2) :: x
  454. x(1) = x1(t)
  455. x(2) = x2(t)
  456. fim = -kappa * BESSEL_J1(kappa * module(x, yStar)) * ((x1(t) - yStar(1)) * x2d(t) - (x2(t) - yStar(2)) * x1d(t)) / &
  457. (2 * module(x, yStar));
  458. end
  459.  
  460. real(kind=8) function module(x, y)
  461. real(kind=8), dimension(:) :: x,y
  462. module = Sqrt((x(1) - y(1)) * (x(1) - y(1)) + (x(2) - y(2)) * (x(2) - y(2)))
  463. return
  464. end
  465.  
  466. real(kind=8) function resultRe(x, density)
  467. real(kind=8), dimension(:) :: x,density
  468. real(kind=8), dimension(2) :: temp
  469. real(kind=8) sum ,h
  470. integer j
  471. sum = 0
  472. do j=0,2*n-1
  473. temp(1) = x1(t(j))
  474. temp(2) = x2(t(j))
  475. h = ((x1(t(j)) - x(1)) * x2d(t(j)) - (x2(t(j)) - x(2)) * x1d(t(j))) / (module(x, temp) * derModule(t(j)));
  476. sum = sum + (density(j+1) * (kappa * BESSEL_Y1(kappa * module(x, temp)) * h + eta * BESSEL_J0(kappa * &
  477. module(x, temp))) - density(j + 2 * n+1) * (eta * BESSEL_Y0(kappa * module(x, temp)) - &
  478. kappa * BESSEL_J1(kappa * module(x, temp)) * h)) * derModule(t(j))
  479. end do
  480. resultRe = pi * sum / (4 * n)
  481. return
  482. end
  483.  
  484. real(kind=8) function exactRe(x)
  485. real(kind=8), dimension(:) :: x
  486. exactRe =-BESSEL_Y0(kappa*module(x, yStar))/4
  487. return
  488. end
  489.  
  490. real(kind=8) function resultIm(x, density)
  491. real(kind=8), dimension(:) :: x,density
  492. real(kind=8), dimension(2) :: temp
  493. real(kind=8) sum ,h
  494. integer j
  495. sum = 0
  496. do j=0,2*n-1
  497. temp(1) = x1(t(j))
  498. temp(2) = x2(t(j))
  499. h = ((x1(t(j)) - x(1)) * x2d(t(j)) - (x2(t(j)) - x(2)) * x1d(t(j))) / (module(x, temp) * derModule(t(j)));
  500. sum = sum + (density(j+1) * (eta * BESSEL_Y0(kappa * module(x, temp)) - &
  501. kappa * BESSEL_J1(kappa * module(x, temp)) * h) &
  502. + density(j + 2*n + 1) * (kappa * BESSEL_Y1(kappa * module(x, temp)) * h + &
  503. eta * BESSEL_J0(k * module(x, temp)))) * derModule(t(j))
  504. end do
  505. resultIm = pi * sum / (4 * n)
  506. return
  507. end
  508.  
  509. real(kind=8) function exactIm(x)
  510. real(kind=8), dimension(:) :: x
  511. exactIm = BESSEL_J0(kappa*module(x, yStar))/4
  512. return
  513. end
  514.  
  515.  
  516.  
  517. ! PARA
  518.  
  519.  
  520. function solve(A, b) result(x)
  521.  
  522. integer, parameter :: master_id = 0
  523. integer :: id, tag, numprocs
  524. integer :: mpi_error_code
  525. integer, dimension(MPI_STATUS_SIZE) :: status
  526.  
  527. real(kind=8), dimension(:, :) :: A
  528. real(kind=8), dimension(:) :: b
  529.  
  530. real(kind=8), dimension(:), allocatable :: x
  531.  
  532. real(kind=8) :: upi
  533. integer :: i,j,n,p
  534. real(kind=8), dimension(:, :), allocatable :: u(:,:)
  535.  
  536. n = size(A, 1)
  537.  
  538. u = reshape( [A, b], [n,n+1] )
  539. do j = 1, n
  540.  
  541. if (id == master_id) then
  542. p = maxloc(abs(u(j:n, j)),1) + j - 1 ! maxloc returns indices between (1,n-j+1)
  543. if (p /= j) u([p, j],j) = u([j, p],j)
  544. u(j + 1:,j) = u(j + 1:,j)/u(j, j)
  545. end if
  546.  
  547. !call MPI_comm_rank(MPI_COMM_WORLD, id, mpi_error_code)
  548. !call MPI_comm_size(MPI_COMM_WORLD, numprocs, mpi_error_code)
  549.  
  550. !call MPI_bcast (u(j + 1:n, j), n - j, MPI_DOUBLE_PRECISION, master_id, &
  551. ! MPI_COMM_WORLD, mpi_error_code)
  552. !if (id == master_id) then
  553. !do i = j + 1, n
  554. !tag = i
  555. !call MPI_send(u(j + 1:n, i), n - j, MPI_DOUBLE_PRECISION, &
  556. ! mod(i - j, numprocs - 1) + 1, tag, &
  557. ! MPI_COMM_WORLD, mpi_error_code)
  558. !write(*, *) "SENT", id, u(j + 1:n, i)
  559. !tag = i + n
  560. !call MPI_recv(u(j + 1:n, i), n - j, MPI_DOUBLE_PRECISION, &
  561. ! mod(i - j, numprocs - 1) + 1, tag, &
  562. ! MPI_COMM_WORLD, status, mpi_error_code)
  563. !write(*, *) "RECV", id, u(j + 1:n, i)
  564. !end do
  565. !else
  566. do i = j + 1, n !j + id, n, numprocs - 1
  567. !tag = i
  568. !call MPI_recv(u(j + 1:n, i), n - j, MPI_DOUBLE_PRECISION, master_id, tag, &
  569. ! MPI_COMM_WORLD, status, mpi_error_code)
  570. !write(*, *) "RECV", id, u(j + 1:n, i)
  571. upi = u(p, i)
  572. if (p /= j) u([p, j],i) = u([j, p],i)
  573. u(j + 1:n,i) = u(j + 1:n, i) - upi * u(j + 1:n, j)
  574. !tag = i + n
  575. !call MPI_send(u(j + 1:n, i), n - j, MPI_DOUBLE_PRECISION, master_id, tag, &
  576. ! MPI_COMM_WORLD, mpi_error_code)
  577. !write(*, *) "SENT", id, u(j + 1:n, i)
  578. end do
  579. !end if
  580.  
  581. end do
  582.  
  583. allocate(x(n))
  584.  
  585. !if (id == master_id) then
  586. do i = n, 1, -1
  587. x(i) = (u(i,n+1) - sum(u(i,i+1:n)*x(i+1:n))) / u(i,i)
  588. end do ! i
  589. !end if
  590. !call MPI_bcast (x, n, MPI_DOUBLE_PRECISION, master_id, &
  591. ! MPI_COMM_WORLD, mpi_error_code)
  592. end
  593.  
  594.  
  595. ! PARA
  596.  
  597. end program Helmholtz
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!