diff --git a/src/CPFEM.f90 b/src/CPFEM.f90
index 42c41c50c..d34a79bf7 100644
--- a/src/CPFEM.f90
+++ b/src/CPFEM.f90
@@ -259,7 +259,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
    restartWrite
  use math, only: &
    math_identity2nd, &
-   math_mul33x33, &
    math_det33, &
    math_delta, &
    math_sym3333to66, &
@@ -557,7 +556,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
        endif
 
        ! translate from P to CS
-       Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
+       Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
        J_inverse  = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
        CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
 
diff --git a/src/DAMASK_grid.f90 b/src/DAMASK_grid.f90
index a1241f97f..663921084 100644
--- a/src/DAMASK_grid.f90
+++ b/src/DAMASK_grid.f90
@@ -337,7 +337,7 @@ program DAMASK_spectral
        endif
        enddo; write(6,'(/)',advance='no')
      enddo
-    if (any(abs(math_mul33x33(newLoadCase%rotation, &
+    if (any(abs(matmul(newLoadCase%rotation, &
                 transpose(newLoadCase%rotation))-math_I3) > &
                 reshape(spread(tol_math_check,1,9),[ 3,3]))&
                 .or. abs(math_det33(newLoadCase%rotation)) > &
diff --git a/src/constitutive.f90 b/src/constitutive.f90
index 96b29d846..4df97b240 100644
--- a/src/constitutive.f90
+++ b/src/constitutive.f90
@@ -381,8 +381,6 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
                                          S, Fi, ipc, ip, el)
  use prec, only: &
    pReal
- use math, only: &
-   math_mul33x33
  use material, only: &
    phasememberAt, &
    phase_plasticity, &
@@ -439,7 +437,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
  ho = material_homogenizationAt(el)
  tme = thermalMapping(ho)%p(ip,el)
 
- Mp  = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
+ Mp  = matmul(matmul(transpose(Fi),Fi),S)
 
  plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
 
@@ -483,9 +481,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
 #else
  do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
 #endif
-   dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
-                          math_mul33x33(math_mul33x33(Fi,dLp_dMp(i,j,1:3,1:3)),S)
-   dLp_dS(i,j,1:3,1:3)  = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3))       ! ToDo: @PS: why not:   dLp_dMp:(FiT Fi)
+   dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
+                          matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
+   dLp_dS(i,j,1:3,1:3)  = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3))       ! ToDo: @PS: why not:   dLp_dMp:(FiT Fi)
 #if defined(__INTEL_COMPILER) || defined(__PGI)
  end forall
 #else
@@ -506,8 +504,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
  use math, only: &
    math_I3, &
    math_inv33, &
-   math_det33, &
-   math_mul33x33
+   math_det33
  use material, only: &
    phasememberAt, &
    phase_plasticity, &
@@ -591,11 +588,11 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
 
  FiInv = math_inv33(Fi)
  detFi = math_det33(Fi)
- Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi                                               !< push forward to intermediate configuration
- temp_33 = math_mul33x33(FiInv,Li)
+ Li = matmul(matmul(Fi,Li),FiInv)*detFi                                               !< push forward to intermediate configuration
+ temp_33 = matmul(FiInv,Li)
 
  do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
-   dLi_dS(1:3,1:3,i,j)  = math_mul33x33(math_mul33x33(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
+   dLi_dS(1:3,1:3,i,j)  = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
    dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
    dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
  end do; end do
@@ -689,7 +686,6 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
  use prec, only: &
    pReal
  use math, only : &
-   math_mul33x33, &
    math_mul3333xx33, &
    math_66toSym3333, &
    math_I3
@@ -733,14 +729,14 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
    end select degradationType
  enddo DegradationLoop
 
- E = 0.5_pReal*(math_mul33x33(transpose(Fe),Fe)-math_I3)                                            !< Green-Lagrange strain in unloaded configuration
- S = math_mul3333xx33(C,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi))                           !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
+ E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3)                                            !< Green-Lagrange strain in unloaded configuration
+ S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi))                           !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
 
  dS_dFe = 0.0_pReal
  forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt)
    dS_dFe(i,j,1:3,1:3) = &
-     math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi)))                !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
-   dS_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3))                !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
+     matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi)))                !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
+   dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3))                !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
  end forall
 
 end subroutine constitutive_hooke_SandItsTangents
@@ -756,9 +752,6 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
    debug_level, &
    debug_constitutive, &
    debug_levelBasic
- use math, only: &
-   math_mul33x33, &
-   math_mul33x33
  use mesh, only: &
    theMesh
  use material, only: &
@@ -829,7 +822,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
  ho = material_homogenizationAt(el)
  tme = thermalMapping(ho)%p(ip,el)
 
- Mp  = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
+ Mp  = matmul(matmul(transpose(Fi),Fi),S)
 
  plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
 
@@ -897,8 +890,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
    debug_level, &
    debug_constitutive, &
    debug_levelBasic
- use math, only: &
-   math_mul33x33
  use material, only: &
    phasememberAt, &
    phase_plasticityInstance, &
@@ -931,7 +922,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
    i, &
    instance, of
 
- Mp  = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
+ Mp  = matmul(matmul(transpose(Fi),Fi),S)
 
  plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
 
@@ -966,8 +957,6 @@ end subroutine constitutive_collectDeltaState
 function constitutive_postResults(S, Fi, ipc, ip, el)
  use prec, only: &
    pReal
- use math, only: &
-  math_mul33x33
  use material, only: &
    phasememberAt, &
    phase_plasticityInstance, &
@@ -1035,7 +1024,7 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
 
  constitutive_postResults = 0.0_pReal
 
- Mp  = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
+ Mp  = matmul(matmul(transpose(Fi),Fi),S)
 
  ho = material_homogenizationAt(el)
  tme = thermalMapping(ho)%p(ip,el)
diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 55cc553ea..305043606 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -144,8 +144,7 @@ subroutine crystallite_init
  use math, only: &
    math_I3, &
    math_EulerToR, &
-   math_inv33, &
-   math_mul33x33
+   math_inv33
  use mesh, only: &
    theMesh, &
    mesh_element
@@ -353,7 +352,7 @@ subroutine crystallite_init
      crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e)
      crystallite_F0(1:3,1:3,c,i,e)  = math_I3
      crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phase(c,i,e))
-     crystallite_Fe(1:3,1:3,c,i,e)  = math_inv33(math_mul33x33(crystallite_Fi0(1:3,1:3,c,i,e), &
+     crystallite_Fe(1:3,1:3,c,i,e)  = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), &
                                                                crystallite_Fp0(1:3,1:3,c,i,e)))     ! assuming that euler angles are given in internal strain free configuration
      crystallite_Fp(1:3,1:3,c,i,e)  = crystallite_Fp0(1:3,1:3,c,i,e)
      crystallite_Fi(1:3,1:3,c,i,e)  = crystallite_Fi0(1:3,1:3,c,i,e)
@@ -430,8 +429,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
    IO_warning, &
    IO_error
  use math, only: &
-   math_inv33, &
-   math_mul33x33
+   math_inv33
  use mesh, only: &
    theMesh, &
    mesh_element
@@ -602,7 +600,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
            crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) &
                                            + crystallite_subStep(c,i,e) * (crystallite_partionedF (1:3,1:3,c,i,e) &
                                                                          - crystallite_partionedF0(1:3,1:3,c,i,e))
-           crystallite_Fe(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,c,i,e), &
+           crystallite_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), &
                                                                        crystallite_invFp(1:3,1:3,c,i,e)), &
                                                                        crystallite_invFi(1:3,1:3,c,i,e))
            crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e)
@@ -691,7 +689,6 @@ subroutine crystallite_stressTangent()
  use math, only: &
    math_inv33, &
    math_identity2nd, &
-   math_mul33x33, &
    math_3333to99, &
    math_99to3333, &
    math_I3, &
@@ -753,11 +750,11 @@ subroutine crystallite_stressTangent()
          lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
          do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
            lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
-                                 + crystallite_subdt(c,i,e)*math_mul33x33(invSubFi0,dLidFi(1:3,1:3,o,p))
+                                 + crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
            lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
                                  + crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e)
            rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
-                                 - crystallite_subdt(c,i,e)*math_mul33x33(invSubFi0,dLidS(1:3,1:3,o,p))
+                                 - crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
          enddo;enddo
          call math_invert2(temp_99,error,math_3333to99(lhs_3333))
          if (error) then
@@ -777,19 +774,19 @@ subroutine crystallite_stressTangent()
 
 !--------------------------------------------------------------------------------------------------
 ! calculate dSdF
-       temp_33_1 = transpose(math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
+       temp_33_1 = transpose(matmul(crystallite_invFp(1:3,1:3,c,i,e), &
                                            crystallite_invFi(1:3,1:3,c,i,e)))
-       temp_33_2 = math_mul33x33(           crystallite_subF  (1:3,1:3,c,i,e), &
+       temp_33_2 = matmul(           crystallite_subF  (1:3,1:3,c,i,e), &
                                  math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)))
-       temp_33_3 = math_mul33x33(math_mul33x33(crystallite_subF  (1:3,1:3,c,i,e), &
+       temp_33_3 = matmul(matmul(crystallite_subF  (1:3,1:3,c,i,e), &
                                                crystallite_invFp (1:3,1:3,c,i,e)), &
                                     math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)))
 
        forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) 
-         rhs_3333(p,o,1:3,1:3)  = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33_1)
-         temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33_2,dLpdS(1:3,1:3,p,o)), &
+         rhs_3333(p,o,1:3,1:3)  = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
+         temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), &
                                                 crystallite_invFi(1:3,1:3,c,i,e)) &
-                                + math_mul33x33(temp_33_3,dLidS(1:3,1:3,p,o))
+                                + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
        end forall
        lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) &
                 + math_mul3333xx3333(dSdFi,dFidS)
@@ -809,20 +806,20 @@ subroutine crystallite_stressTangent()
        forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt)
          dFpinvdF(1:3,1:3,p,o) &
            = -crystallite_subdt(c,i,e) &
-           * math_mul33x33(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), &
-                           math_mul33x33(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e)))
+           * matmul(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), &
+                           matmul(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e)))
        end forall
 
 !--------------------------------------------------------------------------------------------------
 ! assemble dPdF
-       temp_33_1 = math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
-                                 math_mul33x33(crystallite_S(1:3,1:3,c,i,e), &
+       temp_33_1 = matmul(crystallite_invFp(1:3,1:3,c,i,e), &
+                                 matmul(crystallite_S(1:3,1:3,c,i,e), &
                                                transpose(crystallite_invFp(1:3,1:3,c,i,e))))
-       temp_33_2 = math_mul33x33(crystallite_S(1:3,1:3,c,i,e), &
+       temp_33_2 = matmul(crystallite_S(1:3,1:3,c,i,e), &
                                  transpose(crystallite_invFp(1:3,1:3,c,i,e)))
-       temp_33_3 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
+       temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e), &
                                  crystallite_invFp(1:3,1:3,c,i,e))
-       temp_33_4 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
+       temp_33_4 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
                                                crystallite_invFp(1:3,1:3,c,i,e)), &
                                               crystallite_S(1:3,1:3,c,i,e))
 
@@ -832,9 +829,9 @@ subroutine crystallite_stressTangent()
        enddo
        forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt)
          crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
-           math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33_2) + &
-           math_mul33x33(math_mul33x33(temp_33_3,dSdF(1:3,1:3,p,o)),transpose(crystallite_invFp(1:3,1:3,c,i,e))) + &
-           math_mul33x33(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
+           matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33_2) + &
+           matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)),transpose(crystallite_invFp(1:3,1:3,c,i,e))) + &
+           matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
        end forall
 
    enddo; enddo
@@ -895,7 +892,6 @@ end subroutine crystallite_orientations
 !--------------------------------------------------------------------------------------------------
 function crystallite_push33ToRef(ipc,ip,el, tensor33)
  use math, only: &
-  math_mul33x33, &
   math_inv33, &
   math_EulerToR
  use material, only: &
@@ -910,9 +906,9 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33)
    ip, &                      ! integration point index
    ipc                         ! grain index
 
- T = math_mul33x33(math_EulerToR(material_EulerAngles(1:3,ipc,ip,el)), &
+ T = matmul(math_EulerToR(material_EulerAngles(1:3,ipc,ip,el)), &
                    transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el))))
- crystallite_push33ToRef = math_mul33x33(transpose(T),math_mul33x33(tensor33,T))
+ crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
 
 end function crystallite_push33ToRef
 
@@ -924,7 +920,6 @@ function crystallite_postResults(ipc, ip, el)
  use math, only: &
    math_qToEuler, &
    math_qToEulerAxisAngle, &
-   math_mul33x33, &
    math_det33, &
    math_I3, &
    inDeg
@@ -1093,11 +1088,7 @@ logical function integrateStress(&
  use constitutive, only: constitutive_LpAndItsTangents, &
                          constitutive_LiAndItsTangents, &
                          constitutive_SandItsTangents
- use math, only:         math_mul33x33, &
-#ifdef __PGI
-   norm2, &
-#endif
-                         math_mul33xx33, &
+ use math, only:         math_mul33xx33, &
                          math_mul3333xx3333, &
                          math_inv33, &
                          math_det33, &
@@ -1203,7 +1194,7 @@ logical function integrateStress(&
 #endif
    return
  endif failedInversionFp
- A = math_mul33x33(Fg_new,invFp_current)                                                            ! intermediate tensor needed later to calculate dFe_dLp
+ A = matmul(Fg_new,invFp_current)                                                            ! intermediate tensor needed later to calculate dFe_dLp
 
  invFi_current = math_inv33(crystallite_subFi0(1:3,1:3,ipc,ip,el))
  failedInversionFi: if (all(dEq0(invFi_current))) then
@@ -1235,7 +1226,7 @@ logical function integrateStress(&
      return
    endif LiLoopLimit
 
-   invFi_new = math_mul33x33(invFi_current,math_I3 - dt*Liguess)
+   invFi_new = matmul(invFi_current,math_I3 - dt*Liguess)
    Fi_new    = math_inv33(invFi_new)
    detInvFi  = math_det33(invFi_new)
 
@@ -1260,7 +1251,7 @@ logical function integrateStress(&
      !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
 
      B  = math_I3 - dt*Lpguess
-     Fe = math_mul33x33(math_mul33x33(A,B), invFi_new)
+     Fe = matmul(matmul(A,B), invFi_new)
      call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, &
                                        Fe, Fi_new, ipc, ip, el)                                     ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
 
@@ -1406,13 +1397,13 @@ logical function integrateStress(&
 
    !* calculate Jacobian for correction term
    if (mod(jacoCounterLi, iJacoLpresiduum) == 0_pInt) then
-     temp_33     = math_mul33x33(math_mul33x33(A,B),invFi_current)
+     temp_33     = matmul(matmul(A,B),invFi_current)
      forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt)
        dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33                                          ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
        dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current
      end forall
      forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) &
-       dFi_dLi(1:3,1:3,o,p) = math_mul33x33(math_mul33x33(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new)
+       dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new)
 
      dRLi_dLi  = math_identity2nd(9_pInt) &
                - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) + &
@@ -1449,7 +1440,7 @@ logical function integrateStress(&
  enddo LiLoop
 
  !* calculate new plastic and elastic deformation gradient
- invFp_new = math_mul33x33(invFp_current,B)
+ invFp_new = matmul(invFp_current,B)
  invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal)                               ! regularize
  Fp_new = math_inv33(invFp_new)
  failedInversionInvFp: if (all(dEq0(Fp_new))) then
@@ -1465,13 +1456,13 @@ logical function integrateStress(&
 #endif
    return
  endif failedInversionInvFp
- Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi_new)
+ Fe_new = matmul(matmul(Fg_new,invFp_new),invFi_new)
 
 !--------------------------------------------------------------------------------------------------
 ! stress integration was successful
  integrateStress = .true.
- crystallite_P    (1:3,1:3,ipc,ip,el) = math_mul33x33(math_mul33x33(Fg_new,invFp_new), &
-                                                      math_mul33x33(S,transpose(invFp_new)))
+ crystallite_P    (1:3,1:3,ipc,ip,el) = matmul(matmul(Fg_new,invFp_new), &
+                                                      matmul(S,transpose(invFp_new)))
  crystallite_S    (1:3,1:3,ipc,ip,el) = S
  crystallite_Lp   (1:3,1:3,ipc,ip,el) = Lpguess
  crystallite_Li   (1:3,1:3,ipc,ip,el) = Liguess
@@ -1489,9 +1480,9 @@ logical function integrateStress(&
    write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> P / MPa', &
               transpose(crystallite_P(1:3,1:3,ipc,ip,el))*1.0e-6_pReal
    write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Cauchy / MPa', &
-              math_mul33x33(crystallite_P(1:3,1:3,ipc,ip,el), transpose(Fg_new)) * 1.0e-6_pReal / math_det33(Fg_new)
+              matmul(crystallite_P(1:3,1:3,ipc,ip,el), transpose(Fg_new)) * 1.0e-6_pReal / math_det33(Fg_new)
    write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fe Lp Fe^-1', &
-              transpose(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,ipc,ip,el), math_inv33(Fe_new))))
+              transpose(matmul(Fe_new, matmul(crystallite_Lp(1:3,1:3,ipc,ip,el), math_inv33(Fe_new))))
    write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fp',transpose(crystallite_Fp(1:3,1:3,ipc,ip,el))
    write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fi',transpose(crystallite_Fi(1:3,1:3,ipc,ip,el))
  endif
diff --git a/src/grid_damage_spectral.f90 b/src/grid_damage_spectral.f90
index 1ce6e0c45..2019664e2 100644
--- a/src/grid_damage_spectral.f90
+++ b/src/grid_damage_spectral.f90
@@ -286,8 +286,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   use mesh, only: &
     grid, &
     grid3
-  use math, only: &
-    math_mul33x3
   use spectral_utilities, only: &
     scalarField_real, &
     vectorField_real, &
@@ -328,7 +326,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   cell = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     cell = cell + 1
-    vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
+    vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
                                                vectorField_real(1:3,i,j,k))
   enddo; enddo; enddo
   call utilities_FFTvectorForward
diff --git a/src/grid_mech_FEM.f90 b/src/grid_mech_FEM.f90
index 099d71d33..82273c8f1 100644
--- a/src/grid_mech_FEM.f90
+++ b/src/grid_mech_FEM.f90
@@ -316,7 +316,6 @@ end function grid_mech_FEM_solution
 !--------------------------------------------------------------------------------------------------
 subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
   use math, only: &
-    math_mul33x33 ,&
     math_rotate_backward33
   use numerics, only: &
     worldrank
@@ -402,7 +401,7 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
     ! calculate rate for aim
     if     (deformation_BC%myType=='l') then                                                          ! calculate F_aimDot from given L and current F
       F_aimDot = &
-      F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
+      F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
     elseif(deformation_BC%myType=='fdot') then                                                        ! F_aimDot is prescribed
       F_aimDot = &
       F_aimDot + deformation_BC%maskFloat * deformation_BC%values
diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90
index f17f2f8fd..99839e50f 100644
--- a/src/grid_mech_spectral_basic.f90
+++ b/src/grid_mech_spectral_basic.f90
@@ -285,7 +285,6 @@ end function grid_mech_spectral_basic_solution
 !--------------------------------------------------------------------------------------------------
 subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
   use math, only: &
-    math_mul33x33 ,&
     math_rotate_backward33
   use numerics, only: &
     worldrank
@@ -370,7 +369,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
     ! calculate rate for aim
     if     (deformation_BC%myType=='l') then                                                          ! calculate F_aimDot from given L and current F
       F_aimDot = &
-      F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
+      F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
     elseif(deformation_BC%myType=='fdot') then                                                        ! F_aimDot is prescribed
       F_aimDot = &
       F_aimDot + deformation_BC%maskFloat * deformation_BC%values
diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90
index 0a5501e98..aff4913b1 100644
--- a/src/grid_mech_spectral_polarisation.f90
+++ b/src/grid_mech_spectral_polarisation.f90
@@ -308,7 +308,6 @@ end function grid_mech_spectral_polarisation_solution
 !--------------------------------------------------------------------------------------------------
 subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
  use math, only: &
-   math_mul33x33, &
    math_mul3333xx33, &
    math_rotate_backward33
  use numerics, only: &
@@ -402,7 +401,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
     ! calculate rate for aim
     if     (deformation_BC%myType=='l') then                                                          ! calculate F_aimDot from given L and current F
       F_aimDot = &
-      F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
+      F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
     elseif(deformation_BC%myType=='fdot') then                                                        ! F_aimDot is prescribed
       F_aimDot = &
       F_aimDot + deformation_BC%maskFloat * deformation_BC%values
@@ -435,9 +434,9 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
   else
    do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
       F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
-      F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
+      F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, &
                                   math_mul3333xx33(C_scale,&
-                                                   math_mul33x33(transpose(F_lambda33),&
+                                                   matmul(transpose(F_lambda33),&
                                                                  F_lambda33)-math_I3))*0.5_pReal)&
                               + math_I3
       F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
@@ -528,8 +527,7 @@ subroutine formResidual(in, FandF_tau, &
    math_rotate_forward33, &
    math_rotate_backward33, &
    math_mul3333xx33, &
-   math_invSym3333, &
-   math_mul33x33
+   math_invSym3333
  use debug, only: &
    debug_level, &
    debug_spectral, &
@@ -605,7 +603,7 @@ subroutine formResidual(in, FandF_tau, &
  do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
    tensorField_real(1:3,1:3,i,j,k) = &
      polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
-     polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), &
+     polarAlpha*matmul(F(1:3,1:3,i,j,k), &
                         math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
  enddo; enddo; enddo
  
@@ -644,7 +642,7 @@ subroutine formResidual(in, FandF_tau, &
    e = e + 1
    residual_F(1:3,1:3,i,j,k) = &
      math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
-                      residual_F(1:3,1:3,i,j,k) - math_mul33x33(F(1:3,1:3,i,j,k), &
+                      residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
                       math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
                       + residual_F_tau(1:3,1:3,i,j,k)
  enddo; enddo; enddo
diff --git a/src/grid_thermal_spectral.f90 b/src/grid_thermal_spectral.f90
index 69e23c86b..adaf0d429 100644
--- a/src/grid_thermal_spectral.f90
+++ b/src/grid_thermal_spectral.f90
@@ -295,8 +295,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   use mesh, only: &
     grid, &
     grid3
-  use math, only: &
-    math_mul33x3
   use spectral_utilities, only: &
     scalarField_real, &
     vectorField_real, &
@@ -338,7 +336,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   cell = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     cell = cell + 1
-    vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, &
+    vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity33(1,cell) - D_ref, &
                                                vectorField_real(1:3,i,j,k))
   enddo; enddo; enddo
   call utilities_FFTvectorForward
diff --git a/src/homogenization_RGC.f90 b/src/homogenization_RGC.f90
index 1a170c66c..6a513193b 100644
--- a/src/homogenization_RGC.f90
+++ b/src/homogenization_RGC.f90
@@ -932,8 +932,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  !--------------------------------------------------------------------------------------------------
  function surfaceCorrection(avgF,instance,of)
   use math, only: &
-    math_invert33, &
-    math_mul33x33
+    math_invert33
   
   implicit none
   real(pReal), dimension(3)               :: surfaceCorrection
@@ -947,7 +946,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
   integer(pInt) :: i,j,iBase
   logical     ::  error
  
-  call math_invert33(math_mul33x33(transpose(avgF),avgF),invC,detF,error)
+  call math_invert33(matmul(transpose(avgF),avgF),invC,detF,error)
  
   surfaceCorrection = 0.0_pReal
   do iBase = 1_pInt,3_pInt
@@ -1139,8 +1138,6 @@ end function relaxationVector
 !> @brief identify the normal of an interface 
 !--------------------------------------------------------------------------------------------------
 pure function interfaceNormal(intFace,instance,of)
- use math, only: &
-   math_mul33x3
  
  implicit none
  real(pReal),   dimension (3)  ::            interfaceNormal
@@ -1156,7 +1153,7 @@ pure function interfaceNormal(intFace,instance,of)
  nPos = abs(intFace(1))                                                                             ! identify the position of the interface in global state array
  interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal)                                     ! get the normal vector w.r.t. cluster axis
 
- interfaceNormal = math_mul33x3(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal)   ! map the normal vector into sample coordinate system (basis)
+ interfaceNormal = matmul(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal)   ! map the normal vector into sample coordinate system (basis)
 
 end function interfaceNormal
 
diff --git a/src/lattice.f90 b/src/lattice.f90
index d3d2b3ce5..6c5a709e4 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -655,7 +655,6 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
  use prec, only: &
   tol_math_check
  use math, only: &
-   math_mul33x33, &
    math_sym3333to66, &
    math_Voigt66to3333, &
    math_cross
@@ -1141,8 +1140,7 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
    math_axisAngleToR, &
    math_sym3333to66, &
    math_66toSym3333, &
-   math_rotate_forward3333, &
-   math_mul33x33
+   math_rotate_forward3333
 
  implicit none
  integer,     dimension(:),             intent(in) :: Ntrans                                        !< number of active twin systems per family
@@ -1210,7 +1208,6 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
    INRAD, &
    math_outer, &
    math_cross, &
-   math_mul33x3, &
    math_axisAngleToR
  implicit none
  integer,     dimension(:),                intent(in) :: Nslip                                      !< number of active slip systems per family
@@ -1232,7 +1229,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
  do i = 1,sum(Nslip)
    direction = coordinateSystem(1:3,1,i)
    normal    = coordinateSystem(1:3,2,i)
-   np = math_mul33x3(math_axisAngleToR(direction,60.0_pReal*INRAD), normal)
+   np = matmul(math_axisAngleToR(direction,60.0_pReal*INRAD), normal)
    if (size(nonSchmidCoefficients)>0) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
      + nonSchmidCoefficients(1) * math_outer(direction, np)
    if (size(nonSchmidCoefficients)>1) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
@@ -2401,8 +2398,6 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
  use math, only: &
    math_cross, &
    math_outer, &
-   math_mul33x33, &
-   math_mul33x3, &
    math_axisAngleToR, &
    INRAD, &
    MATH_I3
@@ -2508,8 +2503,8 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
      U = (a_bcc/a_fcc)*math_outer(x,x) &
        + (a_bcc/a_fcc)*math_outer(y,y) * sqrt(2.0_pReal) &
        + (a_bcc/a_fcc)*math_outer(z,z) * sqrt(2.0_pReal)
-     Q(1:3,1:3,i) = math_mul33x33(R,B)
-     S(1:3,1:3,i) = math_mul33x33(R,U) - MATH_I3
+     Q(1:3,1:3,i) = matmul(R,B)
+     S(1:3,1:3,i) = matmul(R,U) - MATH_I3
    enddo
  elseif (cOverA > 0.0_pReal .and. dEq0(a_bcc)) then                                                 ! fcc -> hex transformation
    ss      = MATH_I3
@@ -2525,7 +2520,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
      Q(1:3,1,i) = x
      Q(1:3,2,i) = y
      Q(1:3,3,i) = z
-     S(1:3,1:3,i) = math_mul33x33(Q(1:3,1:3,i), math_mul33x33(math_mul33x33(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3  ! ToDo: This is of interest for the Schmid matrix only
+     S(1:3,1:3,i) = matmul(Q(1:3,1:3,i), matmul(matmul(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3  ! ToDo: This is of interest for the Schmid matrix only
    enddo
  else
    call IO_error(0) !ToDo: define error
diff --git a/src/math.f90 b/src/math.f90
index 8c0020799..5ee50cbfc 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -277,7 +277,7 @@ subroutine math_check
  ! +++ check rotation sense of q and R +++
  v = halton([2,8,5])     ! random vector
  R = math_qToR(q)
- if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then
+ if (any(abs(matmul(R,v) - math_qRot(q,v)) > tol_math_check)) then
    write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v'
    call IO_error(401,ext_msg=error_msg)
  endif
@@ -700,7 +700,7 @@ pure function math_exp33(A,n)
 
  do i = 1, merge(n,5,present(n))
    invFac = invFac/real(i,pReal)                                                                    ! invfac = 1/i!
-   B = math_mul33x33(B,A)
+   B = matmul(B,A)
    math_exp33 = math_exp33 + invFac*B                                                               ! exp = SUM (A^i)/i!
  enddo
 
@@ -1754,7 +1754,7 @@ real(pReal) pure function math_EulerMisorientation(EulerA,EulerB)
  real(pReal), dimension(3), intent(in) :: EulerA,EulerB
  real(pReal) :: cosTheta
 
- cosTheta = (math_trace33(math_mul33x33(math_EulerToR(EulerB), &
+ cosTheta = (math_trace33(matmul(math_EulerToR(EulerB), &
                               transpose(math_EulerToR(EulerA)))) - 1.0_pReal) * 0.5_pReal
 
  math_EulerMisorientation = acos(math_clip(cosTheta,-1.0_pReal,1.0_pReal))
@@ -1807,7 +1807,7 @@ function math_sampleGaussOri(center,FWHM)
      angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R))
      if (rnd(4) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) exit                       ! rejection sampling (Gaussian)
    enddo GaussConvolution
-   math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center)))
+   math_sampleGaussOri = math_RtoEuler(matmul(R,math_EulerToR(center)))
  endif
 
 end function math_sampleGaussOri
@@ -1842,7 +1842,7 @@ function math_sampleFiberOri(alpha,beta,FWHM)
  R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS)))             !< rotation to align fiber axis in crystal and sample system
 
  rnd = halton([7,10,3])
- R = math_mul33x33(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI))                              !< additional rotation (0..360deg) perpendicular to fiber axis
+ R = matmul(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI))                              !< additional rotation (0..360deg) perpendicular to fiber axis
 
  if (FWHM > 0.1_pReal*INRAD) then
    reducedTo2D: do i=1,3
@@ -1863,7 +1863,7 @@ function math_sampleFiberOri(alpha,beta,FWHM)
      u(j)    = fInS(j)
 
      rejectionSampling: if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) then
-       R = math_mul33x33(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle))                 ! tilt around direction of smallest component
+       R = matmul(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle))                 ! tilt around direction of smallest component
        exit
      endif rejectionSampling
      rnd = halton([7,10,3])
@@ -2079,23 +2079,23 @@ pure function math_eigenvectorBasisSym33(m)
    N(1:3,1:3,2) = m-values(2)*math_I3
    N(1:3,1:3,3) = m-values(3)*math_I3
    twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then 
-     EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
+     EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
                                                ((values(3)-values(1))*(values(3)-values(2)))
      EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
    elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
-     EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
+     EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
                                                ((values(1)-values(2))*(values(1)-values(3)))
      EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
    elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues 
-     EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
+     EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
                                                ((values(2)-values(1))*(values(2)-values(3)))
      EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
    else twoSimilarEigenvalues
-     EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
+     EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
                                                ((values(1)-values(2))*(values(1)-values(3)))
-     EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
+     EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
                                                ((values(2)-values(1))*(values(2)-values(3)))
-     EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
+     EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
                                                ((values(3)-values(1))*(values(3)-values(2)))
    endif twoSimilarEigenvalues
  endif threeSimilarEigenvalues
@@ -2144,23 +2144,23 @@ pure function math_eigenvectorBasisSym33_log(m)
    N(1:3,1:3,2) = m-values(2)*math_I3
    N(1:3,1:3,3) = m-values(3)*math_I3
    twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then 
-     EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
+     EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
                                                ((values(3)-values(1))*(values(3)-values(2)))
      EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
    elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
-     EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
+     EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
                                                ((values(1)-values(2))*(values(1)-values(3)))
      EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
    elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues 
-     EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
+     EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
                                                ((values(2)-values(1))*(values(2)-values(3)))
      EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
    else twoSimilarEigenvalues
-     EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
+     EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
                                                ((values(1)-values(2))*(values(1)-values(3)))
-     EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
+     EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
                                                ((values(2)-values(1))*(values(2)-values(3)))
-     EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
+     EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
                                                ((values(3)-values(1))*(values(3)-values(2)))
    endif twoSimilarEigenvalues
  endif threeSimilarEigenvalues
@@ -2186,14 +2186,14 @@ function math_rotationalPart33(m)
  real(pReal), dimension(3,3) :: math_rotationalPart33
  real(pReal), dimension(3,3) :: U , Uinv
 
- U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m))
+ U = math_eigenvectorBasisSym33(matmul(transpose(m),m))
  Uinv = math_inv33(U)
 
  inversionFailed: if (all(dEq0(Uinv))) then
    math_rotationalPart33 = math_I3
    call IO_warning(650)
  else inversionFailed
-   math_rotationalPart33 = math_mul33x33(m,Uinv)
+   math_rotationalPart33 = matmul(m,Uinv)
  endif inversionFailed
 
 end function math_rotationalPart33
@@ -2586,7 +2586,7 @@ 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
 
- math_rotate_forward33 = math_mul33x33(rot_tensor,math_mul33x33(tensor,transpose(rot_tensor)))
+ math_rotate_forward33 = matmul(rot_tensor,matmul(tensor,transpose(rot_tensor)))
 
 end function math_rotate_forward33
 
@@ -2600,7 +2600,7 @@ pure function math_rotate_backward33(tensor,rot_tensor)
  real(pReal), dimension(3,3) ::  math_rotate_backward33
  real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
 
- math_rotate_backward33 = math_mul33x33(transpose(rot_tensor),math_mul33x33(tensor,rot_tensor))
+ math_rotate_backward33 = matmul(transpose(rot_tensor),matmul(tensor,rot_tensor))
 
 end function math_rotate_backward33
 
diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90
index 38abad1aa..95d8e5b27 100644
--- a/src/mesh_grid.f90
+++ b/src/mesh_grid.f90
@@ -561,8 +561,6 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
    debug_mesh, &
    debug_level, &
    debug_levelBasic
- use math, only: &
-   math_mul33x3
 
  implicit none
  real(pReal), intent(in), dimension(:,:,:,:) :: &
@@ -624,7 +622,7 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
          lookup = me-diag+shift*iRes
          wrappedCentres(1:3,i+1_pInt,        j+1_pInt,        k+1_pInt) = &
                 centres(1:3,lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt) &
-                - math_mul33x3(Favg, real(shift,pReal)*gDim)
+                - matmul(Favg, real(shift,pReal)*gDim)
        endif
  enddo; enddo; enddo
 
diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90
index a0a996dd6..cb13265b4 100644
--- a/src/plastic_dislotwin.f90
+++ b/src/plastic_dislotwin.f90
@@ -651,8 +651,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
    math_eigenValuesVectorsSym, &
    math_outer, &
    math_symmetric33, &
-   math_mul33xx33, &
-   math_mul33x3
+   math_mul33xx33
  
  implicit none
  real(pReal), dimension(3,3),     intent(out) :: Lp
@@ -723,8 +722,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
    call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error)
 
    do i = 1,6
-     P_sb = 0.5_pReal * math_outer(math_mul33x3(eigVectors,sb_sComposition(1:3,i)),&
-                                               math_mul33x3(eigVectors,sb_mComposition(1:3,i)))
+     P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
+                                   matmul(eigVectors,sb_mComposition(1:3,i)))
      tau = math_mul33xx33(Mp,P_sb)
    
      significantShearBandStress: if (abs(tau) > tol_math_check) then
diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90
index f0b28d711..94f0fde04 100644
--- a/src/plastic_nonlocal.f90
+++ b/src/plastic_nonlocal.f90
@@ -838,8 +838,7 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
     IO_error
   use math, only: &
     PI, &
-    math_mul33x3, &
-    math_mul3x3, &
+    math_inner, &
     math_inv33
 #ifdef DEBUG
   use debug, only: &
@@ -1004,10 +1003,10 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
             neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor(:,scr)),2)
   
             connection_latticeConf(1:3,n) = &
-              math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) &
+              matmul(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) &
                                   - mesh_ipCoordinates(1:3,ip,el))
-            normal_latticeConf = math_mul33x3(transpose(invFp), mesh_ipAreaNormal(1:3,n,ip,el))
-            if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) &        ! neighboring connection points in opposite direction to face normal: must be periodic image
+            normal_latticeConf = matmul(transpose(invFp), mesh_ipAreaNormal(1:3,n,ip,el))
+            if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) &        ! neighboring connection points in opposite direction to face normal: must be periodic image
               connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el)/mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
         else
           ! local neighbor or different lattice structure or different constitution instance -> use central values instead
@@ -1047,7 +1046,7 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
         invConnections = math_inv33(connections)
         if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error')
   
-        rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), math_mul33x3(invConnections,rhoExcessDifferences))
+        rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences))
       enddo
         
         ! ... plus gradient from deads ...
@@ -1528,13 +1527,8 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
     debug_e
 #endif
   use math, only: &
-#ifdef __PGI
-    norm2, &
-#endif
-    math_mul3x3, &
-    math_mul33x3, &
+    math_inner, &
     math_mul33xx33, &
-    math_mul33x33, &
     math_inv33, &
     math_det33, &
     PI
@@ -1756,7 +1750,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
     m(1:3,1:ns,4) =  prm%slip_transverse
     
     my_Fe = Fe(1:3,1:3,1,ip,el)
-    my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,1,ip,el))
+    my_F = matmul(my_Fe, Fp(1:3,1:3,1,ip,el))
     
     neighbors: do n = 1,theMesh%elem%nIPneighbors
   
@@ -1774,7 +1768,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
       if (neighbor_n > 0) then                                                                      ! if neighbor exists, average deformation gradient
         neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el))
         neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el)
-        neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el))
+        neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el))
         Favg = 0.5_pReal * (my_F + neighbor_F)
       else                                                                                          ! if no neighbor, take my value as average
         Favg = my_F
@@ -1809,9 +1803,9 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
         where (neighbor_rhoSgl * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN &
           .or. neighbor_rhoSgl < prm%significantRho) &
           neighbor_rhoSgl = 0.0_pReal
-        normal_neighbor2me_defConf = math_det33(Favg) * math_mul33x3(math_inv33(transpose(Favg)), &
+        normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
                                      mesh_ipAreaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el))     ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
-        normal_neighbor2me = math_mul33x3(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
+        normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
                            / math_det33(neighbor_Fe)                                                ! interface normal in the lattice configuration of my neighbor
         area = mesh_ipArea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me)
         normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me)                         ! normalize the surface normal to unit length
@@ -1819,10 +1813,10 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
           do t = 1,4
             c = (t + 1) / 2
             topp = t + mod(t,2) - mod(t+1,2)
-            if (neighbor_v(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal &         ! flux from my neighbor to me == entering flux for me
+            if (neighbor_v(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal &         ! flux from my neighbor to me == entering flux for me
                 .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then                                  ! ... only if no sign change in flux density  
               lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) &
-                         * math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area                       ! positive line length that wants to enter through this interface
+                         * math_inner(m(1:3,s,t), normal_neighbor2me) * area                       ! positive line length that wants to enter through this interface
               where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) &                                 ! positive compatibility...
                 rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) &
                                    + lineLength / mesh_ipVolume(ip,el) &                            ! ... transferring to equally signed mobile dislocation type
@@ -1856,23 +1850,23 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
         my_v = v
         
         normal_me2neighbor_defConf = math_det33(Favg) &
-                                   * math_mul33x3(math_inv33(transpose(Favg)), & 
+                                   * matmul(math_inv33(transpose(Favg)), & 
                                                              mesh_ipAreaNormal(1:3,n,ip,el))        ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!)
-        normal_me2neighbor = math_mul33x3(transpose(my_Fe), normal_me2neighbor_defConf) &
+        normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) &
                            / math_det33(my_Fe)                                                      ! interface normal in my lattice configuration
         area = mesh_ipArea(n,ip,el) * norm2(normal_me2neighbor)
         normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor)                         ! normalize the surface normal to unit length    
         do s = 1,ns
           do t = 1,4
             c = (t + 1) / 2
-            if (my_v(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then          ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
+            if (my_v(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then          ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
               if (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then                                    ! no sign change in flux density
                 transmissivity = sum(compatibility(c,1:ns,s,n,ip,el)**2.0_pReal)                    ! overall transmissivity from this slip system to my neighbor
               else                                                                                  ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
                 transmissivity = 0.0_pReal
               endif
               lineLength = my_rhoSgl(s,t) * my_v(s,t) &
-                         * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area                       ! positive line length of mobiles that wants to leave through this interface
+                         * math_inner(m(1:3,s,t), normal_me2neighbor) * area                       ! positive line length of mobiles that wants to leave through this interface
               rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el)                 ! subtract dislocation flux from current type
               rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
                                      + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) &
@@ -2017,7 +2011,7 @@ end subroutine plastic_nonlocal_dotState
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) 
   use math, only: &
-    math_mul3x3, &
+    math_inner, &
     math_qRot
   use rotations, only: &
     rotation
@@ -2134,13 +2128,13 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
       absoluteMisorientation = rot%asQuaternion()
       mySlipSystems: do s1 = 1,ns
         neighborSlipSystems: do s2 = 1,ns
-          my_compatibility(1,s2,s1,n) =  math_mul3x3(prm%slip_normal(1:3,s1), &
+          my_compatibility(1,s2,s1,n) =  math_inner(prm%slip_normal(1:3,s1), &
           math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2))) &
-                                  * abs(math_mul3x3(prm%slip_direction(1:3,s1), &
+                                  * abs(math_inner(prm%slip_direction(1:3,s1), &
                                   math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2))))
-          my_compatibility(2,s2,s1,n) = abs(math_mul3x3(prm%slip_normal(1:3,s1), &
+          my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
           math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2)))) &
-                                  * abs(math_mul3x3(prm%slip_direction(1:3,s1), &
+                                  * abs(math_inner(prm%slip_direction(1:3,s1), &
                                   math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2))))
         enddo neighborSlipSystems
         
diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90
index d6ee268a3..2d13277c7 100644
--- a/src/source_damage_isoBrittle.f90
+++ b/src/source_damage_isoBrittle.f90
@@ -174,8 +174,6 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
    sourceState
  use math, only : &
    math_sym33to6, &
-   math_mul33x33, &
-   math_mul66x6, &
    math_I3
 
  implicit none
@@ -200,9 +198,10 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
  sourceOffset = source_damage_isoBrittle_offset(phase)
 
                                    
- strain = 0.5_pReal*math_sym33to6(math_mul33x33(transpose(Fe),Fe)-math_I3)
+ strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
 
- strainenergy = 2.0_pReal*sum(strain*math_mul66x6(C,strain))/param(instance)%critStrainEnergy
+ strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/param(instance)%critStrainEnergy
+ ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/param(instance)%critStrainEnergy
  
  if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
    sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90
index 26c1cf9fc..dd79e67e2 100644
--- a/src/spectral_utilities.f90
+++ b/src/spectral_utilities.f90
@@ -610,7 +610,6 @@ end subroutine utilities_fourierGammaConvolution
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT)
   use math, only: &
-    math_mul33x3, &
     PI
   use mesh, only: &
     grid, &
@@ -1158,8 +1157,6 @@ subroutine utilities_updateIPcoords(F)
     cNeq
   use IO, only: &
     IO_error
-  use math, only: &
-    math_mul33x3
   use mesh, only: &
     grid, &
     grid3, &
@@ -1200,12 +1197,12 @@ subroutine utilities_updateIPcoords(F)
   if (grid3Offset == 0) offset_coords = vectorField_real(1:3,1,1,1)
   call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
   if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords')
-  offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords
+  offset_coords = matmul(Favg,step/2.0_pReal) - offset_coords
   m = 1
   do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1)
     mesh_ipCoordinates(1:3,1,m) = vectorField_real(1:3,i,j,k) &
                                 + offset_coords &
-                                + math_mul33x3(Favg,step*real([i,j,k+grid3Offset]-1,pReal))
+                                + matmul(Favg,step*real([i,j,k+grid3Offset]-1,pReal))
     m = m+1
   enddo; enddo; enddo