diff --git a/ctest/test_gridsys.cpp b/ctest/test_gridsys.cpp
index 2d1e8306..fffec772 100644
--- a/ctest/test_gridsys.cpp
+++ b/ctest/test_gridsys.cpp
@@ -2,6 +2,7 @@
 
 extern "C" {
 #include "gridsys.h"
+#include "utils.h"
 }
 
 /**
@@ -11,13 +12,13 @@ extern "C" {
 TEST(test_gridsys, test_gridsys_get_all_grid_addresses) {
     long(*gr_grid_addresses)[3];
     long D_diag[3] = {3, 4, 5};
-    long n, i, j, k, grid_index;
+    long n, i, j, k;
+    long grid_index = 0;
 
     n = D_diag[0] * D_diag[1] * D_diag[2];
     gr_grid_addresses = (long(*)[3])malloc(sizeof(long[3]) * n);
     gridsys_get_all_grid_addresses(gr_grid_addresses, D_diag);
 
-    grid_index = 0;
     for (k = 0; k < D_diag[2]; k++) {
         for (j = 0; j < D_diag[1]; j++) {
             for (i = 0; i < D_diag[0]; i++) {
@@ -60,7 +61,8 @@ TEST(test_gridsys, test_gridsys_get_double_grid_address) {
 TEST(test_gridsys, test_gridsys_get_grid_address_from_index) {
     long address[3];
     long D_diag[3] = {3, 4, 5};
-    long i, j, k, l, grid_index;
+    long i, j, k, l;
+    long grid_index = 0;
 
     for (k = 0; k < D_diag[2]; k++) {
         for (j = 0; j < D_diag[1]; j++) {
@@ -78,6 +80,135 @@ TEST(test_gridsys, test_gridsys_get_grid_address_from_index) {
     }
 }
 
+/**
+ * @brief gridsys_get_double_grid_index
+ * Return grid point index corresponding to double grid address in GR-grid.
+ */
+TEST(test_gridsys, test_gridsys_get_double_grid_index) {
+    long address_double[3], address[3];
+    long D_diag[3] = {3, 4, 5};
+    long PS[8][3] = {{0, 0, 0}, {0, 0, 3}, {0, 3, 0}, {0, 3, 3},
+                     {3, 0, 0}, {3, 0, 3}, {3, 3, 0}, {3, 3, 3}};
+    long i, j, k, l, grid_index;
+
+    for (l = 0; l < 8; l++) {
+        grid_index = 0;
+        for (k = 0; k < D_diag[2]; k++) {
+            address[2] = k;
+            for (j = 0; j < D_diag[1]; j++) {
+                address[1] = j;
+                for (i = 0; i < D_diag[0]; i++) {
+                    address[0] = i;
+                    gridsys_get_double_grid_address(address_double, address,
+                                                    PS[l]);
+                    ASSERT_EQ(grid_index, gridsys_get_double_grid_index(
+                                              address_double, D_diag, PS[l]));
+                    grid_index++;
+                }
+            }
+        }
+    }
+}
+
+/**
+ * @brief gridsys_get_grid_index_from_address
+ * Return grid point index corresponding to grid address in GR-grid.
+ */
+TEST(test_gridsys, test_gridsys_get_grid_index_from_address) {
+    long address[3];
+    long D_diag[3] = {3, 4, 5};
+    long i, j, k;
+    long grid_index = 0;
+
+    for (k = 0; k < D_diag[2]; k++) {
+        address[2] = k;
+        for (j = 0; j < D_diag[1]; j++) {
+            address[1] = j;
+            for (i = 0; i < D_diag[0]; i++) {
+                address[0] = i;
+                ASSERT_EQ(grid_index,
+                          gridsys_get_grid_index_from_address(address, D_diag));
+                grid_index++;
+            }
+        }
+    }
+}
+
+/**
+ * @brief gridsys_rotate_grid_index
+ * Return grid point index of rotated address of given grid point index.
+ */
+TEST(test_gridsys, test_gridsys_rotate_grid_index) {
+    long address[3], rot_address[3];
+    long D_diag[3] = {3, 4, 5};
+    long PS[8][3] = {{0, 0, 0}, {0, 0, 3}, {0, 3, 0}, {0, 3, 3},
+                     {3, 0, 0}, {3, 0, 3}, {3, 3, 0}, {3, 3, 3}};
+    long i, j, k, l;
+    long grid_index = 0;
+    long rot_grid_index;
+    long rot_grid_indices[8][60] = {
+        {0,  9,  6,  1,  10, 7,  2,  11, 8,  0,  9,  6,  48, 57, 54,
+         49, 58, 55, 50, 59, 56, 48, 57, 54, 36, 45, 42, 37, 46, 43,
+         38, 47, 44, 36, 45, 42, 24, 33, 30, 25, 34, 31, 26, 35, 32,
+         24, 33, 30, 12, 21, 18, 13, 22, 19, 14, 23, 20, 12, 21, 18},
+        {24, 33, 30, 25, 34, 31, 26, 35, 32, 24, 33, 30, 12, 21, 18,
+         13, 22, 19, 14, 23, 20, 12, 21, 18, 0,  9,  6,  1,  10, 7,
+         2,  11, 8,  0,  9,  6,  48, 57, 54, 49, 58, 55, 50, 59, 56,
+         48, 57, 54, 36, 45, 42, 37, 46, 43, 38, 47, 44, 36, 45, 42},
+        {10, 7,  4,  11, 8,  5,  9,  6,  3,  10, 7,  4,  58, 55, 52,
+         59, 56, 53, 57, 54, 51, 58, 55, 52, 46, 43, 40, 47, 44, 41,
+         45, 42, 39, 46, 43, 40, 34, 31, 28, 35, 32, 29, 33, 30, 27,
+         34, 31, 28, 22, 19, 16, 23, 20, 17, 21, 18, 15, 22, 19, 16},
+        {34, 31, 28, 35, 32, 29, 33, 30, 27, 34, 31, 28, 22, 19, 16,
+         23, 20, 17, 21, 18, 15, 22, 19, 16, 10, 7,  4,  11, 8,  5,
+         9,  6,  3,  10, 7,  4,  58, 55, 52, 59, 56, 53, 57, 54, 51,
+         58, 55, 52, 46, 43, 40, 47, 44, 41, 45, 42, 39, 46, 43, 40},
+        {11, 8,  5,  9,  6,  3,  9,  6,  3,  10, 7,  4,  59, 56, 53,
+         57, 54, 51, 57, 54, 51, 58, 55, 52, 47, 44, 41, 45, 42, 39,
+         45, 42, 39, 46, 43, 40, 35, 32, 29, 33, 30, 27, 33, 30, 27,
+         34, 31, 28, 23, 20, 17, 21, 18, 15, 21, 18, 15, 22, 19, 16},
+        {35, 32, 29, 33, 30, 27, 33, 30, 27, 34, 31, 28, 23, 20, 17,
+         21, 18, 15, 21, 18, 15, 22, 19, 16, 11, 8,  5,  9,  6,  3,
+         9,  6,  3,  10, 7,  4,  59, 56, 53, 57, 54, 51, 57, 54, 51,
+         58, 55, 52, 47, 44, 41, 45, 42, 39, 45, 42, 39, 46, 43, 40},
+        {3,  0,  9,  4,  1,  10, 5,  2,  11, 3,  0,  9,  51, 48, 57,
+         52, 49, 58, 53, 50, 59, 51, 48, 57, 39, 36, 45, 40, 37, 46,
+         41, 38, 47, 39, 36, 45, 27, 24, 33, 28, 25, 34, 29, 26, 35,
+         27, 24, 33, 15, 12, 21, 16, 13, 22, 17, 14, 23, 15, 12, 21},
+        {27, 24, 33, 28, 25, 34, 29, 26, 35, 27, 24, 33, 15, 12, 21,
+         16, 13, 22, 17, 14, 23, 15, 12, 21, 3,  0,  9,  4,  1,  10,
+         5,  2,  11, 3,  0,  9,  51, 48, 57, 52, 49, 58, 53, 50, 59,
+         51, 48, 57, 39, 36, 45, 40, 37, 46, 41, 38, 47, 39, 36, 45}};
+
+    // Rutile R^T of a screw operation.
+    long rotation[3][3] = {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}};
+
+    grid_index = 0;
+    for (k = 0; k < D_diag[2]; k++) {
+        address[2] = k;
+        for (j = 0; j < D_diag[1]; j++) {
+            address[1] = j;
+            for (i = 0; i < D_diag[0]; i++) {
+                address[0] = i;
+                lagmat_multiply_matrix_vector_l3(rot_address, rotation,
+                                                 address);
+                ASSERT_EQ(
+                    (gridsys_get_grid_index_from_address(rot_address, D_diag)),
+                    (gridsys_rotate_grid_index(grid_index, rotation, D_diag,
+                                               PS[0])));
+                grid_index++;
+            }
+        }
+    }
+    for (l = 0; l < 8; l++) {
+        for (grid_index = 0; grid_index < 60; grid_index++) {
+            ASSERT_EQ(
+                rot_grid_indices[l][grid_index],
+                gridsys_rotate_grid_index(grid_index, rotation, D_diag, PS[l]));
+        }
+    }
+}
+
 /* TEST(test_gridsys, test_gridsys_get_all_grid_addresses) {
     double lattice[3][3] = {{4, 0, 0}, {0, 4, 0}, {0, 0, 3}};
     double position[][3] = {
diff --git a/ctest/utils.c b/ctest/utils.c
new file mode 100644
index 00000000..7297de24
--- /dev/null
+++ b/ctest/utils.c
@@ -0,0 +1,37 @@
+#include <stdio.h>
+
+void lagmat_multiply_matrix_vector_l3(long v[3], const long a[3][3],
+                                      const long b[3]) {
+    long i;
+    long c[3];
+    for (i = 0; i < 3; i++) {
+        c[i] = a[i][0] * b[0] + a[i][1] * b[1] + a[i][2] * b[2];
+    }
+    for (i = 0; i < 3; i++) {
+        v[i] = c[i];
+    }
+}
+
+void lagmat_multiply_matrix_l3(long m[3][3], const long a[3][3],
+                               const long b[3][3]) {
+    long i, j; /* a_ij */
+    long c[3][3];
+    for (i = 0; i < 3; i++) {
+        for (j = 0; j < 3; j++) {
+            c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
+        }
+    }
+    lagmat_copy_matrix_l3(m, c);
+}
+
+void lagmat_copy_matrix_l3(long a[3][3], const long b[3][3]) {
+    a[0][0] = b[0][0];
+    a[0][1] = b[0][1];
+    a[0][2] = b[0][2];
+    a[1][0] = b[1][0];
+    a[1][1] = b[1][1];
+    a[1][2] = b[1][2];
+    a[2][0] = b[2][0];
+    a[2][1] = b[2][1];
+    a[2][2] = b[2][2];
+}
diff --git a/ctest/utils.h b/ctest/utils.h
new file mode 100644
index 00000000..a44f77cd
--- /dev/null
+++ b/ctest/utils.h
@@ -0,0 +1,10 @@
+#ifndef __test_utils_H__
+#define __test_utils_H__
+
+void lagmat_multiply_matrix_vector_l3(long v[3], const long a[3][3],
+                                      const long b[3]);
+void lagmat_multiply_matrix_l3(long m[3][3], const long a[3][3],
+                               const long b[3][3]);
+void lagmat_copy_matrix_l3(long a[3][3], const long b[3][3]);
+
+#endif  // __test_utils_H__