From 4f2b1b37fa259a4afcc3e301b1792df2c9760afe Mon Sep 17 00:00:00 2001
From: schnetter <schnetter@f69c4107-0314-4c4f-9ad4-17e986b73f4a>
Date: Tue, 8 Apr 2008 23:15:49 +0000
Subject: Parallelise with OpenMP.  Remove function outer_prod_3d, ad write its
 content explicitly, since this is much easier to parallelise (and also
 shorter).

git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@99 f69c4107-0314-4c4f-9ad4-17e986b73f4a
---
 src/set_norm_mask.F90 | 34 ++++++++++++----------------------
 1 file changed, 12 insertions(+), 22 deletions(-)

diff --git a/src/set_norm_mask.F90 b/src/set_norm_mask.F90
index 267d7f8..a7eab31 100644
--- a/src/set_norm_mask.F90
+++ b/src/set_norm_mask.F90
@@ -44,6 +44,7 @@ subroutine SBP_SetNormMask (CCTK_ARGUMENTS)
                                                      
   CCTK_REAL, dimension(8) :: bmask
   CCTK_INT :: ni, nj, nk
+  CCTK_INT :: i, j, k
 
   ni = cctk_lsh(1); nj = cctk_lsh(2); nk = cctk_lsh(3)
 
@@ -51,7 +52,7 @@ subroutine SBP_SetNormMask (CCTK_ARGUMENTS)
 
 
   if (any (onesided/=0)) then
-    allocate ( mask_x(ni), mask_y(nj), mask_z(nk) )
+     allocate ( mask_x(ni), mask_y(nj), mask_z(nk) )
     mask_x = 1.0d0; mask_y = 1.0d0; mask_z = 1.0d0
 
     if ( CCTK_EQUALS(norm_type,'Diagonal') ) then
@@ -85,29 +86,18 @@ subroutine SBP_SetNormMask (CCTK_ARGUMENTS)
     if (onesided(5)/=0) mask_z(1:np) = bmask(1:np)
     if (onesided(6)/=0) mask_z(nk:nk-np+1:-1) = bmask(1:np)
 
-    nmask = outer_prod_3d ( mask_x, mask_y, mask_z )
-    deallocate ( mask_x, mask_y, mask_z )
+!$omp parallel workshare
+    forall (i=1:ni, j=1:nj, k=1:nk)
+       nmask(i,j,k) = mask_x(i) * mask_y(j) * mask_z(k)
+    end forall
+!$omp end parallel workshare
+
   else
+
+!$omp parallel workshare
     nmask = 1.0_wp
-  end if
+!$omp end parallel workshare
 
-contains
-
-  function outer_prod_3d ( ax, ay, az )
-    CCTK_REAL, dimension(:), intent(IN) :: ax
-    CCTK_REAL, dimension(:), intent(IN) :: ay
-    CCTK_REAL, dimension(:), intent(IN) :: az
-    CCTK_REAL, dimension(size(ax),size(ay),size(az)) :: outer_prod_3d
-    CCTK_INT :: nx, ny, nz, i, j, k
-  
-    nx = size(ax); ny = size(ay); nz = size(az)
-    do k = 1, nz
-      do j = 1, ny
-        do i = 1, nx
-          outer_prod_3d(i,j,k) = ax(i) * ay(j) * az(k)
-        end do
-      end do
-    end do
-  end function outer_prod_3d
+  end if
 
 end subroutine SBP_SetNormMask
-- 
cgit v1.2.3