diff --git a/.github/workflows/run_test_vs.yml b/.github/workflows/run_test_vs.yml
index 97f7dc2..19e9be4 100644
--- a/.github/workflows/run_test_vs.yml
+++ b/.github/workflows/run_test_vs.yml
@@ -2,7 +2,7 @@ name: Run test VS
 
 on:
   push:
-    branches: [ develop ]
+    branches: [ develop, feature/* ]
 
 jobs:
   build:
diff --git a/python_numpy/python_numpy_linalg_solver.hpp b/python_numpy/python_numpy_linalg_solver.hpp
index 7c8b964..d4fc598 100644
--- a/python_numpy/python_numpy_linalg_solver.hpp
+++ b/python_numpy/python_numpy_linalg_solver.hpp
@@ -164,8 +164,9 @@ class LinalgSolver {
         division_min(std::move(other.division_min)), rho(std::move(other.rho)),
         rep_num(std::move(other.rep_num)) {}
 
-  LinalgSolver<T, M, K, SparseAvailable_A, SparseAvailable_B> &operator=(
-      LinalgSolver<T, M, K, SparseAvailable_A, SparseAvailable_B> &&other) {
+  LinalgSolver<T, M, K, SparseAvailable_A, SparseAvailable_B> &
+  operator=(LinalgSolver<T, M, K, SparseAvailable_A, SparseAvailable_B>
+                &&other) noexcept {
     if (this != &other) {
       this->X_1 = std::move(other.X_1);
       this->decay_rate = std::move(other.decay_rate);
@@ -626,7 +627,7 @@ class LinalgPartitionSolver {
 public:
   /* Constructor */
   LinalgPartitionSolver()
-      : decay_rate(static_cast<Value_Type>(0)),
+      : X_1(), decay_rate(static_cast<Value_Type>(0)),
         division_min(
             static_cast<Value_Type>(DEFAULT_DIVISION_MIN_LINALG_SOLVER)),
         rho({}), rep_num({}) {}
@@ -634,13 +635,15 @@ class LinalgPartitionSolver {
   /* Copy Constructor */
   LinalgPartitionSolver(const LinalgPartitionSolver<T, M, K, SparseAvailable_A,
                                                     SparseAvailable_B> &other)
-      : decay_rate(other.decay_rate), division_min(other.division_min),
-        rho(other.rho), rep_num(other.rep_num) {}
+      : X_1(other.X_1), decay_rate(other.decay_rate),
+        division_min(other.division_min), rho(other.rho),
+        rep_num(other.rep_num) {}
 
   LinalgPartitionSolver<T, M, K, SparseAvailable_A, SparseAvailable_B> &
   operator=(const LinalgPartitionSolver<T, M, K, SparseAvailable_A,
                                         SparseAvailable_B> &other) {
     if (this != &other) {
+      this->X_1 = other.X_1;
       this->decay_rate = other.decay_rate;
       this->division_min = other.division_min;
       this->rho = other.rho;
@@ -653,7 +656,7 @@ class LinalgPartitionSolver {
   LinalgPartitionSolver(
       LinalgPartitionSolver<T, M, K, SparseAvailable_A, SparseAvailable_B>
           &&other) noexcept
-      : decay_rate(std::move(other.decay_rate)),
+      : X_1(std::move(other.X_1)), decay_rate(std::move(other.decay_rate)),
         division_min(std::move(other.division_min)), rho(std::move(other.rho)),
         rep_num(std::move(other.rep_num)) {}
 
@@ -661,6 +664,7 @@ class LinalgPartitionSolver {
   operator=(LinalgPartitionSolver<T, M, K, SparseAvailable_A, SparseAvailable_B>
                 &&other) noexcept {
     if (this != &other) {
+      this->X_1 = std::move(other.X_1);
       this->decay_rate = std::move(other.decay_rate);
       this->division_min = std::move(other.division_min);
       this->rho = std::move(other.rho);
@@ -674,41 +678,89 @@ class LinalgPartitionSolver {
                     const Matrix<DefDense, T, M, K> &B, std::size_t matrix_size)
       -> Matrix<DefDense, T, M, K> {
 
-    Base::Matrix::Matrix<T, M, K> X_1;
+    if (matrix_size > M) {
+      matrix_size = M;
+    }
+
+    Base::Matrix::gmres_k_partition_matrix(
+        A.matrix, B.matrix, this->X_1, this->decay_rate, this->division_min,
+        this->rho, this->rep_num, matrix_size);
+
+    return Matrix<DefDense, T, M, K>(this->X_1);
+  }
+
+  inline auto cold_solve(const Matrix<DefDense, T, M, M> &A,
+                         const Matrix<DefDense, T, M, K> &B,
+                         std::size_t matrix_size) -> Matrix<DefDense, T, M, K> {
+
+    Base::Matrix::Matrix<T, M, K> X_1_temporary;
 
     if (matrix_size > M) {
       matrix_size = M;
     }
 
     Base::Matrix::gmres_k_partition_matrix(
-        A.matrix, B.matrix, X_1, this->decay_rate, this->division_min,
+        A.matrix, B.matrix, X_1_temporary, this->decay_rate, this->division_min,
         this->rho, this->rep_num, matrix_size);
 
-    return Matrix<DefDense, T, M, K>(X_1);
+    return Matrix<DefDense, T, M, K>(X_1_temporary);
   }
 
   inline auto solve(const Matrix<DefDense, T, M, M> &A,
                     const Matrix<DefDiag, T, M> &B, std::size_t matrix_size)
       -> Matrix<DefDense, T, M, M> {
 
-    Base::Matrix::Matrix<T, M, K> X_1;
+    if (matrix_size > M) {
+      matrix_size = M;
+    }
+
+    Base::Matrix::gmres_k_partition_matrix(
+        A.matrix, B.matrix, this->X_1, this->decay_rate, this->division_min,
+        this->rho, this->rep_num, matrix_size);
+
+    return Matrix<DefDense, T, M, M>(this->X_1);
+  }
+
+  inline auto cold_solve(const Matrix<DefDense, T, M, M> &A,
+                         const Matrix<DefDiag, T, M> &B,
+                         std::size_t matrix_size) -> Matrix<DefDense, T, M, M> {
+
+    Base::Matrix::Matrix<T, M, K> X_1_temporary;
 
     if (matrix_size > M) {
       matrix_size = M;
     }
 
     Base::Matrix::gmres_k_partition_matrix(
-        A.matrix, B.matrix, X_1, this->decay_rate, this->division_min,
+        A.matrix, B.matrix, X_1_temporary, this->decay_rate, this->division_min,
         this->rho, this->rep_num, matrix_size);
 
-    return Matrix<DefDense, T, M, M>(X_1);
+    return Matrix<DefDense, T, M, M>(X_1_temporary);
   }
 
   inline auto solve(const Matrix<DefDense, T, M, M> &A,
                     const Matrix<DefSparse, T, M, K, SparseAvailable_B> &B,
                     std::size_t matrix_size) -> Matrix<DefDense, T, M, K> {
 
-    Base::Matrix::Matrix<T, M, K> X_1;
+    if (matrix_size > M) {
+      matrix_size = M;
+    }
+
+    Base::Matrix::Matrix<T, M, K> B_dense_matrix =
+        Base::Matrix::output_dense_matrix(B.matrix);
+
+    Base::Matrix::gmres_k_partition_matrix(
+        A.matrix, B_dense_matrix, this->X_1, this->decay_rate,
+        this->division_min, this->rho, this->rep_num, matrix_size);
+
+    return Matrix<DefDense, T, M, K>(this->X_1);
+  }
+
+  inline auto cold_solve(const Matrix<DefDense, T, M, M> &A,
+                         const Matrix<DefSparse, T, M, K, SparseAvailable_B> &B,
+                         std::size_t matrix_size) -> Matrix<DefDense, T, M, K> {
+
+    Base::Matrix::Matrix<T, M, K> X_1_temporary;
 
     if (matrix_size > M) {
       matrix_size = M;
@@ -718,10 +770,10 @@ class LinalgPartitionSolver {
         Base::Matrix::output_dense_matrix(B.matrix);
 
     Base::Matrix::gmres_k_partition_matrix(
-        A.matrix, B_dense_matrix, X_1, this->decay_rate, this->division_min,
-        this->rho, this->rep_num, matrix_size);
+        A.matrix, B_dense_matrix, X_1_temporary, this->decay_rate,
+        this->division_min, this->rho, this->rep_num, matrix_size);
 
-    return Matrix<DefDense, T, M, K>(X_1);
+    return Matrix<DefDense, T, M, K>(X_1_temporary);
   }
 
   inline auto solve(const Matrix<DefDiag, T, M> &A,
@@ -732,11 +784,25 @@ class LinalgPartitionSolver {
       matrix_size = M;
     }
 
-    Base::Matrix::Matrix<T, M, K> X_1 =
+    this->X_1 = Base::Matrix::diag_inv_multiply_dense_partition(
+        A.matrix, B.matrix, this->division_min, matrix_size);
+
+    return Matrix<DefDense, T, M, K>(this->X_1);
+  }
+
+  inline auto cold_solve(const Matrix<DefDiag, T, M> &A,
+                         const Matrix<DefDense, T, M, K> &B,
+                         std::size_t matrix_size) -> Matrix<DefDense, T, M, K> {
+
+    if (matrix_size > M) {
+      matrix_size = M;
+    }
+
+    Base::Matrix::Matrix<T, M, K> X_1_temporary =
         Base::Matrix::diag_inv_multiply_dense_partition(
             A.matrix, B.matrix, this->division_min, matrix_size);
 
-    return Matrix<DefDense, T, M, K>(X_1);
+    return Matrix<DefDense, T, M, K>(X_1_temporary);
   }
 
   inline auto solve(const Matrix<DefDiag, T, M> &A,
@@ -758,6 +824,21 @@ class LinalgPartitionSolver {
                     const Matrix<DefSparse, T, M, K, SparseAvailable_B> &B,
                     std::size_t matrix_size) -> Matrix<DefDense, T, M, K> {
 
+    if (matrix_size > M) {
+      matrix_size = M;
+    }
+
+    this->X_1 = Base::Matrix::output_dense_matrix(
+        Base::Matrix::diag_inv_multiply_sparse_partition(
+            A.matrix, B.matrix, this->division_min, matrix_size));
+
+    return Matrix<DefDense, T, M, K>(this->X_1);
+  }
+
+  inline auto cold_solve(const Matrix<DefDiag, T, M> &A,
+                         const Matrix<DefSparse, T, M, K, SparseAvailable_B> &B,
+                         std::size_t matrix_size) -> Matrix<DefDense, T, M, K> {
+
     using RowIndices_B = RowIndicesFromSparseAvailable<SparseAvailable_B>;
     using RowPointers_B = RowPointersFromSparseAvailable<SparseAvailable_B>;
 
@@ -766,34 +847,68 @@ class LinalgPartitionSolver {
     }
 
     Base::Matrix::CompiledSparseMatrix<T, M, K, RowIndices_B, RowPointers_B>
-        X_1 = Base::Matrix::diag_inv_multiply_sparse_partition(
+        X_1_temporary = Base::Matrix::diag_inv_multiply_sparse_partition(
             A.matrix, B.matrix, this->division_min, matrix_size);
 
-    return Matrix<DefDense, T, M, K>(Base::Matrix::output_dense_matrix(X_1));
+    return Matrix<DefDense, T, M, K>(
+        Base::Matrix::output_dense_matrix(X_1_temporary));
   }
 
   inline auto solve(const Matrix<DefSparse, T, M, M, SparseAvailable_A> &A,
                     const Matrix<DefDense, T, M, K> &B, std::size_t matrix_size)
       -> Matrix<DefDense, T, M, K> {
 
-    Base::Matrix::Matrix<T, M, K> X_1;
+    if (matrix_size > M) {
+      matrix_size = M;
+    }
+
+    Base::Matrix::sparse_gmres_k_partition_matrix(
+        A.matrix, B.matrix, this->X_1, this->decay_rate, this->division_min,
+        this->rho, this->rep_num, matrix_size);
+
+    return Matrix<DefDense, T, M, K>(this->X_1);
+  }
+
+  inline auto cold_solve(const Matrix<DefSparse, T, M, M, SparseAvailable_A> &A,
+                         const Matrix<DefDense, T, M, K> &B,
+                         std::size_t matrix_size) -> Matrix<DefDense, T, M, K> {
+
+    Base::Matrix::Matrix<T, M, K> X_1_temporary;
 
     if (matrix_size > M) {
       matrix_size = M;
     }
 
     Base::Matrix::sparse_gmres_k_partition_matrix(
-        A.matrix, B.matrix, X_1, this->decay_rate, this->division_min,
+        A.matrix, B.matrix, X_1_temporary, this->decay_rate, this->division_min,
         this->rho, this->rep_num, matrix_size);
 
-    return Matrix<DefDense, T, M, K>(X_1);
+    return Matrix<DefDense, T, M, K>(X_1_temporary);
   }
 
   inline auto solve(const Matrix<DefSparse, T, M, M, SparseAvailable_A> &A,
                     const Matrix<DefDiag, T, M> &B, std::size_t matrix_size)
       -> Matrix<DefDense, T, M, M> {
 
-    Base::Matrix::Matrix<T, M, K> X_1;
+    if (matrix_size > M) {
+      matrix_size = M;
+    }
+
+    Base::Matrix::Matrix<T, M, M> B_dense_matrix =
+        output_dense_matrix(B.matrix);
+
+    Base::Matrix::sparse_gmres_k_partition_matrix(
+        A.matrix, B_dense_matrix, this->X_1, this->decay_rate,
+        this->division_min, this->rho, this->rep_num, matrix_size);
+
+    return Matrix<DefDense, T, M, M>(this->X_1);
+  }
+
+  inline auto cold_solve(const Matrix<DefSparse, T, M, M, SparseAvailable_A> &A,
+                         const Matrix<DefDiag, T, M> &B,
+                         std::size_t matrix_size) -> Matrix<DefDense, T, M, M> {
+
+    Base::Matrix::Matrix<T, M, K> X_1_temporary;
 
     if (matrix_size > M) {
       matrix_size = M;
@@ -803,17 +918,35 @@ class LinalgPartitionSolver {
         output_dense_matrix(B.matrix);
 
     Base::Matrix::sparse_gmres_k_partition_matrix(
-        A.matrix, B_dense_matrix, X_1, this->decay_rate, this->division_min,
-        this->rho, this->rep_num, matrix_size);
+        A.matrix, B_dense_matrix, X_1_temporary, this->decay_rate,
+        this->division_min, this->rho, this->rep_num, matrix_size);
 
-    return Matrix<DefDense, T, M, M>(X_1);
+    return Matrix<DefDense, T, M, M>(X_1_temporary);
   }
 
   inline auto solve(const Matrix<DefSparse, T, M, M, SparseAvailable_A> &A,
                     const Matrix<DefSparse, T, M, K, SparseAvailable_B> &B,
                     std::size_t matrix_size) -> Matrix<DefDense, T, M, K> {
 
-    Base::Matrix::Matrix<T, M, K> X_1;
+    if (matrix_size > M) {
+      matrix_size = M;
+    }
+
+    Base::Matrix::Matrix<T, M, K> B_dense_matrix =
+        Base::Matrix::output_dense_matrix(B.matrix);
+
+    Base::Matrix::sparse_gmres_k_partition_matrix(
+        A.matrix, B_dense_matrix, this->X_1, this->decay_rate,
+        this->division_min, this->rho, this->rep_num, matrix_size);
+
+    return Matrix<DefDense, T, M, K>(this->X_1);
+  }
+
+  inline auto cold_solve(const Matrix<DefSparse, T, M, M, SparseAvailable_A> &A,
+                         const Matrix<DefSparse, T, M, K, SparseAvailable_B> &B,
+                         std::size_t matrix_size) -> Matrix<DefDense, T, M, K> {
+
+    Base::Matrix::Matrix<T, M, K> X_1_temporary;
 
     if (matrix_size > M) {
       matrix_size = M;
@@ -823,10 +956,10 @@ class LinalgPartitionSolver {
         Base::Matrix::output_dense_matrix(B.matrix);
 
     Base::Matrix::sparse_gmres_k_partition_matrix(
-        A.matrix, B_dense_matrix, X_1, this->decay_rate, this->division_min,
-        this->rho, this->rep_num, matrix_size);
+        A.matrix, B_dense_matrix, X_1_temporary, this->decay_rate,
+        this->division_min, this->rho, this->rep_num, matrix_size);
 
-    return Matrix<DefDense, T, M, K>(X_1);
+    return Matrix<DefDense, T, M, K>(X_1_temporary);
   }
 
   inline auto get_answer(void) -> Matrix<DefDense, T, M, K> {
@@ -850,6 +983,8 @@ class LinalgPartitionSolver {
 
 public:
   /* Variable */
+  Base::Matrix::Matrix<T, M, K> X_1;
+
   Value_Type decay_rate;
   Value_Type division_min;
   std::array<Value_Type, K> rho;
diff --git a/test_vs/check_python_numpy.hpp b/test_vs/check_python_numpy.hpp
index 9eb65d7..8b57539 100644
--- a/test_vs/check_python_numpy.hpp
+++ b/test_vs/check_python_numpy.hpp
@@ -1217,6 +1217,11 @@ void CheckPythonNumpy<T>::check_python_numpy_left_divide_and_inv(void) {
     tester.expect_near(A_A_P_x.matrix.data, A_A_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
         "check LinalgPartitionSolver solve Dense and Dense.");
 
+    A_A_P_x = A_A_P_linalg_solver.cold_solve(A_P, A_P, 3);
+
+    tester.expect_near(A_A_P_x.matrix.data, A_A_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
+        "check LinalgPartitionSolver solve Dense and Dense cold solve.");
+
     static auto A_B_P_linalg_solver = make_LinalgPartitionSolver<decltype(A_P), decltype(B_P)>();
 
     auto A_B_P_x = A_B_P_linalg_solver.solve(A_P, B_P, 3);
@@ -1231,6 +1236,11 @@ void CheckPythonNumpy<T>::check_python_numpy_left_divide_and_inv(void) {
     tester.expect_near(A_B_P_x.matrix.data, A_B_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
         "check LinalgPartitionSolver solve Dense and Diag.");
 
+    A_B_P_x = A_B_P_linalg_solver.cold_solve(A_P, B_P, 3);
+
+    tester.expect_near(A_B_P_x.matrix.data, A_B_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
+        "check LinalgPartitionSolver solve Dense and Diag cold solve.");
+
     static auto A_C_P_linalg_solver = make_LinalgPartitionSolver<decltype(A_P), decltype(C_P)>();
 
     auto A_C_P_x = A_C_P_linalg_solver.solve(A_P, C_P, 3);
@@ -1245,6 +1255,11 @@ void CheckPythonNumpy<T>::check_python_numpy_left_divide_and_inv(void) {
     tester.expect_near(A_C_P_x.matrix.data, A_C_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
         "check LinalgPartitionSolver solve Dense and Sparse.");
 
+    A_C_P_x = A_C_P_linalg_solver.cold_solve(A_P, C_P, 3);
+
+    tester.expect_near(A_C_P_x.matrix.data, A_C_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
+        "check LinalgPartitionSolver solve Dense and Sparse cold solve.");
+
     static auto B_A_P_linalg_solver = make_LinalgPartitionSolver<decltype(B_P), decltype(A_P)>();
 
     auto B_A_P_x = B_A_P_linalg_solver.solve(B_P, A_P, 3);
@@ -1259,6 +1274,11 @@ void CheckPythonNumpy<T>::check_python_numpy_left_divide_and_inv(void) {
     tester.expect_near(B_A_P_x.matrix.data, B_A_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
         "check LinalgPartitionSolver solve Diag and Dense.");
 
+    B_A_P_x = B_A_P_linalg_solver.cold_solve(B_P, A_P, 3);
+
+    tester.expect_near(B_A_P_x.matrix.data, B_A_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
+        "check LinalgPartitionSolver solve Diag and Dense cold solve.");
+
     static auto B_B_P_linalg_solver = make_LinalgPartitionSolver<decltype(B_P), decltype(B_P)>();
 
     auto B_B_P_x = B_B_P_linalg_solver.solve(B_P, B_P, 3);
@@ -1288,6 +1308,11 @@ void CheckPythonNumpy<T>::check_python_numpy_left_divide_and_inv(void) {
     tester.expect_near(B_C_P_x.matrix.data, B_C_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
         "check LinalgPartitionSolver solve Diag and Sparse.");
 
+    B_C_P_x = B_C_P_linalg_solver.cold_solve(B_P, C_P, 3);
+
+    tester.expect_near(B_C_P_x.matrix.data, B_C_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
+        "check LinalgPartitionSolver solve Diag and Sparse cold solve.");
+
     static auto C_A_P_linalg_solver = make_LinalgPartitionSolver<decltype(C_P), decltype(A_P)>();
 
     auto C_A_P_x = C_A_P_linalg_solver.solve(C_P, A_P, 3);
@@ -1302,6 +1327,11 @@ void CheckPythonNumpy<T>::check_python_numpy_left_divide_and_inv(void) {
     tester.expect_near(C_A_P_x.matrix.data, C_A_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
         "check LinalgPartitionSolver solve Sparse and Dense.");
 
+    C_A_P_x = C_A_P_linalg_solver.cold_solve(C_P, A_P, 3);
+
+    tester.expect_near(C_A_P_x.matrix.data, C_A_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
+        "check LinalgPartitionSolver solve Sparse and Dense cold solve.");
+
     static auto C_B_P_linalg_solver = make_LinalgPartitionSolver<decltype(C_P), decltype(B_P)>();
 
     auto C_B_P_x = C_B_P_linalg_solver.solve(C_P, B_P, 3);
@@ -1316,6 +1346,11 @@ void CheckPythonNumpy<T>::check_python_numpy_left_divide_and_inv(void) {
     tester.expect_near(C_B_P_x.matrix.data, C_B_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
         "check LinalgPartitionSolver solve Sparse and Diag.");
 
+    C_B_P_x = C_B_P_linalg_solver.cold_solve(C_P, B_P, 3);
+
+    tester.expect_near(C_B_P_x.matrix.data, C_B_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
+        "check LinalgPartitionSolver solve Sparse and Diag cold solve.");
+
     static auto C_C_P_linalg_solver = make_LinalgPartitionSolver<decltype(C_P), decltype(C_P)>();
 
     auto C_C_P_x = C_C_P_linalg_solver.solve(C_P, C_P, 3);
@@ -1330,6 +1365,11 @@ void CheckPythonNumpy<T>::check_python_numpy_left_divide_and_inv(void) {
     tester.expect_near(C_C_P_x.matrix.data, C_C_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
         "check LinalgPartitionSolver solve Sparse and Sparse.");
 
+    C_C_P_x = C_C_P_linalg_solver.cold_solve(C_P, C_P, 3);
+
+    tester.expect_near(C_C_P_x.matrix.data, C_C_P_x_answer.matrix.data, NEAR_LIMIT_STRICT,
+        "check LinalgPartitionSolver solve Sparse and Sparse cold solve.");
+
 
     /* 矩形　左除算 */
     Matrix<DefDense, T, 4, 3> AL({ {1, 2, 3}, {5, 4, 6}, {9, 8, 7}, {2, 2, 3} });