diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp
index f7ef57546f..c46ace3f73 100644
--- a/src/cmfd_solver.cpp
+++ b/src/cmfd_solver.cpp
@@ -101,6 +101,7 @@ int cmfd_linsolver_1g(const double* A_data, const double* b, double* x,
     for (int irb = 0; irb < 2; irb++) {
 
       // Loop around matrix rows
+#pragma omp parallel for
       for (int irow = 0; irow < cmfd::dim; irow++) {
         int g, i, j, k;
         matrix_to_indices(irow, g, i, j, k);
@@ -167,6 +168,7 @@ int cmfd_linsolver_2g(const double* A_data, const double* b, double* x,
     for (int irb = 0; irb < 2; irb++) {
 
       // Loop around matrix rows
+#pragma omp parallel for
       for (int irow = 0; irow < cmfd::dim; irow+=2) {
         int g, i, j, k;
         matrix_to_indices(irow, g, i, j, k);
@@ -255,6 +257,7 @@ int cmfd_linsolver_ng(const double* A_data, const double* b, double* x,
     std::vector<double> tmpx {x, x+cmfd::dim};
 
     // Loop around matrix rows
+#pragma omp parallel for
     for (int irow = 0; irow < cmfd::dim; irow++) {
       // Get index of diagonal for current row
       int didx = get_diagonal_index(irow);