* add OMP pragmas to EOSOmni routines

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/trunk@83 8e189c6b-2ab8-4400-aa02-70a9cfce18b9

2 files changed, 115 insertions, 15 deletions

diff --git a/src/EOS_Omni_MultiVarCalls.F90 b/src/EOS_Omni_MultiVarCalls.F90
index e6a37b3..5db911d 100644
--- a/src/EOS_Omni_MultiVarCalls.F90
+++ b/src/EOS_Omni_MultiVarCalls.F90
@@ -46,7 +46,9 @@ subroutine EOS_Omni_EOS_short(eoskey,keytemp,rf_precision,npoints,&
 
   anyerr    = 0
   keyerr(:) = 0
-  
+
+  !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xent,xcs2,&
+  !$OMP                    xdedt,xdpderho,xdpdrhoe,xmunu)
   do i=1,npoints
 
      xrho = rho(i) * inv_rho_gf
@@ -59,7 +61,9 @@ subroutine EOS_Omni_EOS_short(eoskey,keytemp,rf_precision,npoints,&
           keytemp,keyerr(i),rf_precision)
      
      if(keyerr(i).ne.0) then
+        !$OMP CRITICAL
         anyerr = 1
+        !$OMP END CRITICAL
      endif
      
      if(keytemp.eq.1) then
@@ -80,6 +84,7 @@ subroutine EOS_Omni_EOS_short(eoskey,keytemp,rf_precision,npoints,&
      munu(i)      = xmunu
      
   enddo
+  !$OMP END PARALLEL DO
 
 end subroutine EOS_Omni_EOS_short
 
@@ -111,7 +116,6 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: xprs,xmunu,xcs2
   real*8           :: xdedt,xdpderho,xdpdrhoe
 
-
   anyerr    = 0
   keyerr(:) = 0
 
@@ -119,35 +123,45 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
   case (1)
      ! polytropic EOS                                                                                      
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            dpdrhoe(i) = poly_k *  &
                 poly_gamma *        &
                 rho(i) ** (poly_gamma - 1.d0)
            dpderho(i) = 0.0d0
         enddo
+        !$OMP END PARALLEL DO
   case (2)
      ! gamma-law EOS                                                                                       
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            dpdrhoe(i) = (gl_gamma-1.0d0) * &
                 eps(i)
            dpderho(i) = (gl_gamma - 1.0d0) * &
                 rho(i)
         enddo
+        !$OMP END PARALLEL DO
   case (3)
      ! hybrid EOS                                                                                          
+     !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
+     !$OMP                     hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2)
         do i=1,npoints
            if(rho(i).gt.hybrid_rho_nuc) then
               hybrid_local_gamma = hybrid_gamma2
@@ -173,7 +187,9 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
            dpdrhoe(i) = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2
            dpderho(i) = (hybrid_gamma_th - 1.0d0) * rho(i)
         enddo
+        !$OMP END PARALLEL DO
     case (4)
+       !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xdpderho,xdpdrhoe)
         do i=1,npoints
            xrho = rho(i) * inv_rho_gf
            xtemp = temp(i)
@@ -184,7 +200,9 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
                 keytemp,keyerr(i),rf_precision)
 
            if(keyerr(i).ne.0) then
+              !$OMP CRITICAL
               anyerr = 1
+              !$OMP END CRITICAL
            endif
 
            if(keytemp.eq.1) then
@@ -239,34 +257,45 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
   case (1)
      ! polytropic EOS                                                                                      
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            depsdpress(i) = 1.0d0/(poly_gamma - 1.0d0)/rho(i)
            depsdrho(i) = depsdpress(i) * poly_k *  &
                 poly_gamma *         &
                 rho(i) ** (poly_gamma - 1.d0)
         enddo
+        !$OMP END PARALLEL DO
   case (2)
      ! gamma-law EOS                                                                                       
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            depsdpress(i) = 1.0/( (gl_gamma - 1.0d0) * &
                 rho(i))
            depsdrho(i) = -eps(i)/rho(i)
        enddo
+       !$OMP END PARALLEL DO
   case (3)
      ! hybrid EOS                                                                                          
+     !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
+     !$OMP                    hybrid_dp_poly,hybrid_dp_th1,&
+     !$OMP                    hybrid_dp_th2,xdpdrhoe,xdpderho)
         do i=1,npoints
            if(rho(i).gt.hybrid_rho_nuc) then
               hybrid_local_gamma = hybrid_gamma2
diff --git a/src/EOS_Omni_SingleVarCalls.F90 b/src/EOS_Omni_SingleVarCalls.F90
index fc7a90b..fead6cb 100644
--- a/src/EOS_Omni_SingleVarCalls.F90
+++ b/src/EOS_Omni_SingleVarCalls.F90
@@ -47,31 +47,41 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
      case (1)
         ! polytropic EOS
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) =  poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            press(i) = poly_k * &
                 rho(i)**poly_gamma
         enddo
+        !$OMP END PARALLEL DO
      case (2)
         ! gamma-law EOS
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            press(i) = (gl_gamma - 1.0d0) * rho(i) * eps(i)
         enddo
+        !$OMP END PARALLEL DO
 
      case (3)
         ! hybrid EOS
+        !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
+        !$OMP                     hybrid_p_poly, hybrid_p_th)                          
         do i=1,npoints
            if(rho(i).gt.hybrid_rho_nuc) then
               hybrid_local_gamma = hybrid_gamma2
@@ -92,11 +102,12 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
            hybrid_p_th = max(zero, hybrid_p_th)
            press(i) = hybrid_p_poly + hybrid_p_th
         enddo
+        !$OMP END PARALLEL DO
 
      case (4)
         ! nuc eos
+        !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs)
         do i=1,npoints
-
            xrho = rho(i) * inv_rho_gf
            xtemp = temp(i)
            xye = ye(i)
@@ -105,7 +116,9 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
                 keytemp,keyerr(i),rf_precision)
            
            if(keyerr(i).ne.0) then
+              !$OMP CRITICAL
               anyerr = 1
+              !$OMP END CRITICAL
            endif
 
            if(keytemp.eq.1) then
@@ -115,11 +128,13 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
            endif
 
            press(i) = xprs * press_gf
-
         enddo
+        !$OMP END PARALLEL DO
 
      case (5)
         ! cold tabular EOS with gamma law
+        !$OMP PARALLEL DO PRIVATE(xrho,ir,gamma,eps_cold,eps_th, &
+        !$OMP                     press_cold,press_th)
         do i=1,npoints
            if(rho(i).lt.coldeos_rhomin) then
               press(i) = coldeos_low_kappa * rho(i)**coldeos_low_gamma
@@ -127,7 +142,9 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
               cycle
            else if(rho(i).gt.coldeos_rhomax) then
               keyerr(i) = 103
+              !$OMP CRITICAL
               anyerr = 1
+              !$OMP END CRITICAL
            else
               xrho = log10(rho(i))
               ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
@@ -153,6 +170,7 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
            press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th
            press(i) = press_cold + press_th
         enddo
+        !$OMP END PARALLEL DO
 
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
@@ -201,32 +219,43 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
      case (1)
         ! polytropic EOS
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            dpdrhoe(i) = poly_k *  &
                 poly_gamma *         & 
                 rho(i) ** (poly_gamma - 1.d0)
         enddo
+        !$OMP END PARALLEL DO
      case (2)
         ! gamma-law EOS
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            dpdrhoe(i) = (gl_gamma-1.0d0) * &
                 eps(i)
         enddo
+        !$OMP END PARALLEL DO
      case (3)
         ! hybrid EOS
+        !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
+        !$OMP                     hybrid_dp_poly, hybrid_dp_th1, &
+        !$OMP                     hybrid_dp_th2)
         do i=1,npoints
            if(rho(i).gt.hybrid_rho_nuc) then
               hybrid_local_gamma = hybrid_gamma2
@@ -252,8 +281,10 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
            dpdrhoe(i) = hybrid_dp_poly + max(0.0d0,hybrid_dp_th1 + hybrid_dp_th2)
 
         enddo
-
+        !$OMP END PARALLEL DO
      case (4)
+        !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
+        !$OMP                     xdpderho,xdpdrhoe,xmunu,h)
         do i=1,npoints
            xrho = rho(i) * inv_rho_gf
            xtemp = temp(i)
@@ -264,7 +295,9 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
                 keytemp,keyerr(i),rf_precision)
            
            if(keyerr(i).ne.0) then
+              !$OMP CRITICAL
               anyerr = 1
+              !$OMP END CRITICAL
            endif
 
            if(keytemp.eq.1) then
@@ -275,11 +308,13 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
 
            dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf
         enddo
-        
+        !$OMP END PARALLEL DO
      case (5)
         ! with the cold eos we have to assume P = P(rho), so
         ! by definition dPdrho is at constant internal energy
         ! and entropy (the latter, because T=0)
+        !$OMP PARALLEL DO PRIVATE(xrho,ir,gamma,eps_cold,eps_th, &
+        !$OMP                     cs2,press_cold,press_th)
         do i=1,npoints
            if(rho(i).lt.coldeos_rhomin) then
               dpdrhoe(i) = coldeos_low_kappa * coldeos_low_gamma * &
@@ -287,7 +322,9 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
               cycle
            else if(rho(i).gt.coldeos_rhomax) then
               keyerr(i) = 103
+              !$OMP CRITICAL
               anyerr = 1
+              !$OMP END CRITICAL
            else
               xrho = log10(rho(i))
               ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
@@ -364,38 +401,49 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
      case (1)
         ! polytropic EOS
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            dpdepsrho(i) = 0.0d0
         enddo
+        !$OMP END PARALLEL DO
      case (2)
         ! gamma-law EOS
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            dpdepsrho(i) = (gl_gamma - 1.0d0) * &
                 rho(i)
         enddo
+        !$OMP END PARALLEL DO
      case (3)
         ! hybrid EOS
+        !$OMP PARALLEL DO
         do i=1,npoints
            dpdepsrho(i) = (hybrid_gamma_th - 1.0d0) * rho(i)
         enddo
-
+        !$OMP END PARALLEL DO
      case (4)
         ! nuc_eos
+        !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs,&
+        !$OMP                     xent,xcs2,xdedt,xdpderho,xdpdrhoe,&
+        !$OMP                     xmunu)
         do i=1,npoints
-
            xrho = rho(i) * inv_rho_gf
            xtemp = temp(i)
            xye = ye(i)
@@ -405,7 +453,9 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
                 keytemp,keyerr(i),rf_precision)
            
            if(keyerr(i).ne.0) then
+              !$OMP CRITICAL
               anyerr = 1
+              !$OMP END CRITICAL
            endif
 
            if(keytemp.eq.1) then
@@ -417,17 +467,20 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
            dpdepsrho(i) = xdpderho * press_gf * inv_eps_gf 
 
         enddo
-
+        !$OMP END PARALLEL DO
      case (5)
         ! with the cold eos we have to assume P = P(rho), so
         ! only the gamma law has non-zero dPdeps
+        !$OMP PARALLEL DO PRIVATE(xrho,ir,eps_cold,eps_th)
         do i=1,npoints
            if(rho(i).lt.coldeos_rhomin) then
               dpdepsrho(i) = 0.0d0
               cycle
            else if(rho(i).gt.coldeos_rhomax) then
               keyerr(i) = 103
+              !$OMP CRITICAL
               anyerr = 1
+              !$OMP END CRITICAL
            else
               xrho = log10(rho(i))
               ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
@@ -494,27 +547,34 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
      case (1)
         ! polytropic EOS
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO
         do i=1,npoints
            xpress = poly_k * &
                 rho(i)**(poly_gamma)
            cs2(i) = poly_gamma * xpress / rho(i) / &
                 (1 + eps(i) + xpress/rho(i))
         enddo
+        !$OMP END PARALLEL DO
      case (2)
         ! gamma-law EOS
         if(keytemp.eq.1) then
+           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
+           !$OMP END PARALLEL DO
         endif
+        !$OMP PARALLEL DO PRIVATE(xpress,xdpdrhoe,xdpderho)
         do i=1,npoints
            xpress = (gl_gamma-1.0d0)*rho(i)*eps(i)
            xdpdrhoe = (gl_gamma-1.0d0)*eps(i)
@@ -522,8 +582,12 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
            cs2(i) = (xdpdrhoe + xpress * xdpderho / (rho(i)**2)) / &
                 (1.0d0 + eps(i) + xpress/rho(i))
         enddo
+        !$OMP END PARALLEL DO
      case(3)
         ! hybrid EOS
+        !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
+        !$OMP                     hybrid_p_poly,hybrid_p_th,&
+        !$OMP                     xpress)
         do i=1,npoints
            if(rho(i).gt.hybrid_rho_nuc) then
               hybrid_local_gamma = hybrid_gamma2
@@ -547,11 +611,12 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
            cs2(i) = (hybrid_local_gamma * hybrid_p_poly + hybrid_gamma_th * hybrid_p_th) / &
                 rho(i) / (1.0d0 + eps(i) + xpress/rho(i))
         enddo
+        !$OMP END PARALLEL DO
      case(4)
         ! nuc_eos
-
+        !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs,xent,xcs2,&
+        !$OMP                     xdedt,xdpderho,xdpdrhoe,xmunu)
         do i=1,npoints
-           
            xrho = rho(i) * inv_rho_gf
            xtemp = temp(i)
            xye = ye(i)
@@ -561,7 +626,9 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
                 keytemp,keyerr(i),rf_precision)
            
            if(keyerr(i).ne.0) then
+              !$OMP CRITICAL
               anyerr = 1
+              !$OMP END CRITICAL
            endif
 
            if(keytemp.eq.1) then
@@ -572,12 +639,15 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
 
            cs2(i) = xcs2 * cliteinv2 / &
                 (1.0d0 + eps(i) + xprs * press_gf / rho(i))
-
         enddo
+        !$OMP END PARALLEL DO
      case (5)
         ! with the cold eos we have to assume P = P(rho), so
         ! by definition dPdrho is at constant internal energy
         ! and entropy (the latter, because T=0)
+        !$OMP PARALLEL DO PRIVATE(xprs,xrho,ir,gamma,cs2_cold,&
+        !$OMP                     eps_cold,eps_th,press_cold,press_th,&
+        !$OMP                     xdpdrhoe,xdpderho,cs2_th,h,h_cold)
         do i=1,npoints
            if(rho(i).lt.coldeos_rhomin) then
               xprs = coldeos_low_kappa * rho(i)**coldeos_low_gamma
@@ -588,8 +658,10 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
               cycle
            else if(rho(i).gt.coldeos_rhomax) then
               keyerr(i) = 103
+              !$OMP CRITICAL
               anyerr = 1
-                 else
+              !$OMP END CRITICAL
+           else
               xrho = log10(rho(i))
               ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
            endif
@@ -614,7 +686,6 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
               eps_th = 0.0d0
               eps(i) = eps_cold
            endif
-
            press_cold = coldeos_kappa * rho(i)**gamma
            press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th
 
@@ -627,7 +698,7 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
            
            cs2(i) = (cs2_cold * h_cold + cs2_th) / h
         enddo
-           
+        !$OMP END PARALLEL DO
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)