diff --git a/src/add_covariances.cpp b/src/add_covariances.cpp
index 7435983..2c51d0f 100644
--- a/src/add_covariances.cpp
+++ b/src/add_covariances.cpp
@@ -50,8 +50,10 @@ void add_covariances(std::vector<std::string> variables,
   bool covariance_exists = false;
   for(unsigned int i = 0; i < exogenous.size()-1; i++){
     for(unsigned int j = i+1; j < exogenous.size(); j++){
+
       covariance_exists = false;
       for(unsigned int k = 0; k < pt.lhs.size(); k++){
+
         if((pt.lhs.at(k).compare(exogenous.at(i)) == 0) &&
            (pt.op.at(k).compare("~~") == 0) &&
            (pt.rhs.at(k).compare(exogenous.at(j)) == 0)
@@ -59,6 +61,14 @@ void add_covariances(std::vector<std::string> variables,
           covariance_exists = true;
           break;
         }
+        // the order does not matter:
+        if((pt.lhs.at(k).compare(exogenous.at(j)) == 0) &&
+           (pt.op.at(k).compare("~~") == 0) &&
+           (pt.rhs.at(k).compare(exogenous.at(i)) == 0)
+        ){
+          covariance_exists = true;
+          break;
+        }
       }
       if(!covariance_exists){
         pt.add_line();
diff --git a/tests/testthat/test-labeled-covariances.R b/tests/testthat/test-labeled-covariances.R
new file mode 100644
index 0000000..6a14f9d
--- /dev/null
+++ b/tests/testthat/test-labeled-covariances.R
@@ -0,0 +1,43 @@
+test_that("labeled covariances", {
+  library(mxsem)
+
+  model <- '
+  # latent variable definitions
+     ind60 =~ x1 + x2 + x3
+     dem60 =~ y1 + a1*y2 + b*y3 + c1*y4
+     dem65 =~ y5 + a2*y6 + b*y7
+    ind60 ~y8
+    ind60 ~~ v1 * ind60
+    dem60 ~~ v2 * ind60 + v3*dem60
+    dem65 ~~ v4 * ind60 + v5*dem60 + v6*dem65
+'
+
+  fit <- mxsem(model = model,
+               data  = OpenMx::Bollen) |>
+    mxTryHard()
+
+  for(i in paste0("v", 1:6)){
+    testthat::expect_true(i %in% names(omxGetParameters(fit)))
+  }
+
+
+  model <- '
+  # latent variable definitions
+     ind60 =~ x1 + x2 + x3
+     dem60 =~ y1 + a1*y2 + b*y3 + c1*y4
+     dem65 =~ y5 + a2*y6 + b*y7
+    ind60 ~y8
+    ind60 ~~ v1 * ind60
+    dem60 ~~ v2 * ind60 + v3*dem60
+    dem65 ~~ v4 * ind60 + v5*dem60 + v6*dem65
+
+    x1 ~ int1*1; x2 ~ int1*1
+'
+
+  fit <- mxsem(model = model,
+               data  = OpenMx::Bollen) |>
+    mxTryHard()
+
+  testthat::expect_true("int1" %in% names(omxGetParameters(fit)))
+
+})