diff --git a/R/lav_mvnorm.R b/R/lav_mvnorm.R
index dd820e00..9d15d7b3 100644
--- a/R/lav_mvnorm.R
+++ b/R/lav_mvnorm.R
@@ -743,11 +743,19 @@ lav_mvnorm_information_expected <- function(Y = NULL, # unused!
     )
   }
 
-  if(correlation) {
-    I22 <- 0.5 * lav_matrix_duplication_cor_pre_post(Sigma.inv %x%
-                                                     Sigma.inv)
+  if (lav_use_lavaanC()) {
+    if (correlation) {
+      I22 <- lavaanC::m_kronecker_dup_cor_pre_post(Sigma.inv, 
+                                                   multiplicator = 0.5)
+    } else {
+      I22 <- lavaanC::m_kronecker_dup_pre_post(Sigma.inv, multiplicator = 0.5)
+    }
+  } else {
+    if (correlation) {
+      I22 <- 0.5 * lav_matrix_duplication_cor_pre_post(Sigma.inv %x% Sigma.inv)
     } else {
       I22 <- 0.5 * lav_matrix_duplication_pre_post(Sigma.inv %x% Sigma.inv)
+    }
   }
 
   # fixed.x?
@@ -835,7 +843,11 @@ lav_mvnorm_information_observed_samplestats <- function(
   }
 
   AAA <- Sigma.inv %*% (2 * W.tilde - Sigma) %*% Sigma.inv
-  I22 <- (1 / 2) * lav_matrix_duplication_pre_post(Sigma.inv %x% AAA)
+  if (lav_use_lavaanC()) {
+    I22 <- lavaanC::m_kronecker_dup_pre_post(Sigma.inv, AAA, 0.5)
+  } else {
+    I22 <- (1 / 2) * lav_matrix_duplication_pre_post(Sigma.inv %x% AAA)
+  }
 
   if (meanstructure) {
     out <- rbind(
@@ -937,8 +949,13 @@ lav_mvnorm_inverted_information_expected <- function(Y = NULL, # unused!
     # take difference
     R <- Sigma - YbarX.aug
 
-    SS <- 2 * lav_matrix_duplication_ginv_pre_post(Sigma %x% Sigma)
-    RR <- 2 * lav_matrix_duplication_ginv_pre_post(R %x% R)
+    if (lav_use_lavaanC()) {
+      SS <- lavaanC::m_kronecker_dup_ginv_pre_post(Sigma, multiplicator = 2.0)
+      RR <- lavaanC::m_kronecker_dup_ginv_pre_post(R, multiplicator = 2.0)
+    } else {
+      SS <- 2 * lav_matrix_duplication_ginv_pre_post(Sigma %x% Sigma)
+      RR <- 2 * lav_matrix_duplication_ginv_pre_post(R %x% R)
+    }
     I22 <- SS - RR
 
     if (meanstructure) {
@@ -948,7 +965,11 @@ lav_mvnorm_inverted_information_expected <- function(Y = NULL, # unused!
       out <- I22
     }
   } else {
-    I22 <- 2 * lav_matrix_duplication_ginv_pre_post(Sigma %x% Sigma)
+    if (lav_use_lavaanC()) {
+      I22 <- lavaanC::m_kronecker_dup_ginv_pre_post(Sigma, multiplicator = 2.0)
+    } else {
+      I22 <- 2 * lav_matrix_duplication_ginv_pre_post(Sigma %x% Sigma)
+    }
     if (meanstructure) {
       I11 <- Sigma
       out <- lav_matrix_bdiag(I11, I22)
diff --git a/R/lav_mvnorm_cluster.R b/R/lav_mvnorm_cluster.R
index 6b40e48d..8ea8779c 100644
--- a/R/lav_mvnorm_cluster.R
+++ b/R/lav_mvnorm_cluster.R
@@ -1012,7 +1012,11 @@ lav_mvnorm_cluster_information_expected <- function(Lp = NULL,
     omega.j.inv <- solve(omega.j)
 
     I11.j <- omega.j.inv
-    I22.j <- 0.5 * lav_matrix_duplication_pre_post(omega.j.inv %x% omega.j.inv)
+    if (lav_use_lavaanC()) {
+      I22.j <- lavaanC::m_kronecker_dup_pre_post(omega.j.inv, multiplicator = 0.5)
+    } else {
+      I22.j <- 0.5 * lav_matrix_duplication_pre_post(omega.j.inv %x% omega.j.inv)
+    }
     I.j <- lav_matrix_bdiag(I11.j, I22.j)
     info.j <- t(Delta.j) %*% I.j %*% Delta.j
 
@@ -1028,7 +1032,11 @@ lav_mvnorm_cluster_information_expected <- function(Lp = NULL,
   Sigma.W.inv.tilde[ov.idx[[1]], ov.idx[[1]]] <- Sigma.W.inv
 
   I11.w <- Sigma.W.inv.tilde
-  I22.w <- 0.5 * lav_matrix_duplication_pre_post(Sigma.W.inv.tilde %x% Sigma.W.inv.tilde)
+  if (lav_use_lavaanC()) {
+    I22.W <- lavaanC::m_kronecker_dup_pre_post(Sigma.W.inv.tilde, multiplicator = 0.5)
+  } else {
+    I22.w <- 0.5 * lav_matrix_duplication_pre_post(Sigma.W.inv.tilde %x% Sigma.W.inv.tilde)
+  }
   I.w <- lav_matrix_bdiag(I11.w, I22.w)
   information.w <- (nobs - nclusters) *
     (t(Delta.W.tilde) %*% I.w %*% Delta.W.tilde)
@@ -1160,7 +1168,11 @@ lav_mvnorm_cluster_information_expected_delta <- function(Lp = NULL,
     omega.j.inv <- solve(omega.j)
 
     I11.j <- omega.j.inv
-    I22.j <- 0.5 * lav_matrix_duplication_pre_post(omega.j.inv %x% omega.j.inv)
+    if (lav_use_lavaanC()) {
+      I22.j <- lavaanC::m_kronecker_dup_pre_post(omega.j.inv, multiplicator = 0.5)
+    } else {
+      I22.j <- 0.5 * lav_matrix_duplication_pre_post(omega.j.inv %x% omega.j.inv)
+    }
     I.j <- lav_matrix_bdiag(I11.j, I22.j)
     info.j <- t(Delta.j) %*% I.j %*% Delta.j
 
@@ -1173,7 +1185,11 @@ lav_mvnorm_cluster_information_expected_delta <- function(Lp = NULL,
     Sinv.method = Sinv.method
   )
   I11.w <- Sigma.W.inv
-  I22.w <- 0.5 * lav_matrix_duplication_pre_post(Sigma.W.inv %x% Sigma.W.inv)
+  if (lav_use_lavaanC()) {
+    I22.w <- lavaanC::m_kronecker_dup_pre_post(Sigma.W.inv, multiplicator = 0.5)
+  } else {
+    I22.w <- 0.5 * lav_matrix_duplication_pre_post(Sigma.W.inv %x% Sigma.W.inv)
+  }
   I.w <- lav_matrix_bdiag(I11.w, I22.w)
 
   # force zero for means both.idx in within part
diff --git a/R/lav_mvnorm_h1.R b/R/lav_mvnorm_h1.R
index c50590c7..8e268940 100644
--- a/R/lav_mvnorm_h1.R
+++ b/R/lav_mvnorm_h1.R
@@ -262,11 +262,21 @@ lav_mvnorm_h1_information_expected <- function(
 
   I11 <- sample.cov.inv
   if(correlation) {
-    I22 <- 0.5 * lav_matrix_duplication_cor_pre_post(sample.cov.inv %x%
+    if (lav_use_lavaanC()) {
+      I22 <- lavaanC::m_kronecker_dup_cor_pre_post(sample.cov.inv, 
+                                                   multiplicator = 0.5)
+    } else {
+      I22 <- 0.5 * lav_matrix_duplication_cor_pre_post(sample.cov.inv %x%
                                                      sample.cov.inv)
+    }
   } else {
-    I22 <- 0.5 * lav_matrix_duplication_pre_post(sample.cov.inv %x%
+    if (lav_use_lavaanC()) {
+      I22 <- lavaanC::m_kronecker_dup_pre_post(sample.cov.inv,
+                                               multiplicator = 0.5)
+    } else {
+      I22 <- 0.5 * lav_matrix_duplication_pre_post(sample.cov.inv %x%
                                                    sample.cov.inv)
+    }
   }
 
   # fixed.x?
@@ -333,9 +343,14 @@ lav_mvnorm_h1_information_observed_samplestats <- function(
     I11[, x.idx] <- 0
   }
 
-  I22 <- 0.5 * lav_matrix_duplication_pre_post(sample.cov.inv %x%
+  if (lav_use_lavaanC()) {
+    I22 <- lavaanC::m_kronecker_dup_pre_post(sample.cov.inv, 
+                                             multiplicator = 0.5)
+  } else {
+    I22 <- 0.5 * lav_matrix_duplication_pre_post(sample.cov.inv %x%
     sample.cov.inv)
-
+  }
+  
   # fixed.x?
   if (length(x.idx) > 0L) {
     pstar.x <- lav_matrix_vech_which_idx(
@@ -470,8 +485,12 @@ lav_mvnorm_h1_inverted_information_observed <- function(
     )
   } else {
     I11 <- sample.cov
-    I22 <- 2 * lav_matrix_duplication_ginv_pre_post(sample.cov %x% sample.cov)
-
+    if (lav_use_lavaanC()) {
+      I22 <- lavaanC::m_kronecker_dup_ginv_pre_post(sample.cov, 
+                                                    multiplicator = 2.0)
+    } else {
+      I22 <- 2 * lav_matrix_duplication_ginv_pre_post(sample.cov %x% sample.cov)
+    }
     Gamma.NT <- lav_matrix_bdiag(I11, I22)
   }
 
diff --git a/R/lav_mvnorm_missing.R b/R/lav_mvnorm_missing.R
index 2af9d744..ac39b851 100644
--- a/R/lav_mvnorm_missing.R
+++ b/R/lav_mvnorm_missing.R
@@ -1037,8 +1037,11 @@ lav_mvnorm_missing_logl_hessian_samplestats <-
 
       i11 <- S.inv
       i21 <- lav_matrix_duplication_pre(tmp21 %x% S.inv)
-      i22 <- (1 / 2) * lav_matrix_duplication_pre_post(S.inv %x% tmp22)
-
+      if (lav_use_lavaanC()) {
+        i22 <- lavaanC::m_kronecker_dup_pre_post(S.inv, tmp22, 0.5)
+      } else {
+        i22 <- (1 / 2) * lav_matrix_duplication_pre_post(S.inv %x% tmp22)
+      }
       H11 <- H11 + pat.freq * i11
       H21 <- H21 + pat.freq * i21
       H22 <- H22 + pat.freq * i22
@@ -1128,8 +1131,12 @@ lav_mvnorm_missing_information_expected <- function(Y = NULL,
     S.inv <- matrix(0, P, P)
     S.inv[var.idx, var.idx] <- sigma.inv
 
-    S2.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
-
+    if (lav_use_lavaanC()) {
+      S2.inv <- lavaanC::m_kronecker_dup_pre_post(S.inv, multiplicator = 0.5)
+    } else {
+      S2.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
+    }
+    
     if (!is.null(wt)) {
       FREQ <- sum(wt[Mp$case.idx[[p]]])
     } else {
@@ -1360,8 +1367,12 @@ lav_mvnorm_missing_information_both <- function(Y = NULL,
     }
 
     if (information == "expected") {
-      S2.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
-
+      if (lav_use_lavaanC()) {
+        S2.inv <- lavaanC::m_kronecker_dup_pre_post(S.inv, multiplicator = 0.5)
+      } else {
+        S2.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
+      } 
+      
       I11 <- I11 + FREQ * S.inv
       I22 <- I22 + FREQ * S2.inv
     } else {
@@ -1379,8 +1390,12 @@ lav_mvnorm_missing_information_both <- function(Y = NULL,
 
       i11 <- S.inv
       i21 <- lav_matrix_duplication_pre(tmp21 %x% S.inv)
-      i22 <- (1 / 2) * lav_matrix_duplication_pre_post(S.inv %x% tmp22)
-
+      if (lav_use_lavaanC()) {
+        i22 <- lavaanC::m_kronecker_dup_pre_post(S.inv, tmp22, 0.5)
+      } else {
+        i22 <- (1 / 2) * lav_matrix_duplication_pre_post(S.inv %x% tmp22)
+      }
+      
       I11 <- I11 + pat.freq * i11
       I21 <- I21 + pat.freq * i21
       I22 <- I22 + pat.freq * i22
diff --git a/R/lav_mvreg.R b/R/lav_mvreg.R
index 644fa2c1..9608aca9 100644
--- a/R/lav_mvreg.R
+++ b/R/lav_mvreg.R
@@ -517,8 +517,12 @@ lav_mvreg_information_expected <- function(Y = NULL, # not used
 
   # expected information
   I11 <- res.cov.inv %x% sample.xx
-  I22 <- 0.5 * lav_matrix_duplication_pre_post(res.cov.inv %x% res.cov.inv)
-
+  if (lav_use_lavaanC()) {
+    I22 <- lavaanC::m_kronecker_dup_pre_post(res.cov.inv, multiplicator = 0.5)
+  } else {
+    I22 <- 0.5 * lav_matrix_duplication_pre_post(res.cov.inv %x% res.cov.inv)
+  }
+  
   lav_matrix_bdiag(I11, I22)
 }
 
@@ -605,8 +609,12 @@ lav_mvreg_information_observed_samplestats <-
     H12 <- t(H21)
 
     AAA <- res.cov.inv %*% (2 * W.tilde - res.cov) %*% res.cov.inv
-    H22 <- (1 / 2) * lav_matrix_duplication_pre_post(res.cov.inv %x% AAA)
-
+    if (lav_use_lavaanC()) {
+      H22 <- lavaanC::m_kronecker_dup_pre_post(res.cov.inv, AAA, 0.5)
+    } else {
+      H22 <- (1 / 2) * lav_matrix_duplication_pre_post(res.cov.inv %x% AAA)
+    }
+    
     out <- rbind(
       cbind(H11, H12),
       cbind(H21, H22)
diff --git a/R/lav_samplestats_gamma.R b/R/lav_samplestats_gamma.R
index 224fcb08..68e694cc 100644
--- a/R/lav_samplestats_gamma.R
+++ b/R/lav_samplestats_gamma.R
@@ -323,7 +323,11 @@ lav_samplestats_Gamma_NT <- function(Y = NULL, # should include
   if (!conditional.x) {
     # unconditional - stochastic x
     if (!fixed.x) {
-      Gamma <- 2 * lav_matrix_duplication_ginv_pre_post(S %x% S)
+      if (lav_use_lavaanC()) {
+        Gamma <- lavaanC::m_kronecker_dup_ginv_pre_post(S, multiplicator = 2.0)
+      } else {
+        Gamma <- 2 * lav_matrix_duplication_ginv_pre_post(S %x% S)
+      }
       if (meanstructure) {
         Gamma <- lav_matrix_bdiag(S, Gamma)
       }
@@ -346,8 +350,13 @@ lav_samplestats_Gamma_NT <- function(Y = NULL, # should include
       # take difference
       R <- S - YbarX.aug
 
-      Gamma.S <- 2 * lav_matrix_duplication_ginv_pre_post(S %x% S)
-      Gamma.R <- 2 * lav_matrix_duplication_ginv_pre_post(R %x% R)
+      if (lav_use_lavaanC()) {
+        Gamma.S <- lavaanC::m_kronecker_dup_ginv_pre_post(S, multiplicator = 2.0)
+        Gamma.R <- lavaanC::m_kronecker_dup_ginv_pre_post(R, multiplicator = 2.0)
+      } else {
+        Gamma.S <- 2 * lav_matrix_duplication_ginv_pre_post(S %x% S)
+        Gamma.R <- 2 * lav_matrix_duplication_ginv_pre_post(R %x% R)
+      }
       Gamma <- Gamma.S - Gamma.R
 
       if (meanstructure) {
@@ -368,8 +377,12 @@ lav_samplestats_Gamma_NT <- function(Y = NULL, # should include
     B <- S[-x.idx, x.idx, drop = FALSE]
     C <- S[x.idx, x.idx, drop = FALSE]
     Cov.YbarX <- A - B %*% solve(C) %*% t(B)
-    Gamma <- 2 * lav_matrix_duplication_ginv_pre_post(Cov.YbarX %x% Cov.YbarX)
-
+    if (lav_use_lavaanC()) {
+      Gamma <- lavaanC::m_kronecker_dup_ginv_pre_post(Cov.YbarX, multiplicator = 2.0)
+    } else {
+      Gamma <- 2 * lav_matrix_duplication_ginv_pre_post(Cov.YbarX %x% Cov.YbarX)
+    }
+    
     if (meanstructure || slopestructure) {
       MY <- M[-x.idx]
       MX <- M[x.idx]
@@ -716,8 +729,13 @@ lav_samplestats_Gamma <- function(Y, # Y+X if cond!
   # unbiased?
   if (unbiased) {
     # normal-theory Gamma (cov only)
-    GammaNT.cov <- 2 * lav_matrix_duplication_ginv_pre_post(COV %x% COV)
-
+    if (lav_use_lavaanC()) {
+      GammaNT.cov <- lavaanC::m_kronecker_dup_ginv_pre_post(COV, 
+                                                            multiplicator = 2.0)
+    } else {
+      GammaNT.cov <- 2 * lav_matrix_duplication_ginv_pre_post(COV %x% COV)
+    }
+    
     if (meanstructure) {
       Gamma.cov <- Gamma[-(1:p), -(1:p), drop = FALSE]
       Gamma.mean.cov <- Gamma[1:p, -(1:p), drop = FALSE]
diff --git a/R/lav_samplestats_igamma.R b/R/lav_samplestats_igamma.R
index 6cb0954c..08de5525 100644
--- a/R/lav_samplestats_igamma.R
+++ b/R/lav_samplestats_igamma.R
@@ -75,7 +75,11 @@ lav_samplestats_Gamma_inverse_NT <- function(Y = NULL,
   if (!conditional.x) {
     # unconditional - stochastic x
     if (!fixed.x) {
-      Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
+      if (lav_use_lavaanC()) {
+        Gamma.inv <- lavaanC::m_kronecker_dup_pre_post(S.inv, multiplicator = 0.5)
+      } else {
+        Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
+      }
       if (meanstructure) {
         Gamma.inv <- lav_matrix_bdiag(S.inv, Gamma.inv)
       }
@@ -83,8 +87,12 @@ lav_samplestats_Gamma_inverse_NT <- function(Y = NULL,
       # unconditional - fixed x
     } else {
       # handle fixed.x = TRUE
-      Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
-
+      if (lav_use_lavaanC()) {
+        Gamma.inv <- lavaanC::m_kronecker_dup_pre_post(S.inv, multiplicator = 0.5)
+      } else {
+        Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
+      }
+      
       # zero rows/cols corresponding with x/x combinations
       nvar <- NROW(ICOV)
       pstar <- nvar * (nvar + 1) / 2
@@ -112,8 +120,12 @@ lav_samplestats_Gamma_inverse_NT <- function(Y = NULL,
 
     S11 <- S.inv[-x.idx, -x.idx, drop = FALSE]
 
-    Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S11 %x% S11)
-
+    if (lav_use_lavaanC()) {
+      Gamma.inv <- lavaanC::m_kronecker_dup_pre_post(S11, multiplicator = 0.5)
+    } else {
+      Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S11 %x% S11)
+    }
+    
     if (meanstructure || slopestructure) {
       C <- S[x.idx, x.idx, drop = FALSE]
       MY <- M[-x.idx]