diff --git a/anharmonic/BTE_RTA.py b/anharmonic/BTE_RTA.py
index 2e83db19..32386260 100644
--- a/anharmonic/BTE_RTA.py
+++ b/anharmonic/BTE_RTA.py
@@ -82,8 +82,10 @@ class BTE_RTA:
              self._grid_address) = spg.get_ir_reciprocal_mesh(self._mesh,
                                                               self._primitive)
             dense_grid_points = np.unique(grid_mapping_table)
-            dense_grid_weights = [np.sum(grid_mapping_table == g)
-                                  for g in dense_grid_points]
+            weights = np.zeros_like(grid_mapping_table)
+            for g in grid_mapping_table:
+                weights[g] += 1
+            dense_grid_weights = weights[dense_grid_points]
 
             self._grid_points, self._grid_weights = reduce_grid_points(
                 self._mesh_divisors,
diff --git a/anharmonic/triplets.py b/anharmonic/triplets.py
index a741d413..b2e3e630 100644
--- a/anharmonic/triplets.py
+++ b/anharmonic/triplets.py
@@ -44,9 +44,9 @@ def get_grid_address(mesh):
     for i in range(mesh[2]):
         for j in range(mesh[1]):
             for k in range(mesh[0]):
-                grid_address[count] = [k - (k > (mesh[2] // 2)) * mesh[0],
+                grid_address[count] = [k - (k > (mesh[0] // 2)) * mesh[0],
                                        j - (j > (mesh[1] // 2)) * mesh[1],
-                                       i - (i > (mesh[0] // 2)) * mesh[2]]
+                                       i - (i > (mesh[2] // 2)) * mesh[2]]
 
                 count += 1
     
@@ -77,9 +77,11 @@ def reduce_grid_points(mesh_divisors,
     else:
         grid_points = []
         grid_weights = []
-        for i, dgp in enumerate(dense_grid_points):
-            if (grid_address[dgp] % divisors == 0).all():
-                grid_points.append(dgp)
+
+        for i, modulo in enumerate(
+            grid_address[dense_grid_points] % divisors):
+            if (modulo == 0).all():
+                grid_points.append(dense_grid_points[i])
                 if dense_grid_weights is not None:
                     grid_weights.append(dense_grid_weights[i])
 
diff --git a/c/spglib/kpoint.c b/c/spglib/kpoint.c
index 5ed5fe59..bbb25c16 100644
--- a/c/spglib/kpoint.c
+++ b/c/spglib/kpoint.c
@@ -780,10 +780,10 @@ static int get_ir_triplets_at_q(int weights[],
 				const int mesh[3],
 				SPGCONST PointSymmetry * pointgroup)
 {
-  int i, j, num_grid, weight_q, q_2, num_ir_q, num_ir_triplets;
+  int i, j, num_grid, q_2, num_ir_q, num_ir_triplets, ir_address;
   int mesh_double[3], is_shift[3];
   int grid_double0[3], grid_double1[3], grid_double2[3];
-  int *map_q, *ir_addresses;
+  int *map_q, *ir_addresses, *weight_q;
   double tolerance;
   double stabilizer_q[1][3];
   PointSymmetry pointgroup_q;
@@ -825,14 +825,21 @@ static int get_ir_triplets_at_q(int weights[],
 #endif
 
   ir_addresses = (int*) malloc(sizeof(int) * num_ir_q);
+  weight_q = (int*) malloc(sizeof(int) * num_grid);
   num_ir_q = 0;
   for (i = 0; i < num_grid; i++) {
     if (map_q[i] == i) {
       ir_addresses[num_ir_q] = i;
       num_ir_q++;
     }
-    weights[i] = 0;
+    weight_q[i] = 0;
+  }
+
+#pragma omp parallel for
+  for (i = 0; i < num_grid; i++) {
+    weight_q[map_q[i]]++;
     third_q[i] = -1;
+    weights[i] = 0;
   }
 
 #pragma omp parallel for private(j, grid_double1, grid_double2)
@@ -847,26 +854,20 @@ static int get_ir_triplets_at_q(int weights[],
 
   num_ir_triplets = 0;
   for (i = 0; i < num_ir_q; i++) {
-    q_2 = third_q[ir_addresses[i]];
-    weight_q = 0;
-    
-#pragma omp parallel for reduction(+:weight_q)
-    for (j = 0; j < num_grid; j++) {
-      if (ir_addresses[i] == map_q[j]) {
-	weight_q++;
-      }
-    }
-
+    ir_address = ir_addresses[i];
+    q_2 = third_q[ir_address];
     if (weights[map_q[q_2]]) {
-      weights[map_q[q_2]] += weight_q;
+      weights[map_q[q_2]] += weight_q[ir_address];
     } else {
-      weights[ir_addresses[i]] = weight_q;
+      weights[ir_address] = weight_q[ir_address];
       num_ir_triplets++;
     }
   }
 
   free(map_q);
   map_q = NULL;
+  free(weight_q);
+  weight_q = NULL;
   free(ir_addresses);
   ir_addresses = NULL;