diff --git a/rust/Cargo.lock b/rust/Cargo.lock
index 20af5cf..c5cc059 100644
--- a/rust/Cargo.lock
+++ b/rust/Cargo.lock
@@ -5,3 +5,57 @@ version = 4
 [[package]]
 name = "NeuralNetwork"
 version = "0.1.0"
+dependencies = [
+ "rayon",
+]
+
+[[package]]
+name = "crossbeam-deque"
+version = "0.8.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51"
+dependencies = [
+ "crossbeam-epoch",
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-epoch"
+version = "0.9.18"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e"
+dependencies = [
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-utils"
+version = "0.8.21"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28"
+
+[[package]]
+name = "either"
+version = "1.16.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "91622ff5e7162018101f2fea40d6ebf4a78bbe5a49736a2020649edf9693679e"
+
+[[package]]
+name = "rayon"
+version = "1.12.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fb39b166781f92d482534ef4b4b1b2568f42613b53e5b6c160e24cfbfa30926d"
+dependencies = [
+ "either",
+ "rayon-core",
+]
+
+[[package]]
+name = "rayon-core"
+version = "1.13.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91"
+dependencies = [
+ "crossbeam-deque",
+ "crossbeam-utils",
+]
diff --git a/rust/Cargo.toml b/rust/Cargo.toml
index 8766a95..eb08f2e 100644
--- a/rust/Cargo.toml
+++ b/rust/Cargo.toml
@@ -4,3 +4,4 @@ version = "0.1.0"
 edition = "2024"
 
 [dependencies]
+rayon = "1.12.0"
diff --git a/rust/docs/matrix/add.md b/rust/docs/matrix/add.md
index 689661a..533f445 100644
--- a/rust/docs/matrix/add.md
+++ b/rust/docs/matrix/add.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの加算を行います。
 * サイズが異なる行列同士の加算はエラーを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/add_assign.md b/rust/docs/matrix/add_assign.md
index 9495337..c18048b 100644
--- a/rust/docs/matrix/add_assign.md
+++ b/rust/docs/matrix/add_assign.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの加算を行います。
 * サイズが異なる行列同士の加算はパニックを引き起こします。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/add_trait.md b/rust/docs/matrix/add_trait.md
index 29937e2..122bc24 100644
--- a/rust/docs/matrix/add_trait.md
+++ b/rust/docs/matrix/add_trait.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの加算を行います。
 * サイズが異なる行列同士の加算はNoneを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/add_with.md b/rust/docs/matrix/add_with.md
index 2962f80..2a35e17 100644
--- a/rust/docs/matrix/add_with.md
+++ b/rust/docs/matrix/add_with.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの加算を行います。
 * サイズが異なる行列同士の加算はエラーを返します。
+
 # Examples
-``` 
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/col_iter.md b/rust/docs/matrix/col_iter.md
index 025c591..0df5566 100644
--- a/rust/docs/matrix/col_iter.md
+++ b/rust/docs/matrix/col_iter.md
@@ -1,6 +1,8 @@
 指定した列の要素への参照を返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m: Matrix<i32> = Matrix::with_size(2, 2);
 assert_eq!(m.col_iter(0).unwrap().collect::<Vec<&i32>>(), vec![&0, &0]);
diff --git a/rust/docs/matrix/col_mut_iter.md b/rust/docs/matrix/col_iter_mut.md
similarity index 82%
rename from rust/docs/matrix/col_mut_iter.md
rename to rust/docs/matrix/col_iter_mut.md
index 3bffebd..4446070 100644
--- a/rust/docs/matrix/col_mut_iter.md
+++ b/rust/docs/matrix/col_iter_mut.md
@@ -1,9 +1,11 @@
 指定した列の要素への可変参照を返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m: Matrix<i32> = Matrix::with_size(2, 2);
-if let Some(col_iter) = m.col_mut_iter(0) {
+if let Some(col_iter) = m.col_iter_mut(0) {
    col_iter.for_each(|val| *val = 42);
 }
 assert_eq!(m.col_iter(0).unwrap().collect::<Vec<&i32>>(), vec![&42, &42]);
diff --git a/rust/docs/matrix/col_par_iter.md b/rust/docs/matrix/col_par_iter.md
new file mode 100644
index 0000000..357c926
--- /dev/null
+++ b/rust/docs/matrix/col_par_iter.md
@@ -0,0 +1,13 @@
+指定した列の要素への並列参照を返します。
+このイテレータは、列の要素への可変参照を返します。
+
+# Examples
+
+```rust
+use matrix::Matrix;
+let mut m: Matrix<i32> = Matrix::with_size(2, 2);
+if let Some(col_iter) = m.col_par_iter_mut(0) {
+   col_iter.for_each(|val| *val = 42);
+}
+assert_eq!(m.col_iter(0).unwrap().collect::<Vec<&i32>>(), vec![&42, &42]);
+```
\ No newline at end of file
diff --git a/rust/docs/matrix/col_par_iter_mut.md b/rust/docs/matrix/col_par_iter_mut.md
new file mode 100644
index 0000000..a57d106
--- /dev/null
+++ b/rust/docs/matrix/col_par_iter_mut.md
@@ -0,0 +1,12 @@
+指定した列の要素への並列可変参照を返します。
+このイテレータは、列の要素への可変参照を返します。イテレータは、列の要素を順番に返しますが、同時に複数の要素にアクセスすることができます。
+
+# Examples
+
+```rust
+use matrix::Matrix;
+let mut m: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
+m.col_par_iter_mut(0).unwrap().for_each(|x| *x *= 2);
+assert_eq!(m[(0, 0)], 2);
+assert_eq!(m[(1, 0)], 6);
+```
\ No newline at end of file
diff --git a/rust/docs/matrix/cols.md b/rust/docs/matrix/cols.md
index 791bf2e..ff2c20c 100644
--- a/rust/docs/matrix/cols.md
+++ b/rust/docs/matrix/cols.md
@@ -1,6 +1,8 @@
 列数を返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m: Matrix<i32> = Matrix::with_size(3, 4); // [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
 assert_eq!(m.cols(), 4);
diff --git a/rust/docs/matrix/div.md b/rust/docs/matrix/div.md
index 99c4f68..8158768 100644
--- a/rust/docs/matrix/div.md
+++ b/rust/docs/matrix/div.md
@@ -1,8 +1,10 @@
 同じサイズの行列同士の要素ごとの除算を行います。
 * サイズが異なる行列同士の除算はエラーを返します。
 * 0で割る要素がある場合もエラーを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/div_assign.md b/rust/docs/matrix/div_assign.md
index 4990ab1..35f124b 100644
--- a/rust/docs/matrix/div_assign.md
+++ b/rust/docs/matrix/div_assign.md
@@ -1,8 +1,10 @@
 同じサイズの行列同士の要素ごとの除算を行います。
 * サイズが異なる行列同士の除算はパニックを引き起こします。
 * 0で割る要素がある場合もパニックを引き起こします。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[10, 20], [30, 40]]);
 let m2: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
diff --git a/rust/docs/matrix/div_trait.md b/rust/docs/matrix/div_trait.md
index ed74e37..64f6287 100644
--- a/rust/docs/matrix/div_trait.md
+++ b/rust/docs/matrix/div_trait.md
@@ -1,8 +1,10 @@
 同じサイズの行列同士の要素ごとの除算を行います。
 * サイズが異なる行列同士の除算はNoneを返します。
 * 0で割る要素がある場合もNoneを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m2: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m3 = m1 / m2;
diff --git a/rust/docs/matrix/div_with.md b/rust/docs/matrix/div_with.md
index 5cb2fb1..65a936b 100644
--- a/rust/docs/matrix/div_with.md
+++ b/rust/docs/matrix/div_with.md
@@ -1,8 +1,10 @@
 同じサイズの行列同士の要素ごとの除算を行います。
 * サイズが異なる行列同士の除算はエラーを返します。
 * 0で割る要素がある場合もエラーを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/get.md b/rust/docs/matrix/get.md
index 012747b..90061e3 100644
--- a/rust/docs/matrix/get.md
+++ b/rust/docs/matrix/get.md
@@ -1,6 +1,8 @@
 指定した位置の要素への参照を返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m: Matrix<i32> = Matrix::with_size(2, 2);
 assert_eq!(m.get(0, 0), Some(&0));
diff --git a/rust/docs/matrix/get_mut.md b/rust/docs/matrix/get_mut.md
index 34f36da..3b2134a 100644
--- a/rust/docs/matrix/get_mut.md
+++ b/rust/docs/matrix/get_mut.md
@@ -1,6 +1,8 @@
 指定した位置の要素への可変参照を返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m: Matrix<i32> = Matrix::with_size(2, 2);
 if let Some(val) = m.get_mut(0, 0) {
diff --git a/rust/docs/matrix/hadamard_mul.md b/rust/docs/matrix/hadamard_mul.md
index 2ef9f65..53b80e2 100644
--- a/rust/docs/matrix/hadamard_mul.md
+++ b/rust/docs/matrix/hadamard_mul.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの乗算を行います。
 * サイズが異なる行列同士の乗算はエラーを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/hadamard_mul_with.md b/rust/docs/matrix/hadamard_mul_with.md
index 751a6c8..58311ab 100644
--- a/rust/docs/matrix/hadamard_mul_with.md
+++ b/rust/docs/matrix/hadamard_mul_with.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの乗算を行います。
 * サイズが異なる行列同士の乗算はエラーを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/index.md b/rust/docs/matrix/index.md
index d54f408..787ab31 100644
--- a/rust/docs/matrix/index.md
+++ b/rust/docs/matrix/index.md
@@ -1,6 +1,8 @@
 タプルインデックスで要素にアクセスします。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 assert_eq!(m[(0, 0)], 1);
diff --git a/rust/docs/matrix/mul.md b/rust/docs/matrix/mul.md
index 07c2854..0fd6728 100644
--- a/rust/docs/matrix/mul.md
+++ b/rust/docs/matrix/mul.md
@@ -1,7 +1,9 @@
 行列積を行います。
 * 列数と行数が一致しない行列同士の乗算はエラーを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/mul_assign.md b/rust/docs/matrix/mul_assign.md
index 208201f..e970eb4 100644
--- a/rust/docs/matrix/mul_assign.md
+++ b/rust/docs/matrix/mul_assign.md
@@ -1,7 +1,9 @@
 行列積を行います。
 * サイズが異なる行列同士の乗算はパニックを引き起こします。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/mul_trait.md b/rust/docs/matrix/mul_trait.md
index 3e11e85..a161fa6 100644
--- a/rust/docs/matrix/mul_trait.md
+++ b/rust/docs/matrix/mul_trait.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの乗算を行います。
 * サイズが異なる行列同士の乗算はNoneを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/mul_with.md b/rust/docs/matrix/mul_with.md
index 4556afb..f50437b 100644
--- a/rust/docs/matrix/mul_with.md
+++ b/rust/docs/matrix/mul_with.md
@@ -1,7 +1,9 @@
 行列積を行います。
 * 列数と行数が一致しない行列同士の乗算はエラーを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/new.md b/rust/docs/matrix/new.md
index 3259a15..80d7e4e 100644
--- a/rust/docs/matrix/new.md
+++ b/rust/docs/matrix/new.md
@@ -1,7 +1,11 @@
 二次元配列から行列を作成します。
 * TはDefaultとCopyトレイトを実装している必要があります。
+
 # Examples
+
+```rust
 use NeuralNetwork::matrix::Matrix;
+let m: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 assert_eq!(m.get(0, 0), Some(&1));
 assert_eq!(m.get(0, 1), Some(&2));
 assert_eq!(m.get(1, 0), Some(&3));
diff --git a/rust/docs/matrix/row_iter.md b/rust/docs/matrix/row_iter.md
index b372464..e64f6c0 100644
--- a/rust/docs/matrix/row_iter.md
+++ b/rust/docs/matrix/row_iter.md
@@ -1,6 +1,8 @@
 指定した行の要素への参照を返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m: Matrix<i32> = Matrix::with_size(2, 2);
 assert_eq!(m.row_iter(0).unwrap().collect::<Vec<&i32>>(), vec![&0, &0]);
diff --git a/rust/docs/matrix/row_mut_iter.md b/rust/docs/matrix/row_iter_mut.md
similarity index 96%
rename from rust/docs/matrix/row_mut_iter.md
rename to rust/docs/matrix/row_iter_mut.md
index 4f4111e..8f4b252 100644
--- a/rust/docs/matrix/row_mut_iter.md
+++ b/rust/docs/matrix/row_iter_mut.md
@@ -1,6 +1,8 @@
 指定した行の要素への可変参照を返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m: Matrix<i32> = Matrix::with_size(2, 2);
 if let Some(row_iter) = m.row_mut_iter(0) {
diff --git a/rust/docs/matrix/row_par_iter.md b/rust/docs/matrix/row_par_iter.md
new file mode 100644
index 0000000..2b0879f
--- /dev/null
+++ b/rust/docs/matrix/row_par_iter.md
@@ -0,0 +1,11 @@
+並列イテレータを返します。
+このメソッドは、行の要素を並列に処理したい場合に使用します。
+
+# Examples
+
+```rust
+use matrix::Matrix;
+let m: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
+let sum: i32 = m.row_par_iter(0).unwrap().cloned().sum();
+assert_eq!(sum, 3);
+```
\ No newline at end of file
diff --git a/rust/docs/matrix/row_par_iter_mut.md b/rust/docs/matrix/row_par_iter_mut.md
new file mode 100644
index 0000000..00b53af
--- /dev/null
+++ b/rust/docs/matrix/row_par_iter_mut.md
@@ -0,0 +1,11 @@
+並列イテレータを返します。
+このメソッドは、行の要素を並列に処理したい場合に使用します。
+
+# Examples
+
+```rust
+use matrix::Matrix;
+let m: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
+m.row_par_iter_mut(0).unwrap().for_each(|x| *x *= 2);
+assert_eq!(m[(0, 0)], 2);
+```
\ No newline at end of file
diff --git a/rust/docs/matrix/rows.md b/rust/docs/matrix/rows.md
index 11fee2f..5041f53 100644
--- a/rust/docs/matrix/rows.md
+++ b/rust/docs/matrix/rows.md
@@ -1,6 +1,8 @@
 行数を返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m: Matrix<i32> = Matrix::with_size(3, 4); // [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
 assert_eq!(m.rows(), 3);
diff --git a/rust/docs/matrix/sub.md b/rust/docs/matrix/sub.md
index 5d00383..64254a6 100644
--- a/rust/docs/matrix/sub.md
+++ b/rust/docs/matrix/sub.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの減算を行います。
 * サイズが異なる行列同士の減算はエラーを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/sub_assign.md b/rust/docs/matrix/sub_assign.md
index 58c934a..c671e33 100644
--- a/rust/docs/matrix/sub_assign.md
+++ b/rust/docs/matrix/sub_assign.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの減算を行います。
 * サイズが異なる行列同士の減算はパニックを引き起こします。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/sub_trait.md b/rust/docs/matrix/sub_trait.md
index 64b40ea..6693bf1 100644
--- a/rust/docs/matrix/sub_trait.md
+++ b/rust/docs/matrix/sub_trait.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの減算を行います。
 * サイズが異なる行列同士の減算はNoneを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/sub_with.md b/rust/docs/matrix/sub_with.md
index 145b1b4..45c67e9 100644
--- a/rust/docs/matrix/sub_with.md
+++ b/rust/docs/matrix/sub_with.md
@@ -1,7 +1,9 @@
 同じサイズの行列同士の要素ごとの減算を行います。
 * サイズが異なる行列同士の減算はエラーを返します。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let mut m1: Matrix<i32> = Matrix::new([[1, 2], [3, 4]]);
 let m2: Matrix<i32> = Matrix::new([[5, 6], [7, 8]]);
diff --git a/rust/docs/matrix/with_size.md b/rust/docs/matrix/with_size.md
index 3aaacbb..607f6ec 100644
--- a/rust/docs/matrix/with_size.md
+++ b/rust/docs/matrix/with_size.md
@@ -1,8 +1,10 @@
 サイズを指定して行列を作成します。
 * 要素は全てTのデフォルト値で初期化されます。
 * TはDefaultとCopyトレイトを実装している必要があります。
+
 # Examples
-```
+
+```rust
 use matrix::Matrix;
 let m: Matrix<i32> = Matrix::with_size(2, 2);
 assert_eq!(m.get(0, 0), Some(&0));
diff --git a/rust/src/lib.rs b/rust/src/lib.rs
index 1c13d99..cfbe12a 100644
--- a/rust/src/lib.rs
+++ b/rust/src/lib.rs
@@ -1,3 +1,4 @@
+pub mod matrix_element;
 pub mod matrix_layout;
 pub mod matrix_algorithm;
 pub mod matrix;
\ No newline at end of file
diff --git a/rust/src/matrix.rs b/rust/src/matrix.rs
index 3f17c59..dd26b9f 100644
--- a/rust/src/matrix.rs
+++ b/rust/src/matrix.rs
@@ -1,6 +1,7 @@
 use core::fmt;
 use std::{marker::PhantomData, ops::Mul};
-use crate::{matrix_algorithm::{MatrixAlgorithm, NaiveAlgorithm}, matrix_layout::{ MatrixLayout, RowMajor }};
+use rayon::prelude::*;
+use crate::{matrix_algorithm::{MatrixAlgorithm, NaiveAlgorithm}, matrix_element::MatrixElement, matrix_layout::{ MatrixLayout, RowMajor }};
 
 pub struct Matrix<T, L: MatrixLayout = RowMajor> {
     data    : Vec<T>,
@@ -65,28 +66,62 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
         Some(self.data.iter().skip(start).step_by(step).take(self.cols))
     }
 
-    #[doc = include_str!("../docs/matrix/row_mut_iter.md")]
-    pub fn row_mut_iter(&mut self, row: usize) -> Option<impl Iterator<Item = &mut T>> {
+    #[doc = include_str!("../docs/matrix/row_par_iter.md")]
+    pub fn row_par_iter(&self, row: usize) -> Option<impl IndexedParallelIterator<Item = &T>>
+    where T: Sync
+    {
+        let (start, step) = L::row_stride(row, self.rows, self.cols)?;
+        Some(self.data.par_iter().skip(start).step_by(step).take(self.cols))
+    }
+
+    #[doc = include_str!("../docs/matrix/row_iter_mut.md")]
+    pub fn row_iter_mut(&mut self, row: usize) -> Option<impl Iterator<Item = &mut T>> {
         let (start, step) = L::row_stride(row, self.rows, self.cols)?;
         Some(self.data.iter_mut().skip(start).step_by(step).take(self.cols))
     }
 
+    #[doc = include_str!("../docs/matrix/row_par_iter_mut.md")]
+    pub fn row_par_iter_mut(&mut self, row: usize) -> Option<impl IndexedParallelIterator<Item = &mut T>>
+    where T: Sync + Send
+    {
+        let (start, step) = L::row_stride(row, self.rows, self.cols)?;
+        Some(self.data.par_iter_mut().skip(start).step_by(step).take(self.cols))
+    }
+
     #[doc = include_str!("../docs/matrix/col_iter.md")]
     pub fn col_iter(&self, col: usize) -> Option<impl Iterator<Item = &T>> {
         let (start, step) = L::col_stride(col, self.rows, self.cols)?;
         Some(self.data.iter().skip(start).step_by(step).take(self.rows))
     }
 
-    #[doc = include_str!("../docs/matrix/col_mut_iter.md")]
-    pub fn col_mut_iter(&mut self, col: usize) -> Option<impl Iterator<Item = &mut T>> {
+    #[doc = include_str!("../docs/matrix/col_par_iter.md")]
+    pub fn col_par_iter(&self, col: usize) -> Option<impl IndexedParallelIterator<Item = &T>>
+    where 
+        T: Sync 
+    {
+        let (start, step) = L::col_stride(col, self.rows, self.cols)?;
+        Some(self.data.par_iter().skip(start).step_by(step).take(self.rows))
+    }
+
+    #[doc = include_str!("../docs/matrix/col_iter_mut.md")]
+    pub fn col_iter_mut(&mut self, col: usize) -> Option<impl Iterator<Item = &mut T>> {
         let (start, step) = L::col_stride(col, self.rows, self.cols)?;
         Some(self.data.iter_mut().skip(start).step_by(step).take(self.rows))
     }
 
+    #[doc = include_str!("../docs/matrix/col_par_iter_mut.md")]
+    pub fn col_par_iter_mut(&mut self, col: usize) -> Option<impl IndexedParallelIterator<Item = &mut T>>
+    where
+        T: Sync + Send,
+    {
+        let (start, step) = L::col_stride(col, self.rows, self.cols)?;
+        Some(self.data.par_iter_mut().skip(start).step_by(step).take(self.rows))
+    }
+
     #[doc = include_str!("../docs/matrix/add.md")]
     pub fn add(&mut self, other: &Matrix<T, L>) -> Result<(), String> 
     where 
-        T: std::ops::AddAssign + Copy 
+        T: std::ops::AddAssign + MatrixElement
     {
         NaiveAlgorithm::add(self, other)
     }
@@ -94,7 +129,7 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
     #[doc = include_str!("../docs/matrix/add_with.md")]
     pub fn add_with<Calc: MatrixAlgorithm<T, L, OL>, OL: MatrixLayout>(&mut self, other: &Matrix<T, OL>) -> Result<(), String> 
     where 
-        T: std::ops::AddAssign + Copy 
+        T: std::ops::AddAssign + MatrixElement
     {
         Calc::add(self, other)
     }
@@ -102,7 +137,7 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
     #[doc = include_str!("../docs/matrix/sub.md")]
     pub fn sub(&mut self, other: &Matrix<T, L>) -> Result<(), String> 
     where 
-        T: std::ops::SubAssign + Copy 
+        T: std::ops::SubAssign + MatrixElement
     {
         NaiveAlgorithm::sub(self, other)
     }
@@ -110,7 +145,7 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
     #[doc = include_str!("../docs/matrix/sub_with.md")]
     pub fn sub_with<Calc: MatrixAlgorithm<T, L, OL>, OL: MatrixLayout>(&mut self, other: &Matrix<T, OL>) -> Result<(), String> 
     where 
-        T: std::ops::SubAssign + Copy 
+        T: std::ops::SubAssign + MatrixElement
     {
         Calc::sub(self, other)
     }
@@ -118,7 +153,7 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
     #[doc = include_str!("../docs/matrix/mul.md")]
     pub fn mul(&mut self, other: &Matrix<T, L>) -> Result<(), String> 
     where 
-        T: std::ops::Mul<Output = T> + std::ops::AddAssign + Default + Copy 
+        T: std::ops::Mul<Output = T> + std::ops::AddAssign + MatrixElement + std::iter::Sum<T>
     {
         NaiveAlgorithm::mtx_mul(self, other).map(|result| *self = result)
     }
@@ -126,7 +161,7 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
     #[doc = include_str!("../docs/matrix/mul_with.md")]
     pub fn mul_with<Calc: MatrixAlgorithm<T, L, OL>, OL: MatrixLayout>(&mut self, other: &Matrix<T, OL>) -> Result<(), String> 
     where 
-        T: std::ops::Mul<Output = T> + std::ops::AddAssign + Default + Copy 
+        T: std::ops::Mul<Output = T> + std::ops::AddAssign + MatrixElement + std::iter::Sum<T>
     {
         Calc::mtx_mul(self, other).map(|result| *self = result)
     }
@@ -134,7 +169,7 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
     #[doc = include_str!("../docs/matrix/hadamard_mul.md")]
     pub fn hadamard_mul(&mut self, other: &Matrix<T, L>) -> Result<(), String> 
     where 
-        T: std::ops::MulAssign + Copy 
+        T: std::ops::MulAssign + MatrixElement
     {
         NaiveAlgorithm::hadamard_mul(self, other)
     }
@@ -142,7 +177,7 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
     #[doc = include_str!("../docs/matrix/hadamard_mul_with.md")]
     pub fn hadamard_mul_with<Calc: MatrixAlgorithm<T, L, OL>, OL: MatrixLayout>(&mut self, other: &Matrix<T, OL>) -> Result<(), String> 
     where 
-        T: std::ops::MulAssign + Copy 
+        T: std::ops::MulAssign + MatrixElement
     {
         Calc::hadamard_mul(self, other)
     }
@@ -150,7 +185,7 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
     #[doc = include_str!("../docs/matrix/div.md")]
     pub fn div(&mut self, other: &Matrix<T, L>) -> Result<(), String> 
     where 
-        T: std::ops::DivAssign + Copy 
+        T: std::ops::DivAssign + MatrixElement
     {
         NaiveAlgorithm::div(self, other)
     }
@@ -158,7 +193,7 @@ impl<T, L: MatrixLayout> Matrix<T, L> {
     #[doc = include_str!("../docs/matrix/div_with.md")]
     pub fn div_with<Calc: MatrixAlgorithm<T, L, OL>, OL: MatrixLayout>(&mut self, other: &Matrix<T, OL>) -> Result<(), String> 
     where 
-        T: std::ops::DivAssign + Copy 
+        T: std::ops::DivAssign + MatrixElement
     {
         Calc::div(self, other)
     }
@@ -176,7 +211,7 @@ impl<T, L: MatrixLayout> std::ops::Index<(usize, usize)> for Matrix<T, L> {
 
 impl<T, L: MatrixLayout, OL: MatrixLayout> std::ops::Add<Matrix<T, OL>> for Matrix<T, L> 
 where  
-    T: std::ops::AddAssign + Copy
+    T: std::ops::AddAssign + MatrixElement
 {
     type Output = Option<Matrix<T, L>>;
 
@@ -190,7 +225,7 @@ where
 }
 impl<T, L: MatrixLayout, OL: MatrixLayout> std::ops::AddAssign<Matrix<T, OL>> for Matrix<T, L> 
 where  
-    T: std::ops::AddAssign + Copy
+    T: std::ops::AddAssign + MatrixElement
 {
     #[doc = include_str!("../docs/matrix/add_assign.md")]
     fn add_assign(&mut self, other: Matrix<T, OL>){
@@ -201,7 +236,7 @@ where
 }
 impl<T, L: MatrixLayout, OL: MatrixLayout> std::ops::Sub<Matrix<T, OL>> for Matrix<T, L> 
 where 
-    T: std::ops::SubAssign + Copy
+    T: std::ops::SubAssign + MatrixElement
 {
     type Output = Option<Matrix<T, L>>;
 
@@ -215,7 +250,7 @@ where
 }
 impl<T, L: MatrixLayout, OL: MatrixLayout> std::ops::SubAssign<Matrix<T, OL>> for Matrix<T, L> 
 where  
-    T: std::ops::SubAssign + Copy
+    T: std::ops::SubAssign + MatrixElement
 {
     #[doc = include_str!("../docs/matrix/sub_assign.md")]
     fn sub_assign(&mut self, other: Matrix<T, OL>){
@@ -226,7 +261,7 @@ where
 }
 impl<T, L: MatrixLayout, OL: MatrixLayout> std::ops::Mul<Matrix<T, OL>> for Matrix<T, L> 
 where 
-    T: std::ops::MulAssign + Copy + PartialEq + Default + std::ops::Mul<Output = T> + std::ops::AddAssign
+    T: std::ops::MulAssign + MatrixElement + PartialEq + std::ops::Mul<Output = T> + std::ops::AddAssign + std::iter::Sum<T>
 {
     type Output = Option<Matrix<T, L>>;
 
@@ -243,7 +278,7 @@ where
 }
 impl<T, L: MatrixLayout, OL: MatrixLayout> std::ops::MulAssign<Matrix<T, OL>> for Matrix<T, L> 
 where  
-    T: std::ops::MulAssign + Copy + PartialEq + Default + std::ops::Mul<Output = T> + std::ops::AddAssign
+    T: std::ops::MulAssign + MatrixElement + PartialEq + std::ops::Mul<Output = T> + std::ops::AddAssign + std::iter::Sum<T>
 {
     #[doc = include_str!("../docs/matrix/mul_assign.md")]
     fn mul_assign(&mut self, other: Matrix<T, OL>){
@@ -256,7 +291,7 @@ where
 }
 impl<T, L: MatrixLayout, OL: MatrixLayout> std::ops::Div<Matrix<T, OL>> for Matrix<T, L> 
 where 
-    T: std::ops::DivAssign + Copy + PartialEq + Default
+    T: std::ops::DivAssign + MatrixElement + PartialEq
 {
     type Output = Option<Matrix<T, L>>;
 
@@ -272,7 +307,7 @@ where
 }
 impl<T, L: MatrixLayout, OL: MatrixLayout> std::ops::DivAssign<Matrix<T, OL>> for Matrix<T, L> 
 where  
-    T: std::ops::DivAssign + Copy + PartialEq + Default
+    T: std::ops::DivAssign + MatrixElement + PartialEq
 {
     #[doc = include_str!("../docs/matrix/div_assign.md")]
     fn div_assign(&mut self, other: Matrix<T, OL>){
diff --git a/rust/src/matrix_algorithm.rs b/rust/src/matrix_algorithm.rs
index 7698b65..8071bba 100644
--- a/rust/src/matrix_algorithm.rs
+++ b/rust/src/matrix_algorithm.rs
@@ -1,70 +1,72 @@
 use crate::matrix_layout::MatrixLayout;
 use crate::matrix::Matrix;
+use crate::matrix_element::MatrixElement;
+use rayon::prelude::*;
 
 pub trait MatrixAlgorithm<T, L: MatrixLayout, LO: MatrixLayout> {
-    fn add(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> where T: std::ops::AddAssign + Copy;
-    fn sub(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> where T: std::ops::SubAssign + Copy;
-    fn hadamard_mul(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> where T: std::ops::MulAssign + Copy;
-    fn div(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> where T: std::ops::DivAssign + Copy;
-    fn mtx_mul(mtx: &Matrix<T, L>, other: &Matrix<T, LO>) -> Result<Matrix<T, L>, String> where T: std::ops::Mul<Output = T> + std::ops::AddAssign + Default + Copy;
+    fn add(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> where T: std::ops::AddAssign + MatrixElement;
+    fn sub(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> where T: std::ops::SubAssign + MatrixElement;
+    fn hadamard_mul(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> where T: std::ops::MulAssign + MatrixElement;
+    fn div(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> where T: std::ops::DivAssign + MatrixElement;
+    fn mtx_mul(mtx: &Matrix<T, L>, other: &Matrix<T, LO>) -> Result<Matrix<T, L>, String> where T: std::ops::Mul<Output = T> + std::ops::AddAssign + MatrixElement + std::iter::Sum<T>;
 }
 
 pub struct NaiveAlgorithm;
 impl<T, L: MatrixLayout, LO: MatrixLayout> MatrixAlgorithm<T, L, LO> for NaiveAlgorithm {
     fn add(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String>
-    where T: std::ops::AddAssign + Copy
+    where T: std::ops::AddAssign + MatrixElement
     {
         if mtx.rows() != other.rows() || mtx.cols() != other.cols() {
             return Err("行列のサイズが一致しません".to_string());
         }
 
         for i in 0..mtx.rows() {
-            mtx.row_mut_iter(i).unwrap().zip(other.row_iter(i).unwrap()).for_each(|(a, b)| *a += *b);
+            mtx.row_iter_mut(i).unwrap().zip(other.row_iter(i).unwrap()).for_each(|(a, b)| *a += *b);
         }
         Ok(())
     }
 
     fn sub(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> 
-    where T: std::ops::SubAssign + Copy
+    where T: std::ops::SubAssign + MatrixElement
     {
         if mtx.rows() != other.rows() || mtx.cols() != other.cols() {
             return Err("行列のサイズが一致しません".to_string());
         }
 
         for i in 0..mtx.rows() {
-            mtx.row_mut_iter(i).unwrap().zip(other.row_iter(i).unwrap()).for_each(|(a, b)| *a -= *b);
+            mtx.row_iter_mut(i).unwrap().zip(other.row_iter(i).unwrap()).for_each(|(a, b)| *a -= *b);
         }
         Ok(())
     }
 
     fn hadamard_mul(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> 
-    where T: std::ops::MulAssign + Copy
+    where T: std::ops::MulAssign + MatrixElement
     {
         if mtx.rows() != other.rows() || mtx.cols() != other.cols() {
             return Err("行列のサイズが一致しません".to_string());
         }
 
         for i in 0..mtx.rows() {
-            mtx.row_mut_iter(i).unwrap().zip(other.row_iter(i).unwrap()).for_each(|(a, b)| *a *= *b);
+            mtx.row_iter_mut(i).unwrap().zip(other.row_iter(i).unwrap()).for_each(|(a, b)| *a *= *b);
         }
         Ok(())
     }
 
     fn div(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> 
-    where T: std::ops::DivAssign + Copy
+    where T: std::ops::DivAssign + MatrixElement
     {
         if mtx.rows() != other.rows() || mtx.cols() != other.cols() {
             return Err("行列のサイズが一致しません".to_string());
         }
 
         for i in 0..mtx.rows() {
-            mtx.row_mut_iter(i).unwrap().zip(other.row_iter(i).unwrap()).for_each(|(a, b)| *a /= *b);
+            mtx.row_iter_mut(i).unwrap().zip(other.row_iter(i).unwrap()).for_each(|(a, b)| *a /= *b);
         }
         Ok(())
     }
 
     fn mtx_mul(mtx: &Matrix<T, L>, other: &Matrix<T, LO>) -> Result<Matrix<T, L>, String> 
-    where T: Default + Copy + std::ops::Mul<Output = T> + std::ops::AddAssign
+    where T: std::ops::Mul<Output = T> + std::ops::AddAssign + MatrixElement + std::iter::Sum<T>
     {
         if mtx.cols() != other.rows() {
             return Err("行列のサイズが一致しません".to_string());
@@ -73,11 +75,87 @@ impl<T, L: MatrixLayout, LO: MatrixLayout> MatrixAlgorithm<T, L, LO> for NaiveAl
         let mut result = Matrix::with_size(mtx.rows(), other.cols());
         for i in 0..mtx.rows() {
             for j in 0..other.cols() {
-                let mut temp = T::default();
+                let temp = mtx.row_iter(i).unwrap()
+                    .zip(other.col_iter(j).unwrap())
+                    .map(|(&a, &b)| a * b)
+                    .sum::<T>();
+                *result.get_mut(i, j).unwrap() = temp;
+            }
+        }
+        Ok(result)
+    }
+}
+
+pub struct ParallelAlgorithm;
+impl<T, L: MatrixLayout, LO: MatrixLayout> MatrixAlgorithm<T, L, LO> for ParallelAlgorithm {
+    fn add(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String>
+    where
+        T: std::ops::AddAssign + Copy + Sync + Send
+    {
+        if mtx.rows() != other.rows() || mtx.cols() != other.cols() {
+            return Err("行列のサイズが一致しません".to_string());
+        }
+
+        for i in 0..mtx.rows() {
+            mtx.row_par_iter_mut(i).unwrap().zip(other.row_par_iter(i).unwrap()).for_each(|(a, b)| *a += *b);
+        }
+
+        Ok(())
+    }
+
+    fn sub(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> 
+    where T: std::ops::SubAssign + MatrixElement
+    {
+        if mtx.rows() != other.rows() || mtx.cols() != other.cols() {
+            return Err("行列のサイズが一致しません".to_string());
+        }
+
+        for i in 0..mtx.rows() {
+            mtx.row_par_iter_mut(i).unwrap().zip(other.row_par_iter(i).unwrap()).for_each(|(a, b)| *a -= *b);
+        }
+        Ok(())
+    }
+
+    fn hadamard_mul(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> 
+    where T: std::ops::MulAssign + MatrixElement
+    {
+        if mtx.rows() != other.rows() || mtx.cols() != other.cols() {
+            return Err("行列のサイズが一致しません".to_string());
+        }
+
+        for i in 0..mtx.rows() {
+            mtx.row_par_iter_mut(i).unwrap().zip(other.row_par_iter(i).unwrap()).for_each(|(a, b)| *a *= *b);
+        }
+        Ok(())
+    }
+
+    fn div(mtx: &mut Matrix<T, L>, other: &Matrix<T, LO>) -> Result<(), String> 
+    where T: std::ops::DivAssign + MatrixElement
+    {
+        if mtx.rows() != other.rows() || mtx.cols() != other.cols() {
+            return Err("行列のサイズが一致しません".to_string());
+        }
+
+        for i in 0..mtx.rows() {
+            mtx.row_par_iter_mut(i).unwrap().zip(other.row_par_iter(i).unwrap()).for_each(|(a, b)| *a /= *b);
+        }
+        Ok(())
+    }
 
-                mtx.row_iter(i).unwrap().zip(other.col_iter(j).unwrap()).for_each(|(a, b)| {
-                    temp += *a * *b;
-                });
+    fn mtx_mul(mtx: &Matrix<T, L>, other: &Matrix<T, LO>) -> Result<Matrix<T, L>, String> 
+    where T: std::ops::Mul<Output = T> + std::ops::AddAssign + MatrixElement + std::iter::Sum<T>
+    {
+        if mtx.cols() != other.rows() {
+            return Err("行列のサイズが一致しません".to_string());
+        }
+
+        let mut result = Matrix::with_size(mtx.rows(), other.cols());
+        for i in 0..mtx.rows() {
+            for j in 0..other.cols() {
+                let temp: T = mtx.row_par_iter(i).unwrap()
+                    .zip(other.col_par_iter(j).unwrap())
+                    .map(|(&a, &b)| a * b)
+                    .sum();
 
                 *result.get_mut(i, j).unwrap() = temp;
             }
diff --git a/rust/src/matrix_element.rs b/rust/src/matrix_element.rs
new file mode 100644
index 0000000..8023112
--- /dev/null
+++ b/rust/src/matrix_element.rs
@@ -0,0 +1 @@
+pub trait MatrixElement: Default + Copy + Sync + Send {}
\ No newline at end of file
diff --git a/rust/src/matrix_layout.rs b/rust/src/matrix_layout.rs
index 4e7b53c..90a264b 100644
--- a/rust/src/matrix_layout.rs
+++ b/rust/src/matrix_layout.rs
@@ -1,4 +1,4 @@
-pub trait MatrixLayout {
+pub trait MatrixLayout: Sync + Send {
     /// row, colから一次元配列のindexを取得する
     /// # 引数
     /// - `row`: 行番号