diff --git a/Makefile b/Makefile
index 5f3459f1a..ef02313c2 100644
--- a/Makefile
+++ b/Makefile
@@ -65,7 +65,7 @@ LMKL = -L$(MKLROOT)/lib/intel64 -qopt-matmul
 
 MODULES         = $(SWIFTEST_MODULES) $(USER_MODULES) 
 
-.PHONY : all mod lib libdir fast drivers bin clean force 
+.PHONY : all mod lib fast drivers bin clean force 
 
 % : %.f90 force
 	$(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) $< -o $@ \
@@ -118,11 +118,6 @@ lib:
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
 	  ln -s $(SWIFTEST_HOME)/Makefile .; \
 	  make libdir
-	cd $(SWIFTEST_HOME)/src/operators; \
-	  rm -f Makefile.Defines Makefile; \
-	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
-	  ln -s $(SWIFTEST_HOME)/Makefile .; \
-	  make libdir
 	cd $(SWIFTEST_HOME)/src/setup; \
 	  rm -f Makefile.Defines Makefile; \
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
@@ -165,6 +160,11 @@ lib:
 	  make libdir
 
 fast:
+	cd $(SWIFTEST_HOME)/src/operators; \
+	  rm -f Makefile.Defines Makefile; \
+	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
+	  ln -s $(SWIFTEST_HOME)/Makefile .; \
+	  make fastdir
 	cd $(SWIFTEST_HOME)/src/fraggle; \
 	  rm -f Makefile.Defines Makefile; \
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
diff --git a/Makefile.Defines b/Makefile.Defines
index ea2772b04..838038524 100644
--- a/Makefile.Defines
+++ b/Makefile.Defines
@@ -69,8 +69,8 @@ GPRODUCTION = -O2 -ffree-line-length-none $(GPAR)
 
 #FFLAGS  = $(IDEBUG)  #$(SIMDVEC) $(PAR)
 #FFASTFLAGS  = $(IDEBUG) #$(SIMDVEC) $(PAR)
-FFLAGS = $(IPRODUCTION) $(STRICTREAL) #$(ADVIXE_FLAGS)
-FFASTFLAGS = $(IPRODUCTION) -fp-model fast #$(ADVIXE_FLAGS)
+FFLAGS = $(IPRODUCTION) $(STRICTREAL) $(OPTREPORT) #$(ADVIXE_FLAGS)
+FFASTFLAGS = $(IPRODUCTION) -fp-model fast $(OPTREPORT) #$(ADVIXE_FLAGS)
 FORTRAN = ifort
 AR     = xiar
 
diff --git a/src/modules/swiftest_operators.f90 b/src/modules/swiftest_operators.f90
index 2c982f09c..13cb57839 100644
--- a/src/modules/swiftest_operators.f90
+++ b/src/modules/swiftest_operators.f90
@@ -15,42 +15,49 @@ module swiftest_operators
 
    interface operator(.cross.)
       module pure function operator_cross_sp(A, B) result(C)
+         !$omp declare simd(operator_cross_sp)
          implicit none
          real(SP), dimension(:), intent(in) :: A, B
          real(SP), dimension(3) :: C
       end function operator_cross_sp
 
       module pure function operator_cross_dp(A, B) result(C)
+         !$omp declare simd(operator_cross_dp)
          implicit none
          real(DP), dimension(:), intent(in) :: A, B
          real(DP), dimension(3) :: C
       end function operator_cross_dp
 
       module pure function operator_cross_qp(A, B) result(C)
+         !$omp declare simd(operator_cross_qp)
          implicit none
          real(QP), dimension(:), intent(in) :: A, B
          real(QP), dimension(3) :: C
       end function operator_cross_qp
 
       module pure function operator_cross_i1b(A, B) result(C)
+         !$omp declare simd(operator_cross_i1b)
          implicit none
          integer(I1B), dimension(:), intent(in) :: A, B
          integer(I1B), dimension(3) :: C
       end function operator_cross_i1b
 
       module pure function operator_cross_i2b(A, B) result(C)
+         !$omp declare simd(operator_cross_i2b)
          implicit none
          integer(I2B), dimension(:), intent(in) :: A, B
          integer(I2B), dimension(3) :: C
       end function operator_cross_i2b
 
       module pure function operator_cross_i4b(A, B) result(C)
+         !$omp declare simd(operator_cross_i4b)
          implicit none
          integer(I4B), dimension(:), intent(in) :: A, B
          integer(I4B), dimension(3) :: C
       end function operator_cross_i4b
 
       module pure function operator_cross_i8b(A, B) result(C)
+         !$omp declare simd(operator_cross_i8b)
          implicit none
          integer(I8B), dimension(:), intent(in) :: A, B
          integer(I8B), dimension(3) :: C
@@ -105,18 +112,21 @@ end function operator_cross_el_i8b
 
    interface operator(.mag.)
       module pure function operator_mag_sp(A) result(B)
+         !$omp declare simd(operator_mag_sp)
          implicit none
          real(SP), dimension(:), intent(in) :: A
          real(SP)                           :: B
       end function operator_mag_sp
 
       module pure function operator_mag_dp(A) result(B)
+         !$omp declare simd(operator_mag_dp)
          implicit none
          real(DP), dimension(:), intent(in) :: A
          real(DP)                           :: B
       end function operator_mag_dp
 
       module pure function operator_mag_qp(A) result(B)
+         !$omp declare simd(operator_mag_qp)
          implicit none
          real(QP), dimension(:), intent(in) :: A
          real(QP)                           :: B
diff --git a/src/operators/operator_cross.f90 b/src/operators/operator_cross.f90
index 736dc2696..99de730fe 100644
--- a/src/operators/operator_cross.f90
+++ b/src/operators/operator_cross.f90
@@ -8,6 +8,7 @@
 contains
 
    module pure function operator_cross_sp(A, B) result(C)
+      !$omp declare simd(operator_cross_sp)
       implicit none
       real(SP), dimension(:), intent(in) :: A, B
       real(SP), dimension(3) :: C
@@ -18,6 +19,7 @@ module pure function operator_cross_sp(A, B) result(C)
    end function operator_cross_sp
 
    module pure function operator_cross_dp(A, B) result(C)
+      !$omp declare simd(operator_cross_dp)
       implicit none
       real(DP), dimension(:), intent(in) :: A, B
       real(DP), dimension(3) :: C
@@ -28,6 +30,7 @@ module pure function operator_cross_dp(A, B) result(C)
    end function operator_cross_dp
 
    module pure function operator_cross_qp(A, B) result(C)
+      !$omp declare simd(operator_cross_qp)
       implicit none
       real(QP), dimension(:), intent(in) :: A, B
       real(QP), dimension(3) :: C
@@ -38,6 +41,7 @@ module pure function operator_cross_qp(A, B) result(C)
    end function operator_cross_qp
 
    module pure function operator_cross_i1b(A, B) result(C)
+      !$omp declare simd(operator_cross_i1b)
       implicit none
       integer(I1B), dimension(:), intent(in) :: A, B
       integer(I1B), dimension(3) :: C
@@ -48,6 +52,7 @@ module pure function operator_cross_i1b(A, B) result(C)
    end function operator_cross_i1b
 
    module pure function operator_cross_i2b(A, B) result(C)
+      !$omp declare simd(operator_cross_i2b)
       implicit none
       integer(I2B), dimension(:), intent(in) :: A, B
       integer(I2B), dimension(3) :: C
@@ -58,6 +63,7 @@ module pure function operator_cross_i2b(A, B) result(C)
    end function operator_cross_i2b
 
    module pure function operator_cross_i4b(A, B) result(C)
+      !$omp declare simd(operator_cross_i4b)
       implicit none
       integer(I4B), dimension(:), intent(in) :: A, B
       integer(I4B), dimension(3) :: C
@@ -68,6 +74,7 @@ module pure function operator_cross_i4b(A, B) result(C)
    end function operator_cross_i4b
 
    module pure function operator_cross_i8b(A, B) result(C)     
+      !$omp declare simd(operator_cross_i8b)
       implicit none
       integer(I8B), dimension(:), intent(in) :: A, B
       integer(I8B), dimension(3) :: C
@@ -86,9 +93,7 @@ module pure function operator_cross_el_sp(A, B) result(C)
       if (allocated(C)) deallocate(C)
       allocate(C, mold = A)
       do concurrent (i = 1:n) 
-         C(1, i) = A(2, i) * B(3, i) - A(3, i) * B(2, i)
-         C(2, i) = A(3, i) * B(1, i) - A(1, i) * B(3, i)
-         C(3, i) = A(1, i) * B(2, i) - A(2, i) * B(1, i)
+         C(:,i) = operator_cross_sp(A(:,i), B(:,i))
       end do
       return
    end function operator_cross_el_sp
@@ -102,9 +107,7 @@ module pure function operator_cross_el_dp(A, B) result(C)
       if (allocated(C)) deallocate(C)
       allocate(C, mold = A)
       do concurrent (i = 1:n) 
-         C(1, i) = A(2, i) * B(3, i) - A(3, i) * B(2, i)
-         C(2, i) = A(3, i) * B(1, i) - A(1, i) * B(3, i)
-         C(3, i) = A(1, i) * B(2, i) - A(2, i) * B(1, i)
+         C(:,i) = operator_cross_dp(A(:,i), B(:,i))
       end do
       return
    end function operator_cross_el_dp
@@ -118,9 +121,7 @@ module pure function operator_cross_el_qp(A, B) result(C)
       if (allocated(C)) deallocate(C)
       allocate(C, mold = A)
       do concurrent (i = 1:n) 
-         C(1, i) = A(2, i) * B(3, i) - A(3, i) * B(2, i)
-         C(2, i) = A(3, i) * B(1, i) - A(1, i) * B(3, i)
-         C(3, i) = A(1, i) * B(2, i) - A(2, i) * B(1, i)
+         C(:,i) = operator_cross_qp(A(:,i), B(:,i))
       end do
       return
    end function operator_cross_el_qp
@@ -134,9 +135,7 @@ module pure function operator_cross_el_i1b(A, B) result(C)
       if (allocated(C)) deallocate(C)
       allocate(C, mold = A)
       do concurrent (i = 1:n) 
-         C(1, i) = A(2, i) * B(3, i) - A(3, i) * B(2, i)
-         C(2, i) = A(3, i) * B(1, i) - A(1, i) * B(3, i)
-         C(3, i) = A(1, i) * B(2, i) - A(2, i) * B(1, i)
+         C(:,i) = operator_cross_i1b(A(:,i), B(:,i))
       end do
       return
    end function operator_cross_el_i1b
@@ -150,9 +149,7 @@ module pure function operator_cross_el_i2b(A, B) result(C)
       if (allocated(C)) deallocate(C)
       allocate(C, mold = A)
       do concurrent (i = 1:n) 
-         C(1, i) = A(2, i) * B(3, i) - A(3, i) * B(2, i)
-         C(2, i) = A(3, i) * B(1, i) - A(1, i) * B(3, i)
-         C(3, i) = A(1, i) * B(2, i) - A(2, i) * B(1, i)
+         C(:,i) = operator_cross_i2b(A(:,i), B(:,i))
       end do
       return
    end function operator_cross_el_i2b
@@ -166,9 +163,7 @@ module pure function operator_cross_el_i4b(A, B) result(C)
       if (allocated(C)) deallocate(C)
       allocate(C, mold = A)
       do concurrent (i = 1:n) 
-         C(1, i) = A(2, i) * B(3, i) - A(3, i) * B(2, i)
-         C(2, i) = A(3, i) * B(1, i) - A(1, i) * B(3, i)
-         C(3, i) = A(1, i) * B(2, i) - A(2, i) * B(1, i)
+         C(:,i) = operator_cross_i4b(A(:,i), B(:,i))
       end do
       return
    end function operator_cross_el_i4b
@@ -182,9 +177,7 @@ module pure function operator_cross_el_i8b(A, B) result(C)
       if (allocated(C)) deallocate(C)
       allocate(C, mold = A)
       do concurrent (i = 1:n) 
-         C(1, i) = A(2, i) * B(3, i) - A(3, i) * B(2, i)
-         C(2, i) = A(3, i) * B(1, i) - A(1, i) * B(3, i)
-         C(3, i) = A(1, i) * B(2, i) - A(2, i) * B(1, i)
+         C(:,i) = operator_cross_i8b(A(:,i), B(:,i))
       end do
       return
    end function operator_cross_el_i8b
diff --git a/src/operators/operator_mag.f90 b/src/operators/operator_mag.f90
index 5a054d5ce..f74555775 100644
--- a/src/operators/operator_mag.f90
+++ b/src/operators/operator_mag.f90
@@ -7,6 +7,7 @@
    contains
 
    module pure function operator_mag_sp(A) result(B)
+      !$omp declare simd(operator_mag_sp)
       implicit none
       real(SP), dimension(:), intent(in) :: A
       real(SP)                           :: B
@@ -15,6 +16,7 @@ module pure function operator_mag_sp(A) result(B)
    end function operator_mag_sp
 
    module pure function operator_mag_dp(A) result(B)
+      !$omp declare simd(operator_mag_dp)
       implicit none
       real(DP), dimension(:), intent(in) :: A
       real(DP)                           :: B