cleaning

src/math.f90

@@ -6,12 +6,10 @@
 !> @brief Mathematical library, including random number generation and tensor representations
 !--------------------------------------------------------------------------------------------------
 module math
- use prec, only: &
+  use prec
-   pReal
   use future
 
   implicit none
- private
   real(pReal),    parameter, public :: PI = acos(-1.0_pReal)                                        !< ratio of a circle's circumference to its diameter
   real(pReal),    parameter, public :: INDEG = 180.0_pReal/PI                                       !< conversion from radian into degree
   real(pReal),    parameter, public :: INRAD = PI/180.0_pReal                                       !< conversion from degree into radian
@@ -76,79 +74,13 @@ module math
     math_mul3x3
 !---------------------------------------------------------------------------------------------------
 
-
- public :: &
-   math_init, &
-   math_sort, &
-   math_expand, &
-   math_range, &
-   math_identity2nd, &
-   math_identity4th, &
-   math_civita, &
-   math_delta, &
-   math_cross, &
-   math_outer, &
-   math_inner, &
-   math_mul33xx33, &
-   math_mul3333xx33, &
-   math_mul3333xx3333, &
-   math_exp33 , &
-   math_inv33, &
-   math_invert33, &
-   math_invSym3333, &
-   math_invert, &
-   math_invert2, &
-   math_symmetric33, &
-   math_symmetric66, &
-   math_skew33, &
-   math_spherical33, &
-   math_deviatoric33, &
-   math_equivStrain33, &
-   math_equivStress33, &
-   math_trace33, &
-   math_det33, &
-   math_33to9, &
-   math_9to33, &
-   math_sym33to6, &
-   math_6toSym33, &
-   math_3333to99, &
-   math_99to3333, &
-   math_sym3333to66, &
-   math_66toSym3333, &
-   math_Voigt66to3333, &
-   math_qMul, &
-   math_qConj, &
-   math_qRot, &
-   math_RtoEuler, &
-   math_RtoQ, &
-   math_EulerToR, &
-   math_axisAngleToR, &
-   math_qToRodrig, &
-   math_sampleGaussVar, &
-   math_eigenvectorBasisSym33, &
-   math_eigenvectorBasisSym33_log, &
-   math_eigenvectorBasisSym, &
-   math_eigenValuesVectorsSym33, &
-   math_eigenValuesVectorsSym, &
-   math_rotationalPart33, &
-   math_invariantsSym33, &
-   math_eigenvaluesSym33, &
-   math_factorial, &
-   math_binomial, &
-   math_multinomial, &
-   math_volTetrahedron, &
-   math_areaTriangle, &
-   math_rotate_forward33, &
-   math_rotate_backward33, &
-   math_rotate_forward3333, &
-   math_clip
  private :: &
    math_check
 
 contains
 
 !--------------------------------------------------------------------------------------------------
-!> @brief initialization of random seed generator
+!> @brief initialization of random seed generator and internal checks
 !--------------------------------------------------------------------------------------------------
 subroutine math_init
   use numerics, only: &
@@ -184,7 +116,7 @@ subroutine math_init
   write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest
 
   call random_seed(put = randInit)
- call math_check()
+  call math_check
 
 end subroutine math_init
  
@@ -192,7 +124,6 @@ end subroutine math_init
 !> @brief check correctness of (some) math functions
 !--------------------------------------------------------------------------------------------------
 subroutine math_check
- use prec,     only: tol_math_check
   use IO, only: IO_error
 
   character(len=64) :: error_msg
@@ -252,7 +183,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
     call math_sort(a, ipivot+1, e, d)
   endif
 
-!--------------------------------------------------------------------------------------------------
+
   contains
 
   !-------------------------------------------------------------------------------------------------
@@ -304,8 +235,7 @@ pure function math_expand(what,how)
   real(pReal), dimension(sum(how)) ::  math_expand
   integer :: i
  
- if (sum(how) == 0) &
+  if (sum(how) == 0) return
-   return
  
   do i = 1, size(how)
     math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1)
@@ -338,7 +268,9 @@ pure function math_identity2nd(dimen)
   real(pReal), dimension(dimen,dimen) :: math_identity2nd
 
   math_identity2nd = 0.0_pReal
- forall(i=1:dimen) math_identity2nd(i,i) = 1.0_pReal
+  do i=1, dimen
+    math_identity2nd(i,i) = 1.0_pReal
+  enddo
 
 end function math_identity2nd
 
@@ -484,99 +416,6 @@ pure function math_mul3333xx3333(A,B)
 end function math_mul3333xx3333
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication 33x33 = 33
-!--------------------------------------------------------------------------------------------------
-pure function math_mul33x33(A,B)
-
- real(pReal), dimension(3,3) ::  math_mul33x33
- real(pReal), dimension(3,3), intent(in) ::  A,B
- integer :: i,j
-
- forall(i=1:3,j=1:3) math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j)
-
-end function math_mul33x33
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication 66x66 = 66
-!--------------------------------------------------------------------------------------------------
-pure function math_mul66x66(A,B)
-
- real(pReal), dimension(6,6) ::  math_mul66x66
- real(pReal), dimension(6,6), intent(in) ::  A,B
- integer :: i,j
-
- forall(i=1:6,j=1:6) &
-   math_mul66x66(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) &
-                      + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j)
-
-end function math_mul66x66
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication 99x99 = 99
-!--------------------------------------------------------------------------------------------------
-pure function math_mul99x99(A,B)
-
- real(pReal), dimension(9,9) ::  math_mul99x99
- real(pReal), dimension(9,9), intent(in) ::  A,B
- integer  i,j
-
- forall(i=1:9,j=1:9) &
-   math_mul99x99(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) &
-                      + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) &
-                      + A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j)
-
-end function math_mul99x99
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication 33x3 = 3
-!--------------------------------------------------------------------------------------------------
-pure function math_mul33x3(A,B)
-
- real(pReal), dimension(3) ::  math_mul33x3
- real(pReal), dimension(3,3), intent(in) ::  A
- real(pReal), dimension(3),   intent(in) ::  B
- integer :: i 
-
- forall (i=1:3) math_mul33x3(i) = sum(A(i,1:3)*B)
-
-end function math_mul33x3
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication complex(33) x real(3) = complex(3)
-!--------------------------------------------------------------------------------------------------
-pure function math_mul33x3_complex(A,B)
-
- complex(pReal), dimension(3) ::  math_mul33x3_complex
- complex(pReal), dimension(3,3), intent(in) ::  A
- real(pReal),    dimension(3),   intent(in) ::  B
- integer :: i 
-
- forall (i=1:3) math_mul33x3_complex(i) = sum(A(i,1:3)*cmplx(B,0.0_pReal,pReal))
-
-end function math_mul33x3_complex
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication 66x6 = 6
-!--------------------------------------------------------------------------------------------------
-pure function math_mul66x6(A,B)
-
- real(pReal), dimension(6) ::  math_mul66x6
- real(pReal), dimension(6,6), intent(in) ::  A
- real(pReal), dimension(6),   intent(in) ::  B
- integer :: i
-
- forall (i=1:6) math_mul66x6(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) &
-                                          + A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
-
-end function math_mul66x6
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief 3x3 matrix exponential up to series approximation order n (default 5)
 !--------------------------------------------------------------------------------------------------
@@ -607,8 +446,6 @@ end function math_exp33
 ! if determinant is close to zero
 !--------------------------------------------------------------------------------------------------
 pure function math_inv33(A)
- use prec, only: &
-   dNeq0
 
   real(pReal),dimension(3,3),intent(in)  :: A
   real(pReal) :: DetA
@@ -644,8 +481,6 @@ end function math_inv33
 ! ToDo: Output arguments should be first
 !--------------------------------------------------------------------------------------------------
 pure subroutine math_invert33(A, InvA, DetA, error)
- use prec, only: &
-   dEq0
 
   logical, intent(out) :: error
   real(pReal),dimension(3,3),intent(in)  :: A
@@ -907,7 +742,9 @@ pure function math_33to9(m33)
   
   integer :: i
   
- forall(i=1:9) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
+  do i = 1, 9
+    math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
+  enddo
 
 end function math_33to9
 
@@ -922,7 +759,9 @@ pure function math_9to33(v9)
   
   integer :: i
 
- forall(i=1:9) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
+  do i = 1, 9
+    math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
+  enddo
 
 end function math_9to33
 
@@ -948,7 +787,9 @@ pure function math_sym33to6(m33,weighted)
     w = nrmMandel
   endif
 
- forall(i=1:6) math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i))
+  do i = 1, 6
+    math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i))
+  enddo 
 
 end function math_sym33to6
 
@@ -992,8 +833,9 @@ pure function math_3333to99(m3333)
 
   integer :: i,j
 
- forall(i=1:9,j=1:9) &
+  do i=1,9; do j=1,9
     math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j))
+  enddo; enddo
 
 end function math_3333to99
 
@@ -1008,8 +850,9 @@ pure function math_99to3333(m99)
 
   integer :: i,j
 
- forall(i=1:9,j=1:9) &
+  do i=1,9; do j=1,9
     math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
+  enddo; enddo
 
 end function math_99to3333
 
@@ -1035,8 +878,9 @@ pure function math_sym3333to66(m3333,weighted)
     w = nrmMandel
   endif
 
- forall(i=1:6,j=1:6) &
+  do i=1,6; do 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
 
 end function math_sym3333to66
 
@@ -1163,7 +1007,6 @@ pure function math_RtoEuler(R)
   squvw=sqrt(R(1,1)*R(1,1)+R(2,1)*R(2,1)+R(3,1)*R(3,1))
   sqhk =sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3))
 
-! calculate PHI
   math_RtoEuler(2) = acos(math_clip(R(3,3)/sqhkl,-1.0_pReal, 1.0_pReal))
 
   if((math_RtoEuler(2) < 1.0e-8_pReal) .or. (pi-math_RtoEuler(2) < 1.0e-8_pReal)) then
@@ -1316,10 +1159,6 @@ end function math_axisAngleToR
 !> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz)
 !--------------------------------------------------------------------------------------------------
 pure function math_qToRodrig(Q)
- use, intrinsic :: &
-   IEEE_arithmetic
- use prec, only: &
-   tol_math_check
 
   real(pReal), dimension(4), intent(in) :: Q
   real(pReal), dimension(3) :: math_qToRodrig
@@ -1333,8 +1172,6 @@ end function math_qToRodrig
 !> @brief draw a random sample from Gauss variable
 !--------------------------------------------------------------------------------------------------
 real(pReal) function math_sampleGaussVar(meanvalue, stddev, width)
- use prec, only: &
-   tol_math_check
 
   real(pReal), intent(in) ::            meanvalue, &                                                ! meanvalue of gauss distribution
                                         stddev                                                      ! standard deviation of gauss distribution
@@ -1595,8 +1432,6 @@ end function math_eigenvectorBasisSym33_log
 !> @brief rotational part from polar decomposition of 33 tensor m
 !--------------------------------------------------------------------------------------------------
 function math_rotationalPart33(m)
- use prec, only: &
-   dEq0
   use IO, only: &
     IO_warning
 
@@ -1622,8 +1457,6 @@ end function math_rotationalPart33
 ! will return NaN on error
 !--------------------------------------------------------------------------------------------------
 function math_eigenvaluesSym(m)
- use, intrinsic :: &
-   IEEE_arithmetic
  
   real(pReal), dimension(:,:),                  intent(in)  :: m
   real(pReal), dimension(size(m,1))                         :: math_eigenvaluesSym
@@ -1767,7 +1600,6 @@ end function math_areaTriangle
 !--------------------------------------------------------------------------------------------------
 pure function math_rotate_forward33(tensor,rot_tensor)
 
-
   real(pReal), dimension(3,3) ::  math_rotate_forward33
   real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
 
@@ -1815,8 +1647,6 @@ end function math_rotate_forward3333
 ! Will return NaN if left > right
 !--------------------------------------------------------------------------------------------------
 real(pReal) pure elemental function math_clip(a, left, right)
- use, intrinsic :: &
-   IEEE_arithmetic
 
   real(pReal), intent(in) :: a
   real(pReal), intent(in), optional :: left, right