From d19ab4c4f62977bcfafccda895656f35774ba126 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 25 Jul 2021 13:38:56 +0200
Subject: [PATCH] co concurrent possible here

note: do concurrent cannot be used for the double loops for 66 to 3333!
---
 PRIVATE      |  2 +-
 src/math.f90 | 30 +++++++++++++++++++-----------
 2 files changed, 20 insertions(+), 12 deletions(-)

diff --git a/PRIVATE b/PRIVATE
index 174ecac2d..d92b030e5 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378
+Subproject commit d92b030e5777e718c77edc2e1e93abfa0981b024
diff --git a/src/math.f90 b/src/math.f90
index 110af5128..89ee7abfe 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -260,9 +260,9 @@ pure function math_identity4th()
   integer :: i,j,k,l
 
 
-  do i=1,3; do j=1,3; do k=1,3; do l=1,3
+  do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
     math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
-  enddo; enddo; enddo; enddo
+  enddo
 
 end function math_identity4th
 
@@ -329,9 +329,10 @@ pure function math_outer(A,B)
   real(pReal), dimension(size(A,1),size(B,1)) ::  math_outer
   integer :: i,j
 
-  do i=1,size(A,1); do j=1,size(B,1)
+
+  do concurrent(i=1:size(A,1), j=1:size(B,1))
     math_outer(i,j) = A(i)*B(j)
-  enddo; enddo
+  enddo
 
 end function math_outer
 
@@ -371,9 +372,10 @@ pure function math_mul3333xx33(A,B)
   real(pReal), dimension(3,3) :: math_mul3333xx33
   integer :: i,j
 
-  do i=1,3; do j=1,3
+
+  do concurrent(i=1:3, j=1:3)
     math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
-  enddo; enddo
+  enddo
 
 end function math_mul3333xx33
 
@@ -388,9 +390,9 @@ pure function math_mul3333xx3333(A,B)
   real(pReal), dimension(3,3,3,3), intent(in) :: B
   real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333
 
-  do i=1,3; do j=1,3; do k=1,3; do l=1,3
+  do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
     math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
-  enddo; enddo; enddo; enddo
+  enddo
 
 end function math_mul3333xx3333
 
@@ -710,6 +712,7 @@ pure function math_6toSym33(v6,weighted)
   real(pReal), dimension(6) :: w
   integer :: i
 
+
   if(present(weighted)) then
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
@@ -734,6 +737,7 @@ pure function math_3333to99(m3333)
 
   integer :: i,j
 
+
   do concurrent(i=1:9, j=1:9)
     math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
   enddo
@@ -751,6 +755,7 @@ pure function math_99to3333(m99)
 
   integer :: i,j
 
+
   do concurrent(i=1:9, j=1:9)
     math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
   enddo
@@ -773,15 +778,16 @@ pure function math_sym3333to66(m3333,weighted)
   real(pReal), dimension(6) :: w
   integer :: i,j
 
+
   if(present(weighted)) then
     w = merge(NRMMANDEL,1.0_pReal,weighted)
   else
     w = NRMMANDEL
   endif
 
-  do i=1,6; do j=1,6
+  do concurrent(i=1:6, j=1:6)
     math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
-  enddo; enddo
+  enddo
 
 end function math_sym3333to66
 
@@ -801,6 +807,7 @@ pure function math_66toSym3333(m66,weighted)
   real(pReal), dimension(6) :: w
   integer :: i,j
 
+
   if(present(weighted)) then
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
@@ -826,6 +833,7 @@ pure function math_Voigt66to3333(m66)
   real(pReal), dimension(6,6), intent(in) :: m66                                                    !< 6x6 matrix
   integer :: i,j
 
+
   do i=1,6; do j=1, 6
     math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
     math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
@@ -1240,7 +1248,7 @@ subroutine selfTest
   if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
     error stop 'math_det33/math_detSym33'
 
-  if(any(dNeq0(t33+transpose(t33)-math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) &
+  if(any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) &
     error stop 'math_mul3333xx33/math_identity4th'
 
   if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) &