parallelize newton-cotes (untestedgit add quadratures.f90)

1 files changed, 31 insertions, 7 deletions

diff --git a/src/quadratures.f90 b/src/quadratures.f90
index 9be7517..2eac963 100644
--- a/src/quadratures.f90
+++ b/src/quadratures.f90
@@ -47,12 +47,12 @@ CONTAINS
     integer(kind=4), intent(in) :: p
     
     real(kind=8) :: val, subinterval_width, qbeg, qend
-
+    real(kind=8), allocatable :: partval[:]
+    
     procedure(quadrature), pointer :: quad
-    integer(kind=8) :: i
+    integer(kind=8) :: min_i, max_i, i, subintervals_per_thread
+    integer(kind=4) :: im
 
-    subinterval_width = (iend - ibeg) / subintervals
-    
     SELECT CASE (p)
     CASE (:-1)
        STOP "negative interpolationg polynomial order passed"
@@ -66,14 +66,38 @@ CONTAINS
        STOP "Newton-Cotes quadratures only implemented for order < 3"
     END SELECT
 
+    if (this_image() == 1) allocate(partval[*])
+    
+    subintervals_per_thread =                               &
+         (subintervals + num_images() - 1) / num_images()
+    
+    min_i = subintervals_per_thread * (this_image() - 1) + 1
+    max_i = min(subintervals, subintervals_per_thread * this_image())
+
+    subinterval_width = (iend - ibeg) / subintervals
+    
     ! compute integral using quadrature pointed by quad
-    val = 0
+    partval = 0
     
-    DO i = 1, subintervals
+    DO i = min_i, max_i
        qend = ibeg + i * subinterval_width
        qbeg = ibeg + (i - 1) * subinterval_width
-       val = val + quad(qbeg, qend, fun)
+       partval = partval + quad(qbeg, qend, fun)
     END DO
+
+    IF (this_image() == 1 .and. num_images() > 1) THEN
+       sync images([(im, im = 2, num_images())])
+       partval = partval + sum([(partval[im], im = 2, num_images())])
+       sync images([(im, im = 2, num_images())])
+    END IF
+
+    IF (this_image() /= 1) THEN
+       sync images(1)
+       sync images(1)
+    END IF
+
+    val = partval[1]
+    
   END FUNCTION newton_cotes
 
   FUNCTION rectangle(qbeg, qend, fun) result(val)