diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml
index 31000a2..b53e8ed 100644
--- a/.github/workflows/rust.yml
+++ b/.github/workflows/rust.yml
@@ -17,6 +17,6 @@ jobs:
     steps:
     - uses: actions/checkout@v3
     - name: Build
-      run: cargo build --verbose
+      run: cargo build --verbose --all-features
     - name: Run tests
-      run: cargo test --verbose
+      run: cargo test --verbose --all-features
diff --git a/Cargo.lock b/Cargo.lock
index 71411fa..e256324 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -51,9 +51,9 @@ checksum = "b35204fbdc0b3f4446b89fc1ac2cf84a8a68971995d0bf2e925ec7cd960f9cb3"
 
 [[package]]
 name = "cc"
-version = "1.2.37"
+version = "1.2.49"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "65193589c6404eb80b450d618eaf9a2cafaaafd57ecce47370519ef674a7bd44"
+checksum = "90583009037521a116abf44494efecd645ba48b6622457080f080b85544e2215"
 dependencies = [
  "find-msvc-tools",
  "shlex",
@@ -110,9 +110,9 @@ checksum = "877a4ace8713b0bcf2a4e7eec82529c029f1d0619886d18145fea96c3ffe5c0f"
 
 [[package]]
 name = "find-msvc-tools"
-version = "0.1.1"
+version = "0.1.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "7fd99930f64d146689264c637b5af2f0233a933bef0d8570e2526bf9e083192d"
+checksum = "3a3076410a55c90011c298b04d0cfa770b00fa04e1e3c97d3f6c9de105a03844"
 
 [[package]]
 name = "hashbrown"
@@ -144,9 +144,12 @@ dependencies = [
 
 [[package]]
 name = "indoc"
-version = "2.0.6"
+version = "2.0.7"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f4c7245a08504955605670dbf141fceab975f15ca21570696aebe9d2e71576bd"
+checksum = "79cf5c93f93228cf8efb3ba362535fb11199ac548a09ce117c9b1adc3030d706"
+dependencies = [
+ "rustversion",
+]
 
 [[package]]
 name = "js-sys"
@@ -160,9 +163,9 @@ dependencies = [
 
 [[package]]
 name = "libc"
-version = "0.2.175"
+version = "0.2.178"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "6a82ae493e598baaea5209805c49bbf2ea7de956d50d7da0da1164f9c6d28543"
+checksum = "37c93d8daa9d8a012fd8ab92f088405fb202ea0b6ab73ee2482ae66af4f42091"
 
 [[package]]
 name = "libm"
@@ -170,6 +173,16 @@ version = "0.2.15"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "f9fbbcab51052fe104eb5e5d351cf728d30a5be1fe14d9be8a3b097481fb97de"
 
+[[package]]
+name = "matrixmultiply"
+version = "0.3.10"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a06de3016e9fae57a36fd14dba131fccf49f74b40b7fbdb472f96e361ec71a08"
+dependencies = [
+ "autocfg",
+ "rawpointer",
+]
+
 [[package]]
 name = "memoffset"
 version = "0.9.1"
@@ -241,9 +254,9 @@ dependencies = [
 
 [[package]]
 name = "proc-macro2"
-version = "1.0.101"
+version = "1.0.103"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "89ae43fd86e4158d6db51ad8e2b80f313af9cc74f5c0e03ccb87de09998732de"
+checksum = "5ee95bc4ef87b8d5ba32e8b7714ccc834865276eab0aed5c9958d00ec45f49e8"
 dependencies = [
  "unicode-ident",
 ]
@@ -341,13 +354,19 @@ dependencies = [
 
 [[package]]
 name = "quote"
-version = "1.0.40"
+version = "1.0.42"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1885c039570dc00dcb4ff087a89e185fd56bae234ddc7f056a945bf36467248d"
+checksum = "a338cc41d27e6cc6dce6cefc13a0729dfbb81c262b1f519331575dd80ef3067f"
 dependencies = [
  "proc-macro2",
 ]
 
+[[package]]
+name = "rawpointer"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
+
 [[package]]
 name = "rancor"
 version = "0.1.1"
@@ -466,9 +485,9 @@ checksum = "e3a9fe34e3e7a50316060351f37187a3f546bce95496156754b601a5fa71b76e"
 
 [[package]]
 name = "syn"
-version = "2.0.106"
+version = "2.0.111"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "ede7c438028d4436d71104916910f5bb611972c5cfd7f89b8300a8186e6fada6"
+checksum = "390cc9a294ab71bdb1aa2e99d13be9c753cd2d7bd6560c77118597410c4d2e87"
 dependencies = [
  "proc-macro2",
  "quote",
@@ -498,9 +517,9 @@ checksum = "1f3ccbac311fea05f86f61904b462b55fb3df8837a366dfc601a0161d0532f20"
 
 [[package]]
 name = "unicode-ident"
-version = "1.0.19"
+version = "1.0.22"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f63a545481291138910575129486daeaf8ac54aee4387fe7906919f7830c7d9d"
+checksum = "9312f7c4f6ff9069b165498234ce8be658059c6728633667c526e27dc2cf1df5"
 
 [[package]]
 name = "unindent"
@@ -520,10 +539,11 @@ dependencies = [
 
 [[package]]
 name = "vec-utils"
-version = "0.2.5"
+version = "0.3.0"
 dependencies = [
  "assert_float_eq",
  "libm",
+ "matrixmultiply",
  "pretty_assertions",
  "rkyv",
  "serde",
@@ -531,7 +551,7 @@ dependencies = [
 
 [[package]]
 name = "vec-utils-py"
-version = "0.2.5"
+version = "0.3.0"
 dependencies = [
  "ordered-float",
  "pyo3",
diff --git a/README.md b/README.md
index 6eb8aa0..772759a 100644
--- a/README.md
+++ b/README.md
@@ -22,3 +22,5 @@ https://www.cs.utexas.edu/~ear/cs341/automatabook/AutomataTheoryBook.pdf
 https://dl.acm.org/doi/pdf/10.1145/360825.360855
 https://docs.google.com/document/d/1KkKC2-ozJkvbWQIAXeJ1MUGqxjn19c-Mmc7RtxFTA3c/edit?tab=t.0
 https://ntrs.nasa.gov/api/citations/19900013774/downloads/19900013774.pdf
+
+https://www.cs.unc.edu/techreports/96-043.pdf
diff --git a/update_version.py b/update_version.py
index f308ac5..d6f9996 100644
--- a/update_version.py
+++ b/update_version.py
@@ -1,4 +1,4 @@
-VERSION = "0.2.5"
+VERSION = "0.3.0"
 
 locations = [
     "./vec-utils/Cargo.toml",
diff --git a/vec-utils-py/Cargo.toml b/vec-utils-py/Cargo.toml
index 0dc423f..b17a808 100644
--- a/vec-utils-py/Cargo.toml
+++ b/vec-utils-py/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "vec-utils-py"
-version = "0.2.5"
+version = "0.3.0"
 edition = "2024"
 authors = ["Akari Harada <akaritwo0two@gmail.com>"]
 license = "GPL-3"
@@ -17,4 +17,4 @@ pyo3 = { version = "0.27.2", features = [
   "generate-import-lib"
 ] }
 rayon = "1.11.0"
-vec-utils = { path = "../vec-utils" }
+vec-utils = { path = "../vec-utils", features = ["matrix"] }
diff --git a/vec-utils-py/pyproject.toml b/vec-utils-py/pyproject.toml
index 7c3f99c..80d92ae 100644
--- a/vec-utils-py/pyproject.toml
+++ b/vec-utils-py/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "vec-utils-py"
-version = "0.2.5"
+version = "0.3.0"
 requires-python = ">=3.9"
 authors = [
   { name="Akari202", email="akaritwo0two@gmail.com" }
diff --git a/vec-utils-py/src/quat.rs b/vec-utils-py/src/quat.rs
index 282c4af..cb9feb9 100644
--- a/vec-utils-py/src/quat.rs
+++ b/vec-utils-py/src/quat.rs
@@ -104,7 +104,7 @@ impl Quat {
     }
 
     fn to_rotation_matrix(&self) -> [[f64; 3]; 3] {
-        self.inner.to_rotation_matrix()
+        self.inner.to_rotation_matrix().to_nested_arr()
     }
 
     fn rotate(&self, v: &Vec3d) -> Vec3d {
diff --git a/vec-utils/Cargo.toml b/vec-utils/Cargo.toml
index 2a532b0..7dd29c5 100644
--- a/vec-utils/Cargo.toml
+++ b/vec-utils/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "vec-utils"
-version = "0.2.5"
+version = "0.3.0"
 edition = "2024"
 authors = ["Akari Harada <akaritwo0two@gmail.com>"]
 license = "GPL-3"
@@ -10,6 +10,9 @@ name = "vec_utils"
 
 [dependencies]
 libm = { version = "0.2.15" }
+matrixmultiply = { version = "0.3.10", default-features = false, features = [
+  "cgemm"
+], optional = true }
 rkyv = { version = "0.8.12", default-features = false, features = ["bytecheck"], optional = true }
 serde = { version = "1.0.219", default-features = false, features = ["derive"], optional = true }
 
@@ -19,6 +22,7 @@ pretty_assertions = "1.4.1"
 
 [features]
 default = ["std"]
+matrix = ["dep:matrixmultiply"]
 rkyv = ["dep:rkyv"]
 serde = ["dep:serde"]
-std = ["rkyv?/std", "serde?/std"]
+std = ["matrixmultiply?/std", "rkyv?/std", "serde?/std"]
diff --git a/vec-utils/src/complex.rs b/vec-utils/src/complex.rs
index 2f43444..c5dff4f 100644
--- a/vec-utils/src/complex.rs
+++ b/vec-utils/src/complex.rs
@@ -1,3 +1,4 @@
+use core::cmp::Ordering;
 use core::fmt;
 use core::ops::{Add, Div, Index, Mul, Sub};
 
@@ -7,6 +8,7 @@ use crate::{
 };
 
 /// A complex number
+#[repr(C)]
 #[cfg_attr(feature = "serde", derive(serde::Deserialize, serde::Serialize))]
 #[cfg_attr(
     feature = "rkyv",
@@ -241,6 +243,42 @@ impl_single_op_variants_other!(
     "Divide a real number by a complex number"
 );
 
+impl PartialOrd for Complex {
+    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
+        let self_mag_sq = self.real * self.real + self.imaginary * self.imaginary;
+        let other_mag_sq = other.real * other.real + other.imaginary * other.imaginary;
+
+        self_mag_sq.partial_cmp(&other_mag_sq)
+    }
+}
+
+impl PartialEq<f64> for Complex {
+    fn eq(&self, other: &f64) -> bool {
+        self.imaginary == 0.0 && self.real == *other
+    }
+}
+
+impl PartialOrd<f64> for Complex {
+    fn partial_cmp(&self, other: &f64) -> Option<Ordering> {
+        let self_mag_sq = self.real * self.real + self.imaginary * self.imaginary;
+        let other_mag_sq = other * other;
+
+        self_mag_sq.partial_cmp(&other_mag_sq)
+    }
+}
+
+impl PartialEq<Complex> for f64 {
+    fn eq(&self, other: &Complex) -> bool {
+        other.eq(self)
+    }
+}
+
+impl PartialOrd<Complex> for f64 {
+    fn partial_cmp(&self, other: &Complex) -> Option<Ordering> {
+        other.partial_cmp(self).map(Ordering::reverse)
+    }
+}
+
 impl Index<usize> for Complex {
     type Output = f64;
 
diff --git a/vec-utils/src/geometry.rs b/vec-utils/src/geometry.rs
index 8842103..3c2a22f 100644
--- a/vec-utils/src/geometry.rs
+++ b/vec-utils/src/geometry.rs
@@ -1,6 +1,3 @@
-//! Geometry module
-//! This module contains geometric shapes and operations on them.
-
 /// Circles
 pub mod circle;
 /// Intersections
diff --git a/vec-utils/src/lib.rs b/vec-utils/src/lib.rs
index 334ad3b..5ff8679 100644
--- a/vec-utils/src/lib.rs
+++ b/vec-utils/src/lib.rs
@@ -1,12 +1,19 @@
 #![no_std]
-#![deny(missing_docs)]
 #![warn(clippy::pedantic)]
+#![deny(
+    missing_docs,
+    clippy::undocumented_unsafe_blocks,
+    clippy::unnecessary_safety_doc
+)]
 #![allow(
     clippy::must_use_candidate,
     clippy::many_single_char_names,
-    clippy::return_self_not_must_use
+    clippy::return_self_not_must_use,
+    clippy::derive_ord_xor_partial_ord,
+    incomplete_features
 )]
 #![cfg_attr(not(feature = "std"), feature(core_float_math))]
+#![cfg_attr(feature = "matrix", feature(generic_const_exprs))]
 #![doc(test(attr(deny(dead_code))))]
 //! A crate for 3D vector, quaternion, geometry, and matrix operations
 //! plus some miscellaneous other common things.
@@ -24,11 +31,9 @@ pub mod complex;
 pub mod geometry;
 /// Internal macros
 pub(crate) mod macros;
-/// Functions for working with matrices
-/// currently only 2x2, 3x3, and 4x4 matrices are supported
-/// with functions for calculating the determinant, minor, and cofactor
-/// only available on std until i get around to fixing it or <https://github.com/rust-lang/rust/issues/137578> gets merged
-#[cfg(feature = "std")]
+/// Matrix operations and functions.
+/// Requires the "matrix" feature
+#[cfg(feature = "matrix")]
 pub mod matrix;
 /// Quaternion operations and functions
 pub mod quat;
diff --git a/vec-utils/src/matrix.rs b/vec-utils/src/matrix.rs
index a0d652b..0bb1bc0 100644
--- a/vec-utils/src/matrix.rs
+++ b/vec-utils/src/matrix.rs
@@ -1,353 +1,9 @@
-/// Functions for working with 2x2 matrices
-pub mod matrix2x2 {
-    use crate::complex::Complex;
-
-    /// Calculate the determinant of a 2x2 matrix
-    pub fn determinant(matrix: &[[f64; 2]; 2]) -> f64 {
-        matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]
-    }
-
-    /// Calculate the eigenvalues of a 2x2 matrix
-    /// returns a tuple of the eigenvalues as complex numbers
-    pub fn eigenvalues(matrix: &[[f64; 2]; 2]) -> (Complex, Complex) {
-        let mean = f64::midpoint(matrix[0][0], matrix[1][1]);
-        let determinant = determinant(matrix);
-        let discriminant = mean.powi(2) - determinant;
-        let eigenvalue1 = Complex::new(mean, 0.0) + Complex::sqrt(discriminant);
-        let eigenvalue2 = Complex::new(mean, 0.0) - Complex::sqrt(discriminant);
-        (eigenvalue1, eigenvalue2)
-    }
-
-    /// Calculate the eigenvectors of a 2x2 matrix
-    /// returns a tuple of the eigenvectors as 2D arrays
-    pub fn eigenvectors(matrix: &[[f64; 2]; 2]) -> ([f64; 2], [f64; 2]) {
-        let (eigenvalue1, eigenvalue2) = eigenvalues(matrix);
-        let mut eigenvector1 = [0.0; 2];
-        let mut eigenvector2 = [0.0; 2];
-        if eigenvalue1.imaginary == 0.0 {
-            if matrix[0][0] - eigenvalue1.real != 0.0 {
-                eigenvector1[0] = matrix[0][1] / (matrix[0][0] - eigenvalue1.real);
-                eigenvector1[1] = 1.0;
-            } else if matrix[1][0] != 0.0 {
-                eigenvector1[0] = 1.0;
-                eigenvector1[1] = matrix[1][1] / (matrix[1][0] - eigenvalue1.real);
-            }
-        }
-        if eigenvalue2.imaginary == 0.0 {
-            if matrix[0][0] - eigenvalue2.real != 0.0 {
-                eigenvector2[0] = matrix[0][1] / (matrix[0][0] - eigenvalue2.real);
-                eigenvector2[1] = 1.0;
-            } else if matrix[1][0] != 0.0 {
-                eigenvector2[0] = 1.0;
-                eigenvector2[1] = matrix[1][1] / (matrix[1][0] - eigenvalue2.real);
-            }
-        }
-        (eigenvector1, eigenvector2)
-    }
-}
-
-/// Functions for working with 3x3 matrices
-pub mod matrix3x3 {
-    use crate::vec3d::Vec3d;
-
-    /// Calculate the determinant of a 3x3 matrix
-    pub fn determinant(matrix: &[[f64; 3]; 3]) -> f64 {
-        matrix[0][0] * matrix[1][1] * matrix[2][2]
-            + matrix[0][1] * matrix[1][2] * matrix[2][0]
-            + matrix[0][2] * matrix[1][0] * matrix[2][1]
-            - matrix[0][2] * matrix[1][1] * matrix[2][0]
-            - matrix[0][1] * matrix[1][0] * matrix[2][2]
-            - matrix[0][0] * matrix[1][2] * matrix[2][1]
-    }
-
-    /// Calculate the minor of a 3x3 matrix given a row and column index
-    pub fn minor(matrix: &[[f64; 3]; 3], row: usize, col: usize) -> f64 {
-        let mut minor = [[0.0; 2]; 2];
-        for i in 0..3 {
-            for j in 0..3 {
-                if i != row && j != col {
-                    let mut m = i;
-                    let mut n = j;
-                    if i > row {
-                        m -= 1;
-                    }
-                    if j > col {
-                        n -= 1;
-                    }
-                    minor[m][n] = matrix[i][j];
-                }
-            }
-        }
-        super::matrix2x2::determinant(&minor)
-    }
-
-    /// Calculate the cofactor of a 3x3 matrix given a row and column index
-    pub fn cofactor(matrix: &[[f64; 3]; 3], row: usize, col: usize) -> f64 {
-        let minor = minor(matrix, row, col);
-        let factor = [[1.0, -1.0, 1.0], [-1.0, 1.0, -1.0], [1.0, -1.0, 1.0]];
-        factor[row][col] * minor
-    }
-
-    /// Get the cofactor matrix of a 3x3 matrix
-    pub fn cofactor_matrix(matrix: &[[f64; 3]; 3]) -> [[f64; 3]; 3] {
-        let mut cofactor_matrix = [[0.0; 3]; 3];
-        for i in 0..3 {
-            for j in 0..3 {
-                cofactor_matrix[i][j] = cofactor(matrix, i, j);
-            }
-        }
-        cofactor_matrix
-    }
-
-    /// Transpose a 3x3 matrix
-    /// i.e. swap the rows and columns
-    pub fn transpose(matrix: &[[f64; 3]; 3]) -> [[f64; 3]; 3] {
-        let mut transpose = [[0.0; 3]; 3];
-        for i in 0..3 {
-            for j in 0..3 {
-                transpose[i][j] = matrix[j][i];
-            }
-        }
-        transpose
-    }
-
-    /// Calculate the adjoint of a 3x3 matrix
-    /// i.e. the transpose of the cofactor matrix
-    pub fn adjoint(matrix: &[[f64; 3]; 3]) -> [[f64; 3]; 3] {
-        let cofactor_matrix = cofactor_matrix(matrix);
-        let transpose_matrix = transpose(&cofactor_matrix);
-        dbg!(cofactor_matrix);
-        dbg!(transpose_matrix);
-        transpose_matrix
-    }
-
-    /// Vector multiplication of a matrix with a Vec3d
-    pub fn mul(matrix: &[[f64; 3]; 3], vector: &Vec3d) -> Vec3d {
-        let mut result: [f64; 3] = [0.0; 3];
-        for (i, j) in matrix.iter().enumerate() {
-            // result[i] = vector.dot(Vec3d::from_slice(j));
-            result[i] = Vec3d::from_slice(j).dot(vector);
-        }
-        Vec3d::from_slice(&result)
-    }
-
-    // Calculate the eigenvalues of a 3x3 matrix
-    // returns a tuple of the eigenvalues as complex numbers
-    // pub fn eigenvalues(matrix: &[[f64; 3]; 3]) -> (Complex, Complex, Complex) {
-    //
-    // }
-}
-
-/// Functions for working with 4x4 matrices
-pub mod matrix4x4 {
-    /// Calculate the determinant of a 4x4 matrix
-    pub fn determinant(matrix: &[[f64; 4]; 4]) -> f64 {
-        matrix[0][0] * matrix[1][1] * matrix[2][2] * matrix[3][3]
-            + matrix[0][0] * matrix[1][2] * matrix[2][3] * matrix[3][1]
-            + matrix[0][0] * matrix[1][3] * matrix[2][1] * matrix[3][2]
-            + matrix[0][1] * matrix[1][0] * matrix[2][3] * matrix[3][2]
-            + matrix[0][1] * matrix[1][2] * matrix[2][0] * matrix[3][3]
-            + matrix[0][1] * matrix[1][3] * matrix[2][2] * matrix[3][0]
-            + matrix[0][2] * matrix[1][0] * matrix[2][1] * matrix[3][3]
-            + matrix[0][2] * matrix[1][1] * matrix[2][3] * matrix[3][0]
-            + matrix[0][2] * matrix[1][3] * matrix[2][0] * matrix[3][1]
-            + matrix[0][3] * matrix[1][0] * matrix[2][2] * matrix[3][1]
-            + matrix[0][3] * matrix[1][1] * matrix[2][0] * matrix[3][2]
-            + matrix[0][3] * matrix[1][2] * matrix[2][1] * matrix[3][0]
-            - matrix[0][0] * matrix[1][1] * matrix[2][3] * matrix[3][2]
-            - matrix[0][0] * matrix[1][2] * matrix[2][1] * matrix[3][3]
-            - matrix[0][0] * matrix[1][3] * matrix[2][2] * matrix[3][1]
-            - matrix[0][1] * matrix[1][0] * matrix[2][2] * matrix[3][3]
-            - matrix[0][1] * matrix[1][2] * matrix[2][3] * matrix[3][0]
-            - matrix[0][1] * matrix[1][3] * matrix[2][0] * matrix[3][2]
-            - matrix[0][2] * matrix[1][0] * matrix[2][3] * matrix[3][1]
-            - matrix[0][2] * matrix[1][1] * matrix[2][0] * matrix[3][3]
-            - matrix[0][2] * matrix[1][3] * matrix[2][1] * matrix[3][0]
-            - matrix[0][3] * matrix[1][0] * matrix[2][1] * matrix[3][2]
-            - matrix[0][3] * matrix[1][1] * matrix[2][2] * matrix[3][0]
-            - matrix[0][3] * matrix[1][2] * matrix[2][0] * matrix[3][1]
-    }
-
-    /// Calculate the minor of a 4x4 matrix given a row and column index
-    pub fn minor(matrix: &[[f64; 4]; 4], row: usize, col: usize) -> f64 {
-        let mut minor = [[0.0; 3]; 3];
-        for i in 0..4 {
-            for j in 0..4 {
-                if i != row && j != col {
-                    let mut m = i;
-                    let mut n = j;
-                    if i > row {
-                        m -= 1;
-                    }
-                    if j > col {
-                        n -= 1;
-                    }
-                    minor[m][n] = matrix[i][j];
-                }
-            }
-        }
-        super::matrix3x3::determinant(&minor)
-    }
-
-    /// Calculate the cofactor of a 4x4 matrix given a row and column index
-    pub fn cofactor(matrix: &[[f64; 4]; 4], row: usize, col: usize) -> f64 {
-        let minor = minor(matrix, row, col);
-        let factor = [
-            [1.0, -1.0, 1.0, -1.0],
-            [-1.0, 1.0, -1.0, 1.0],
-            [1.0, -1.0, 1.0, -1.0],
-            [-1.0, 1.0, -1.0, 1.0]
-        ];
-        factor[row][col] * minor
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    mod tests2x2 {
-        use super::super::matrix2x2;
-
-        #[test]
-        fn test_matrix2x2_determinant() {
-            let matrix = [[1.0, 2.0], [3.0, 4.0]];
-            assert_eq!(matrix2x2::determinant(&matrix), -2.0);
-        }
-
-        #[test]
-        fn test_matrix2x2_eigenvalues() {
-            let matrix = [[8.0, 4.0], [4.0, 8.0]];
-            let (eigenvalue1, eigenvalue2) = matrix2x2::eigenvalues(&matrix);
-            assert_eq!(eigenvalue1.real, 12.0);
-            assert_eq!(eigenvalue1.imaginary, 0.0);
-            assert_eq!(eigenvalue2.real, 4.0);
-            assert_eq!(eigenvalue2.imaginary, 0.0);
-        }
-
-        #[test]
-        fn test_matrix2x2_eigenvectors() {
-            let matrix = [[8.0, 4.0], [4.0, 8.0]];
-            let (eigenvector1, eigenvector2) = matrix2x2::eigenvectors(&matrix);
-            assert_eq!(eigenvector1[0], -1.0);
-            assert_eq!(eigenvector1[1], 1.0);
-            assert_eq!(eigenvector2[0], 1.0);
-            assert_eq!(eigenvector2[1], 1.0);
-        }
-    }
-
-    mod tests3x3 {
-        use super::super::matrix3x3;
-
-        #[test]
-        fn test_matrix3x3_determinant() {
-            let matrix = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]];
-            assert_eq!(matrix3x3::determinant(&matrix), 0.0);
-        }
-
-        #[test]
-        fn test_matrix3x3_minor() {
-            let matrix = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]];
-            assert_eq!(matrix3x3::minor(&matrix, 0, 0), -3.0);
-        }
-
-        #[test]
-        fn test_matrix3x3_cofactor() {
-            let matrix = [[1.0, 4.0, 7.0], [3.0, 0.0, 5.0], [-1.0, 9.0, 11.0]];
-            assert_eq!(matrix3x3::cofactor(&matrix, 1, 2), -13.0);
-            assert_eq!(matrix3x3::cofactor(&matrix, 2, 0), 20.0);
-        }
-
-        #[test]
-        fn test_matrix3x3_cofactor_matrix() {
-            let matrix = [[1.0, 4.0, 7.0], [3.0, 0.0, 5.0], [-1.0, 9.0, 11.0]];
-            let cofactor = [
-                [-45.0, -38.0, 27.0],
-                [19.0, 18.0, -13.0],
-                [20.0, 16.0, -12.0]
-            ];
-            let cofactor_matrix = matrix3x3::cofactor_matrix(&matrix);
-            assert_eq!(cofactor_matrix[0][0], cofactor[0][0]);
-            assert_eq!(cofactor_matrix[0][1], cofactor[0][1]);
-            assert_eq!(cofactor_matrix[0][2], cofactor[0][2]);
-            assert_eq!(cofactor_matrix[1][0], cofactor[1][0]);
-            assert_eq!(cofactor_matrix[1][1], cofactor[1][1]);
-            assert_eq!(cofactor_matrix[1][2], cofactor[1][2]);
-            assert_eq!(cofactor_matrix[2][0], cofactor[2][0]);
-            assert_eq!(cofactor_matrix[2][1], cofactor[2][1]);
-            assert_eq!(cofactor_matrix[2][2], cofactor[2][2]);
-        }
-
-        #[test]
-        fn test_matrix3x3_transpose() {
-            let matrix = [[1.0, 4.0, 7.0], [3.0, 0.0, 5.0], [-1.0, 9.0, 11.0]];
-            let transpose = [[1.0, 3.0, -1.0], [4.0, 0.0, 9.0], [7.0, 5.0, 11.0]];
-            let transpose_matrix = matrix3x3::transpose(&matrix);
-            assert_eq!(transpose_matrix[0][0], transpose[0][0]);
-            assert_eq!(transpose_matrix[0][1], transpose[0][1]);
-            assert_eq!(transpose_matrix[0][2], transpose[0][2]);
-            assert_eq!(transpose_matrix[1][0], transpose[1][0]);
-            assert_eq!(transpose_matrix[1][1], transpose[1][1]);
-            assert_eq!(transpose_matrix[1][2], transpose[1][2]);
-            assert_eq!(transpose_matrix[2][0], transpose[2][0]);
-            assert_eq!(transpose_matrix[2][1], transpose[2][1]);
-            assert_eq!(transpose_matrix[2][2], transpose[2][2]);
-        }
-
-        #[test]
-        fn test_matrix3x3_adjoint() {
-            let matrix = [[1.0, 4.0, 7.0], [3.0, 0.0, 5.0], [-1.0, 9.0, 11.0]];
-            let adjoint = [
-                [-45.0, 19.0, 20.0],
-                [-38.0, 18.0, 16.0],
-                [27.0, -13.0, -12.0]
-            ];
-            let adjoint_matrix = matrix3x3::adjoint(&matrix);
-            dbg!(adjoint_matrix);
-            assert_eq!(adjoint_matrix[0][0], adjoint[0][0]);
-            assert_eq!(adjoint_matrix[0][1], adjoint[0][1]);
-            assert_eq!(adjoint_matrix[0][2], adjoint[0][2]);
-            assert_eq!(adjoint_matrix[1][0], adjoint[1][0]);
-            assert_eq!(adjoint_matrix[1][1], adjoint[1][1]);
-            assert_eq!(adjoint_matrix[1][2], adjoint[1][2]);
-            assert_eq!(adjoint_matrix[2][0], adjoint[2][0]);
-            assert_eq!(adjoint_matrix[2][1], adjoint[2][1]);
-            assert_eq!(adjoint_matrix[2][2], adjoint[2][2]);
-        }
-    }
-
-    mod tests4x4 {
-        use super::super::matrix4x4;
-
-        #[test]
-        fn test_matrix4x4_determinant() {
-            let matrix = [
-                [1.0, 2.0, 3.0, 4.0],
-                [5.0, 6.0, 7.0, 8.0],
-                [9.0, 10.0, 11.0, 12.0],
-                [13.0, 14.0, 15.0, 16.0]
-            ];
-            assert_eq!(matrix4x4::determinant(&matrix), 0.0);
-        }
-
-        #[test]
-        fn test_matrix4x4_minor() {
-            let matrix = [
-                [1.0, 2.0, 3.0, 4.0],
-                [5.0, 6.0, 7.0, 8.0],
-                [9.0, 10.0, 11.0, 12.0],
-                [13.0, 14.0, 15.0, 16.0]
-            ];
-            assert_eq!(matrix4x4::minor(&matrix, 0, 0), 0.0);
-        }
-
-        #[test]
-        fn test_matrix4x4_cofactor() {
-            let matrix = [
-                [1.0, 2.0, 3.0, 4.0],
-                [5.0, 6.0, 7.0, 8.0],
-                [9.0, 10.0, 11.0, 12.0],
-                [13.0, 14.0, 15.0, 16.0]
-            ];
-            assert_eq!(matrix4x4::cofactor(&matrix, 0, 0), 0.0);
-        }
-    }
-}
+/// Complex valued matracies
+pub mod complex;
+/// Generic matrix implementation.
+/// Not for general use, use the aliases [`crate::matrix::complex::CMatrix`] or [`crate::matrix::real::Matrix`]
+pub mod generic;
+/// Real valued matracies
+pub mod real;
+/// Traits for matrix types
+mod traits;
diff --git a/vec-utils/src/matrix/complex.rs b/vec-utils/src/matrix/complex.rs
new file mode 100644
index 0000000..912de95
--- /dev/null
+++ b/vec-utils/src/matrix/complex.rs
@@ -0,0 +1,106 @@
+use core::ops::Mul;
+
+use matrixmultiply::zgemm;
+
+use crate::complex::Complex;
+#[doc(inline)]
+use crate::matrix::generic::GMatrix;
+
+/// A 2d complex matrix of width R and height C
+pub type CMatrix<const R: usize, const C: usize> = GMatrix<R, C, Complex>;
+
+/// An alias for 2x2 matracies
+pub type CMatrix2x2 = CMatrix<2, 2>;
+/// An alias for 3x3 matracies
+pub type CMatrix3x3 = CMatrix<3, 3>;
+
+impl<const R: usize, const C: usize> CMatrix<R, C> where [Complex; R * C]: Sized {}
+
+impl<const R: usize, const C: usize, const U: usize> Mul<CMatrix<U, C>> for CMatrix<R, U>
+where
+    [Complex; R * C]: Sized,
+    [Complex; R * U]: Sized,
+    [Complex; U * C]: Sized
+{
+    type Output = CMatrix<R, C>;
+
+    fn mul(self, rhs: CMatrix<U, C>) -> Self::Output {
+        let mut result = CMatrix::<R, C>::zeros();
+        // Safety: [Complex; R * C] and [[f64; 2]; R * C] have the exact same memory layout. zgemm
+        // is an unsafe function.
+        unsafe {
+            let lhs_ptr = (&raw const self.values).cast::<[f64; 2]>();
+            let rhs_ptr = (&raw const rhs.values).cast::<[f64; 2]>();
+            let result_ptr = (&raw mut result.values).cast::<[f64; 2]>();
+            zgemm(
+                matrixmultiply::CGemmOption::Standard,
+                matrixmultiply::CGemmOption::Standard,
+                R,
+                U,
+                C,
+                [1.0, 0.0],
+                lhs_ptr,
+                U.cast_signed(),
+                1,
+                rhs_ptr,
+                C.cast_signed(),
+                1,
+                [0.0; 2],
+                result_ptr,
+                C.cast_signed(),
+                1
+            );
+        }
+        result
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use pretty_assertions::assert_eq;
+
+    use super::*;
+
+    #[test]
+    fn test_mul() {
+        let a = CMatrix2x2::from_nested_arr([
+            [Complex::new(1.0, 0.0), Complex::new(0.0, 1.0)],
+            [Complex::new(0.0, 0.0), Complex::new(2.0, 0.0)]
+        ]);
+
+        let b = CMatrix2x2::from_nested_arr([
+            [Complex::new(0.0, 1.0), Complex::new(1.0, 0.0)],
+            [Complex::new(1.0, 0.0), Complex::new(0.0, 1.0)]
+        ]);
+
+        let correct = CMatrix2x2::from_nested_arr([
+            [Complex::new(0.0, 2.0), Complex::new(0.0, 0.0)],
+            [Complex::new(2.0, 0.0), Complex::new(0.0, 2.0)]
+        ]);
+
+        let result = a * b;
+
+        assert_eq!(result, correct);
+        let a = CMatrix3x3::from_nested_arr([
+            [
+                Complex::new(1.0, 2.0),
+                Complex::new(3.0, 4.0),
+                Complex::new(5.0, 6.0)
+            ],
+            [
+                Complex::new(7.0, 8.0),
+                Complex::new(9.0, 0.0),
+                Complex::new(1.0, 1.0)
+            ],
+            [
+                Complex::new(0.0, 0.0),
+                Complex::new(2.0, 2.0),
+                Complex::new(3.0, 3.0)
+            ]
+        ]);
+
+        let zero_mat = CMatrix3x3::zeros();
+        let result = a * zero_mat;
+        assert_eq!(result, zero_mat);
+    }
+}
diff --git a/vec-utils/src/matrix/generic.rs b/vec-utils/src/matrix/generic.rs
new file mode 100644
index 0000000..9090268
--- /dev/null
+++ b/vec-utils/src/matrix/generic.rs
@@ -0,0 +1,663 @@
+use core::fmt::Debug;
+use core::ops::{Add, Div, Index, IndexMut, Mul, Sub};
+#[cfg(feature = "std")]
+use std::vec::Vec;
+
+use crate::matrix::traits::{Oneable, Signed, Zeroable};
+
+/// A generic 2d matrix of width R and height C
+// TODO: I would like to add a generic is row major switch
+#[derive(Debug, Copy, Clone, PartialEq)]
+pub struct GMatrix<const R: usize, const C: usize, T>
+where
+    [T; R * C]: Sized
+{
+    pub(crate) values: [T; R * C]
+}
+
+/// An alias for 2x2 matracies
+pub type GMatrix2x2<T> = GMatrix<2, 2, T>;
+/// An alias for 3x3 matracies
+pub type GMatrix3x3<T> = GMatrix<3, 3, T>;
+
+impl<const R: usize, const C: usize, T> GMatrix<R, C, T>
+where
+    [T; R * C]: Sized,
+    T: Debug
+        + Oneable
+        + Zeroable
+        + Copy
+        + Clone
+        + PartialEq
+        + Signed
+        + PartialOrd
+        + Mul<Output = T>
+        + Div<Output = T>
+        + Sub<Output = T>
+        + Add<Output = T>
+{
+    const IS_2X2: bool = R == 2 && C == 2;
+    const IS_3X3: bool = R == 3 && C == 3;
+    const IS_ONE_DIMM: bool = R == 1 || C == 1;
+    const IS_SQUARE: bool = R == C;
+
+    /// Create a matrix from nested vectors.
+    /// # Panics
+    /// A misshapen vector, i.e. one that's not of length C or that doesnt contain vectors of exclusively length R
+    #[cfg(feature = "std")]
+    pub fn from_nested_vec(values: Vec<Vec<T>>) -> Self {
+        let flattened: Vec<T> = values.into_iter().flatten().collect();
+        let values: [T; R * C] = flattened
+            .try_into()
+            .expect("Input dimensions do not match Matrix size R * C");
+        Self { values }
+    }
+
+    /// Create a matrix from nested arrays.
+    pub fn from_nested_arr(values: [[T; C]; R]) -> Self {
+        // Safety: [[T; C]; R] and [T; R * C] have the exact same memory layout
+        let flat_values = unsafe {
+            let ptr = (&raw const values).cast::<[T; R * C]>();
+            core::ptr::read(ptr)
+        };
+
+        let _ = values;
+
+        Self {
+            values: flat_values
+        }
+    }
+
+    /// Convert a matrix into nested arrays.
+    pub fn to_nested_arr(&self) -> [[T; C]; R] {
+        // Safety: [[T; C]; R] and [T; R * C] have the exact same memory layout
+        let nested_values = unsafe {
+            let ptr = (&raw const self.values).cast::<[[T; C]; R]>();
+            core::ptr::read(ptr)
+        };
+
+        let _ = self.values;
+
+        nested_values
+    }
+
+    /// Create a matrix filled with zeros
+    pub fn zeros() -> Self {
+        Self {
+            values: [T::zero(); R * C]
+        }
+    }
+
+    /// Create a matrix filled with ones
+    pub fn ones() -> Self {
+        Self {
+            values: [T::one(); R * C]
+        }
+    }
+
+    /// Determines if the matrix is square
+    pub fn is_square() -> bool {
+        Self::IS_SQUARE
+    }
+
+    /// Counts the number of nonzero values
+    pub fn count_nonzero(&self) -> usize {
+        self.values
+            .iter()
+            .fold(0, |acc, i| if i.is_zero() { acc } else { acc + 1 })
+    }
+
+    /// Returns the diagonal elements
+    #[cfg(feature = "std")]
+    pub fn diagonals(&self) -> Vec<T> {
+        let min_dimm = R.min(C);
+        (0..min_dimm).map(|i| self.values[i + i * C]).collect()
+    }
+
+    /// Checks if the matrix is upper triangluar
+    /// This does not check if its strictly upper triangluar
+    pub fn is_upper_triangular(&self) -> bool {
+        for row in 1..R {
+            for col in 0..row.min(C) {
+                if !self.values[row * C + col].is_zero() {
+                    return false;
+                }
+            }
+        }
+        true
+    }
+
+    /// Checks if the matrix is lower triangluar
+    /// This does not check if its strictly lower triangluar
+    pub fn is_lower_triangular(&self) -> bool {
+        for row in 0..R.min(C) {
+            for col in (row + 1)..C {
+                if !self.values[row * C + col].is_zero() {
+                    return false;
+                }
+            }
+        }
+        true
+    }
+
+    /// Checks if the matrix is a diagonal matrix
+    pub fn is_diagonal(&self) -> bool {
+        for row in 0..R {
+            for col in 0..row.min(C) {
+                if !self.values[row * C + col].is_zero() {
+                    return false;
+                }
+            }
+            if row + 1 < C {
+                for col in (row + 1)..C {
+                    if !self.values[row * C + col].is_zero() {
+                        return false;
+                    }
+                }
+            }
+        }
+        true
+    }
+
+    /// Iterates over the matrix with enumerated position values
+    pub fn iter_indexed(&self) -> impl Iterator<Item = ((usize, usize), &T)> {
+        self.values.iter().enumerate().map(|(idx, val)| {
+            let r = idx / C;
+            let c = idx % C;
+            ((r, c), val)
+        })
+    }
+
+    /// Iterates over the matrix
+    pub fn iter(&self) -> impl Iterator<Item = &T> {
+        self.values.iter()
+    }
+
+    /// Mutably iterates over the matrix
+    pub fn iter_mut(&mut self) -> impl Iterator<Item = &mut T> {
+        self.values.iter_mut()
+    }
+
+    /// Mutably iterates over the matrix with enumerated position values
+    pub fn iter_indexed_mut(&mut self) -> impl Iterator<Item = ((usize, usize), &mut T)> {
+        self.values.iter_mut().enumerate().map(|(idx, val)| {
+            let r = idx / C;
+            let c = idx % C;
+            ((r, c), val)
+        })
+    }
+
+    /// Calculates the determinant of the matrix
+    /// # Panics
+    /// If the matrix is not square
+    pub fn determinant(&self) -> T {
+        if Self::IS_2X2 {
+            self.values[0] * self.values[3] - self.values[1] * self.values[2]
+        } else if Self::IS_3X3 {
+            self.values[0] * self.values[4] * self.values[8]
+                + self.values[1] * self.values[5] * self.values[6]
+                + self.values[2] * self.values[3] * self.values[7]
+                - self.values[2] * self.values[4] * self.values[6]
+                - self.values[1] * self.values[3] * self.values[8]
+                - self.values[0] * self.values[5] * self.values[7]
+        } else if Self::IS_SQUARE {
+            let mut lu = self.values;
+            let mut sign = T::one();
+
+            for i in 0..R {
+                let mut pivot = i;
+                for j in (i + 1)..R {
+                    if lu[j * C + i].abs() > lu[pivot * C + i].abs() {
+                        pivot = j;
+                    }
+                }
+
+                if lu[pivot * C + i].is_zero() {
+                    return T::zero();
+                }
+
+                if pivot != i {
+                    for k in 0..C {
+                        lu.swap(i * C + k, pivot * C + k);
+                    }
+                    sign.flip();
+                }
+
+                for j in (i + 1)..R {
+                    let factor = lu[j * C + i] / lu[i * C + i];
+                    for k in (i + 1)..C {
+                        let val = lu[i * C + k];
+                        lu[j * C + k] = lu[j * C + k] - factor * val;
+                    }
+                }
+            }
+
+            let mut det = sign;
+            for i in 0..R {
+                det = det * lu[i * C + i];
+            }
+
+            det
+        } else {
+            panic!("Cannot take the determinant of a non square matrix");
+        }
+    }
+
+    /// Returns the minor submatrix by removing the provided row and column
+    /// # Panics
+    /// Will panic if the generic ROUT and COUT are not 1 less than R and C
+    pub fn get_minor_submatrix<const ROUT: usize, const COUT: usize>(
+        &self,
+        row_to_skip: usize,
+        col_to_skip: usize
+    ) -> GMatrix<ROUT, COUT, T>
+    where
+        [T; ROUT * COUT]: Sized
+    {
+        assert_eq!(ROUT, R - 1, "Wrong sized generics. Output rows must be R-1");
+        assert_eq!(
+            COUT,
+            C - 1,
+            "Wrong sized generics. Output columns must be C-1"
+        );
+
+        let mut sub_values = [T::zero(); ROUT * COUT];
+        let mut sub_idx = 0;
+
+        for r in 0..R {
+            if r == row_to_skip {
+                continue;
+            }
+            for c in 0..C {
+                if c == col_to_skip {
+                    continue;
+                }
+                sub_values[sub_idx] = self.values[r * C + c];
+                sub_idx += 1;
+            }
+        }
+
+        GMatrix { values: sub_values }
+    }
+
+    /// Calculates the cofactor of the matrix at the provided row and column
+    pub fn cofactor(&self, row: usize, col: usize) -> T
+    where
+        [T; (R - 1) * (C - 1)]: Sized
+    {
+        let sub = self.get_minor_submatrix::<{ R - 1 }, { C - 1 }>(row, col);
+        let minor_det = sub.determinant();
+
+        if (row + col).is_multiple_of(2) {
+            minor_det
+        } else {
+            T::zero() - minor_det
+        }
+    }
+
+    /// Returns the full matrix of cofactors
+    pub fn cofactor_matrix(&self) -> GMatrix<R, C, T>
+    where
+        [T; (R - 1) * (C - 1)]: Sized
+    {
+        let mut cofactors = [T::zero(); R * C];
+        for r in 0..R {
+            for c in 0..C {
+                cofactors[r * C + c] = self.cofactor(r, c);
+            }
+        }
+        GMatrix { values: cofactors }
+    }
+
+    /// Transposes the matrix.
+    /// For matriacies with a dimension of 1 this is zerocopy.
+    /// Otherwise it has to touch every value.
+    pub fn transpose(&self) -> GMatrix<C, R, T>
+    where
+        [T; C * R]: Sized
+    {
+        if Self::IS_ONE_DIMM {
+            // Safety: Both GMatrix<R, C, T> and GMatrix<C, R, T> have the same size and
+            // both represent a contiguous strip of memory.
+            // transmute_copy is used to move the bits into the new type.
+            unsafe {
+                let result = core::mem::transmute_copy::<Self, GMatrix<C, R, T>>(self);
+                let _ = self;
+                result
+            }
+        } else {
+            // TODO: implement blocking for bigger matracies
+            let mut output = [T::zero(); C * R];
+            for row in 0..R {
+                for col in 0..C {
+                    output[col * R + row] = self.values[row * C + col];
+                }
+            }
+            GMatrix::<C, R, T> { values: output }
+        }
+    }
+
+    /// Calculate the adjoint of a 3x3 matrix
+    /// i.e. the transpose of the cofactor matrix
+    pub fn adjoint(&self) -> GMatrix<C, R, T>
+    where
+        [T; (R - 1) * (C - 1)]: Sized,
+        [T; C * R]: Sized
+    {
+        self.cofactor_matrix().transpose()
+    }
+}
+
+impl<const R: usize, const C: usize, T> Index<[usize; 2]> for GMatrix<R, C, T>
+where
+    [T; R * C]: Sized
+{
+    type Output = T;
+
+    fn index(&self, idx: [usize; 2]) -> &Self::Output {
+        let [row, col] = idx;
+        assert!(
+            row < R && col < C,
+            "Index [{row}, {col}] is out of bounds for matrix of shape [{R}, {C}]"
+        );
+        &self.values[row * C + col]
+    }
+}
+
+impl<const R: usize, const C: usize, T> IndexMut<[usize; 2]> for GMatrix<R, C, T>
+where
+    [T; R * C]: Sized
+{
+    fn index_mut(&mut self, idx: [usize; 2]) -> &mut Self::Output {
+        let [row, col] = idx;
+        assert!(
+            row < R && col < C,
+            "Index [{row}, {col}] is out of bounds for matrix of shape [{R}, {C}]"
+        );
+        &mut self.values[row * C + col]
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use assert_float_eq::assert_f64_near;
+    use pretty_assertions::assert_eq;
+
+    use super::*;
+    use crate::complex::Complex;
+
+    #[test]
+    fn test_constructors_f64() {
+        let zero = GMatrix2x2::<f64>::zeros();
+        let one = GMatrix2x2::<f64>::ones();
+
+        for i in zero.iter() {
+            assert_f64_near!(*i, 0.0);
+        }
+        for i in one.iter() {
+            assert_f64_near!(*i, 1.0);
+        }
+
+        let arr = GMatrix2x2::from_nested_arr([[1.0, 2.0], [3.0, 4.0]]);
+        assert_f64_near!(arr[[0, 0]], 1.0);
+        assert_f64_near!(arr[[1, 1]], 4.0);
+    }
+
+    #[test]
+    fn test_constructors_complex() {
+        let c1 = Complex {
+            real: 1.0,
+            imaginary: 2.0
+        };
+        let c2 = Complex {
+            real: 3.0,
+            imaginary: 4.0
+        };
+        let mat = GMatrix::<1, 2, Complex>::from_nested_arr([[c1, c2]]);
+
+        assert_eq!(mat[[0, 0]], c1);
+        assert_eq!(mat[[0, 1]], c2);
+    }
+
+    #[test]
+    fn test_shape_properties() {
+        assert!(GMatrix2x2::<f64>::is_square());
+        assert!(!GMatrix::<2, 3, f64>::is_square());
+
+        let mut mat = GMatrix2x2::<f64>::zeros();
+        mat[[0, 0]] = 1.0;
+        mat[[1, 1]] = 1.0;
+        assert_eq!(mat.count_nonzero(), 2);
+    }
+
+    #[test]
+    fn test_triangular_checks() {
+        let ut = GMatrix2x2::<f64>::from_nested_arr([[1.0, 2.0], [0.0, 3.0]]);
+        assert!(ut.is_upper_triangular());
+        assert!(!ut.is_lower_triangular());
+
+        let lt = GMatrix2x2::<f64>::from_nested_arr([[1.0, 0.0], [2.0, 3.0]]);
+        assert!(lt.is_lower_triangular());
+        assert!(!lt.is_upper_triangular());
+
+        let diag = GMatrix2x2::<f64>::from_nested_arr([[1.0, 0.0], [0.0, 3.0]]);
+        assert!(diag.is_diagonal());
+    }
+
+    #[test]
+    fn test_iterators() {
+        let mut mat = GMatrix2x2::<f64>::from_nested_arr([[1.0, 2.0], [3.0, 4.0]]);
+
+        let indexed: Vec<((usize, usize), f64)> =
+            mat.iter_indexed().map(|(pos, &val)| (pos, val)).collect();
+        assert_eq!(indexed[1], ((0, 1), 2.0));
+
+        // Mutable iteration
+        for ((r, c), val) in mat.iter_indexed_mut() {
+            if r == c {
+                *val = 0.0;
+            }
+        }
+        assert_f64_near!(mat[[0, 0]], 0.0);
+        assert_f64_near!(mat[[1, 1]], 0.0);
+    }
+
+    #[test]
+    fn test_transpose() {
+        // Test Vector (Zero-Copy path)
+        let vec = GMatrix::<1, 3, f64>::from_nested_arr([[1.0, 2.0, 3.0]]);
+        let vec_t = vec.transpose();
+        assert_f64_near!(vec_t[[0, 0]], 1.0);
+        assert_f64_near!(vec_t[[2, 0]], 3.0);
+
+        // Test Matrix (Standard path)
+        let mat = GMatrix2x2::<f64>::from_nested_arr([[1.0, 2.0], [3.0, 4.0]]);
+        let mat_t = mat.transpose();
+        assert_f64_near!(mat_t[[0, 1]], 3.0);
+        assert_f64_near!(mat_t[[1, 0]], 2.0);
+    }
+
+    #[test]
+    #[should_panic(expected = "out of bounds")]
+    fn test_index_out_of_bounds() {
+        let mat = GMatrix2x2::<f64>::zeros();
+        let _ = mat[[2, 0]];
+    }
+
+    #[test]
+    #[cfg(feature = "std")]
+    fn test_diagonals() {
+        let mat =
+            GMatrix3x3::<f64>::from_nested_arr([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 3.0]]);
+        assert_eq!(mat.diagonals(), vec![1.0, 2.0, 3.0]);
+    }
+
+    #[test]
+    fn test_determinant_2x2() {
+        let m = GMatrix::<2, 2, f64> {
+            values: [1.0, 2.0, 3.0, 4.0]
+        };
+        assert_f64_near!(m.determinant(), -2.0);
+    }
+
+    #[test]
+    fn test_determinant_3x3() {
+        let m = GMatrix::<3, 3, f64> {
+            values: [6.0, 1.0, 1.0, 4.0, -2.0, 5.0, 2.0, 8.0, 7.0]
+        };
+        assert_f64_near!(m.determinant(), -306.0);
+    }
+
+    #[test]
+    fn test_determinant_4x4_lu() {
+        let m = GMatrix::<4, 4, f64> {
+            values: [
+                1.0, 3.0, 5.0, 9.0, 1.0, 3.0, 1.0, 7.0, 4.0, 3.0, 9.0, 7.0, 5.0, 2.0, 0.0, 9.0
+            ]
+        };
+        assert_f64_near!(m.determinant(), -376.0);
+    }
+
+    #[test]
+    fn test_determinant_singular() {
+        let m = GMatrix::<3, 3, f64> {
+            values: [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]
+        };
+        assert_f64_near!(m.determinant(), 0.0);
+    }
+
+    #[test]
+    fn test_determinant_identity() {
+        let m = GMatrix::<4, 4, f64> {
+            values: [
+                1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0
+            ]
+        };
+        assert_f64_near!(m.determinant(), 1.0);
+    }
+
+    #[test]
+    #[should_panic(expected = "Cannot take the determinant of a non square matrix")]
+    fn test_determinant_non_square_panic() {
+        let m = GMatrix::<2, 3, f64> {
+            values: [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
+        };
+        let _ = m.determinant();
+    }
+
+    #[test]
+    fn test_minor_submatrix_extraction() {
+        let m = GMatrix::<3, 3, f64> {
+            values: [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]
+        };
+
+        let sub = m.get_minor_submatrix::<2, 2>(0, 1);
+        let expected = [4.0, 6.0, 7.0, 9.0];
+
+        for (i, j) in sub.values.into_iter().zip(expected.into_iter()) {
+            assert_f64_near!(i, j);
+        }
+    }
+
+    #[test]
+    fn test_cofactor_values() {
+        let m = GMatrix::<3, 3, f64> {
+            values: [1.0, 3.0, 2.0, 4.0, 5.0, 0.0, 2.0, 1.0, 2.0]
+        };
+
+        assert_f64_near!(m.cofactor(1, 1), -2.0);
+        assert_f64_near!(m.cofactor(0, 1), -8.0);
+    }
+
+    #[test]
+    fn test_full_cofactor_matrix() {
+        let m = GMatrix::<3, 3, f64> {
+            values: [1.0, 2.0, 3.0, 0.0, 4.0, 5.0, 1.0, 0.0, 6.0]
+        };
+
+        let cofactors = m.cofactor_matrix();
+        let expected = [24.0, 5.0, -4.0, -12.0, 3.0, 2.0, -2.0, -5.0, 4.0];
+
+        for (i, j) in cofactors.values.into_iter().zip(expected.into_iter()) {
+            assert_f64_near!(i, j);
+        }
+    }
+
+    #[test]
+    fn test_cofactor_determinant_consistency() {
+        let m = GMatrix::<3, 3, f64> {
+            values: [3.0, 5.0, 0.0, 2.0, -1.0, 4.0, 6.0, 0.0, 2.0]
+        };
+
+        let det_direct = m.determinant();
+
+        let det_expansion = m.values[0] * m.cofactor(0, 0)
+            + m.values[1] * m.cofactor(0, 1)
+            + m.values[2] * m.cofactor(0, 2);
+
+        assert_f64_near!(det_direct, det_expansion);
+    }
+
+    #[test]
+    #[should_panic(expected = "Wrong sized generics")]
+    fn test_minor_wrong_dimensions_panic() {
+        let m = GMatrix::<3, 3, f64> { values: [0.0; 9] };
+        let _ = m.get_minor_submatrix::<3, 3>(0, 0);
+    }
+
+    #[test]
+    fn test_adjoint_2x2() {
+        let m = GMatrix::<2, 2, f64> {
+            values: [1.0, 2.0, 3.0, 4.0]
+        };
+        let adj = m.adjoint();
+        let expected = [4.0, -2.0, -3.0, 1.0];
+
+        for (i, j) in adj.values.into_iter().zip(expected.into_iter()) {
+            assert_f64_near!(i, j);
+        }
+    }
+
+    #[test]
+    fn test_adjoint_3x3() {
+        let m = GMatrix::<3, 3, f64> {
+            values: [1.0, 2.0, 3.0, 0.0, 1.0, 4.0, 5.0, 6.0, 0.0]
+        };
+
+        let adj = m.adjoint();
+        let expected = [-24.0, 18.0, 5.0, 20.0, -15.0, -4.0, -5.0, 4.0, 1.0];
+
+        for (i, j) in adj.values.into_iter().zip(expected.into_iter()) {
+            assert_f64_near!(i, j);
+        }
+    }
+
+    #[test]
+    fn test_adjoint_inverse_relationship() {
+        let m = GMatrix::<3, 3, f64> {
+            values: [1.0, 2.0, 3.0, 0.0, 4.0, 5.0, 1.0, 0.0, 6.0]
+        };
+
+        let det = m.determinant();
+        let adj = m.adjoint();
+
+        let row0_col0 = (m.values[0] * adj.values[0]
+            + m.values[1] * adj.values[3]
+            + m.values[2] * adj.values[6]);
+
+        assert_f64_near!(row0_col0, det);
+    }
+
+    #[test]
+    fn test_adjoint_identity() {
+        let m = GMatrix::<3, 3, f64> {
+            values: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
+        };
+        let adj = m.adjoint();
+        let expected = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0];
+
+        for (i, j) in adj.values.into_iter().zip(expected.into_iter()) {
+            assert_f64_near!(i, j);
+        }
+    }
+}
diff --git a/vec-utils/src/matrix/real.rs b/vec-utils/src/matrix/real.rs
new file mode 100644
index 0000000..f5f92c5
--- /dev/null
+++ b/vec-utils/src/matrix/real.rs
@@ -0,0 +1,65 @@
+use core::ops::Mul;
+
+use matrixmultiply::dgemm;
+
+#[doc(inline)]
+use crate::matrix::generic::GMatrix;
+
+/// A generic 2d matrix of width R and height C
+pub type Matrix<const R: usize, const C: usize> = GMatrix<R, C, f64>;
+
+/// An alias for 2x2 matracies
+pub type Matrix2x2 = Matrix<2, 2>;
+/// An alias for 3x3 matracies
+pub type Matrix3x3 = Matrix<3, 3>;
+
+impl<const R: usize, const C: usize> Matrix<R, C> where [f64; R * C]: Sized {}
+
+impl<const R: usize, const C: usize, const U: usize> Mul<Matrix<U, C>> for Matrix<R, U>
+where
+    [f64; R * C]: Sized,
+    [f64; R * U]: Sized,
+    [f64; U * C]: Sized
+{
+    type Output = Matrix<R, C>;
+
+    fn mul(self, rhs: Matrix<U, C>) -> Self::Output {
+        let mut result = Matrix::<R, C>::zeros();
+        // Safety: dgemm is an unsafe function
+        unsafe {
+            dgemm(
+                R,
+                U,
+                C,
+                1.0,
+                self.values.as_ptr(),
+                U.cast_signed(),
+                1,
+                rhs.values.as_ptr(),
+                C.cast_signed(),
+                1,
+                0.0,
+                result.values.as_mut_ptr(),
+                C.cast_signed(),
+                1
+            );
+        }
+        result
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use pretty_assertions::assert_eq;
+
+    use super::*;
+
+    #[test]
+    fn test_mul() {
+        let lhs = Matrix3x3::from_nested_arr([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]);
+        let rhs = Matrix::from_nested_arr([[10.0], [11.0], [12.0]]);
+        let result = lhs * rhs;
+        let correct = Matrix::from_nested_arr([[68.0], [167.0], [266.0]]);
+        assert_eq!(result, correct);
+    }
+}
diff --git a/vec-utils/src/matrix/traits.rs b/vec-utils/src/matrix/traits.rs
new file mode 100644
index 0000000..755346f
--- /dev/null
+++ b/vec-utils/src/matrix/traits.rs
@@ -0,0 +1,88 @@
+use crate::complex::Complex;
+
+pub trait Zeroable {
+    fn zero() -> Self;
+
+    fn is_zero(&self) -> bool;
+}
+
+pub trait Oneable {
+    fn one() -> Self;
+
+    fn is_one(&self) -> bool;
+}
+
+pub trait Signed {
+    fn abs(&self) -> Self;
+
+    fn flip(&mut self);
+}
+
+impl Zeroable for f64 {
+    fn is_zero(&self) -> bool {
+        *self == 0.0
+    }
+
+    fn zero() -> Self {
+        0.0
+    }
+}
+
+impl Oneable for f64 {
+    fn is_one(&self) -> bool {
+        (self - 1.0).abs() < f64::EPSILON
+    }
+
+    fn one() -> Self {
+        1.0
+    }
+}
+
+impl Signed for f64 {
+    fn abs(&self) -> Self {
+        f64::abs(*self)
+    }
+
+    fn flip(&mut self) {
+        *self *= -1.0;
+    }
+}
+
+impl Zeroable for Complex {
+    fn is_zero(&self) -> bool {
+        self.real.abs() < f64::EPSILON && self.imaginary.abs() < f64::EPSILON
+    }
+
+    fn zero() -> Self {
+        Self {
+            real: 0.0,
+            imaginary: 0.0
+        }
+    }
+}
+
+impl Oneable for Complex {
+    fn is_one(&self) -> bool {
+        (self.real - 1.0).abs() < f64::EPSILON && self.imaginary.abs() < f64::EPSILON
+    }
+
+    fn one() -> Self {
+        Self {
+            real: 1.0,
+            imaginary: 0.0
+        }
+    }
+}
+
+impl Signed for Complex {
+    fn abs(&self) -> Self {
+        Self {
+            real: self.magnitude(),
+            imaginary: 0.0
+        }
+    }
+
+    fn flip(&mut self) {
+        *self = *self * -1.0;
+    }
+}
diff --git a/vec-utils/src/quat.rs b/vec-utils/src/quat.rs
index 663bd4c..c61fede 100644
--- a/vec-utils/src/quat.rs
+++ b/vec-utils/src/quat.rs
@@ -4,6 +4,8 @@ use core::ops::{Add, Div, Index, Mul, Sub};
 use std::vec::Vec;
 
 use crate::angle::AngleRadians;
+#[cfg(feature = "matrix")]
+use crate::matrix::real::Matrix3x3;
 use crate::vec3d::Vec3d;
 use crate::{
     impl_dual_op_variants, impl_single_op_comm, impl_single_op_variants,
@@ -62,49 +64,50 @@ impl Quat {
     }
 
     /// Create a new quaternion from a rotation matrix
-    pub fn from_rotation_matrix(m: &[[f64; 3]; 3]) -> Quat {
+    #[cfg(feature = "matrix")]
+    pub fn from_rotation_matrix(m: &Matrix3x3) -> Quat {
         #[cfg(not(feature = "std"))]
-        let w = core::f64::math::sqrt(1.0 + m[0][0] + m[1][1] + m[2][2]) / 2.0;
+        let w = core::f64::math::sqrt(1.0 + m[[0, 0]] + m[[1, 1]] + m[[2, 2]]) / 2.0;
         #[cfg(feature = "std")]
-        let w = (1.0 + m[0][0] + m[1][1] + m[2][2]).sqrt() / 2.0;
+        let w = (1.0 + m[[0, 0]] + m[[1, 1]] + m[[2, 2]]).sqrt() / 2.0;
         #[cfg(not(feature = "std"))]
-        let i = core::f64::math::sqrt(1.0 + m[0][0] - m[1][1] - m[2][2]) / 2.0;
+        let i = core::f64::math::sqrt(1.0 + m[[0, 0]] - m[[1, 1]] - m[[2, 2]]) / 2.0;
         #[cfg(feature = "std")]
-        let i = (1.0 + m[0][0] - m[1][1] - m[2][2]).sqrt() / 2.0;
+        let i = (1.0 + m[[0, 0]] - m[[1, 1]] - m[[2, 2]]).sqrt() / 2.0;
         #[cfg(not(feature = "std"))]
-        let j = core::f64::math::sqrt(1.0 - m[0][0] + m[1][1] - m[2][2]) / 2.0;
+        let j = core::f64::math::sqrt(1.0 - m[[0, 0]] + m[[1, 1]] - m[[2, 2]]) / 2.0;
         #[cfg(feature = "std")]
-        let j = (1.0 - m[0][0] + m[1][1] - m[2][2]).sqrt() / 2.0;
+        let j = (1.0 - m[[0, 0]] + m[[1, 1]] - m[[2, 2]]).sqrt() / 2.0;
         #[cfg(not(feature = "std"))]
-        let k = core::f64::math::sqrt(1.0 - m[0][0] - m[1][1] + m[2][2]) / 2.0;
+        let k = core::f64::math::sqrt(1.0 - m[[0, 0]] - m[[1, 1]] + m[[2, 2]]) / 2.0;
         #[cfg(feature = "std")]
-        let k = (1.0 - m[0][0] - m[1][1] + m[2][2]).sqrt() / 2.0;
+        let k = (1.0 - m[[0, 0]] - m[[1, 1]] + m[[2, 2]]).sqrt() / 2.0;
         if w > i && w > j && w > k {
             Quat {
                 w,
-                i: (m[2][1] - m[1][2]) / (4.0 * w),
-                j: (m[0][2] - m[2][0]) / (4.0 * w),
-                k: (m[1][0] - m[0][1]) / (4.0 * w)
+                i: (m[[2, 1]] - m[[1, 2]]) / (4.0 * w),
+                j: (m[[0, 2]] - m[[2, 0]]) / (4.0 * w),
+                k: (m[[1, 0]] - m[[0, 1]]) / (4.0 * w)
             }
         } else if i > j && i > k {
             Quat {
-                w: (m[2][1] - m[1][2]) / (4.0 * i),
+                w: (m[[2, 1]] - m[[1, 2]]) / (4.0 * i),
                 i,
-                j: (m[0][1] + m[1][0]) / (4.0 * i),
-                k: (m[0][2] + m[2][0]) / (4.0 * i)
+                j: (m[[0, 1]] + m[[1, 0]]) / (4.0 * i),
+                k: (m[[0, 2]] + m[[2, 0]]) / (4.0 * i)
             }
         } else if j > k {
             Quat {
-                w: (m[0][2] - m[2][0]) / (4.0 * j),
-                i: (m[0][1] + m[1][0]) / (4.0 * j),
+                w: (m[[0, 2]] - m[[2, 0]]) / (4.0 * j),
+                i: (m[[0, 1]] + m[[1, 0]]) / (4.0 * j),
                 j,
-                k: (m[1][2] + m[2][1]) / (4.0 * j)
+                k: (m[[1, 2]] + m[[2, 1]]) / (4.0 * j)
             }
         } else {
             Quat {
-                w: (m[1][0] - m[0][1]) / (4.0 * k),
-                i: (m[0][2] + m[2][0]) / (4.0 * k),
-                j: (m[1][2] + m[2][1]) / (4.0 * k),
+                w: (m[[1, 0]] - m[[0, 1]]) / (4.0 * k),
+                i: (m[[0, 2]] + m[[2, 0]]) / (4.0 * k),
+                j: (m[[1, 2]] + m[[2, 1]]) / (4.0 * k),
                 k
             }
         }
@@ -175,8 +178,9 @@ impl Quat {
     }
 
     /// Convert the quaternion to a rotation matrix
-    pub fn to_rotation_matrix(&self) -> [[f64; 3]; 3] {
-        [
+    #[cfg(feature = "matrix")]
+    pub fn to_rotation_matrix(&self) -> Matrix3x3 {
+        Matrix3x3::from_nested_arr([
             [
                 1.0 - 2.0 * (self.j * self.j + self.k * self.k),
                 2.0 * (self.i * self.j - self.k * self.w),
@@ -192,7 +196,7 @@ impl Quat {
                 2.0 * (self.j * self.k + self.i * self.w),
                 1.0 - 2.0 * (self.i * self.i + self.j * self.j)
             ]
-        ]
+        ])
     }
 
     /// Rotate a vector by the quaternion
@@ -340,8 +344,9 @@ mod tests {
     }
 
     #[test]
+    #[cfg(feature = "matrix")]
     fn test_from_rotation_matrix() {
-        let m = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
+        let m = Matrix3x3::from_nested_arr([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]);
         let q = Quat::from_rotation_matrix(&m);
         let good = Quat::new(1.0, 0.0, 0.0, 0.0);
         assert_eq!(q, good);
@@ -385,18 +390,19 @@ mod tests {
     }
 
     #[test]
+    #[cfg(feature = "matrix")]
     fn test_to_rotation_matrix() {
         let q = Quat::new(1.0, 0.0, 0.0, 0.0);
         let m = q.to_rotation_matrix();
-        assert_f64_near!(m[0][0], 1.0);
-        assert_f64_near!(m[0][1], 0.0);
-        assert_f64_near!(m[0][2], 0.0);
-        assert_f64_near!(m[1][0], 0.0);
-        assert_f64_near!(m[1][1], 1.0);
-        assert_f64_near!(m[1][2], 0.0);
-        assert_f64_near!(m[2][0], 0.0);
-        assert_f64_near!(m[2][1], 0.0);
-        assert_f64_near!(m[2][2], 1.0);
+        assert_f64_near!(m[[0, 0]], 1.0);
+        assert_f64_near!(m[[0, 1]], 0.0);
+        assert_f64_near!(m[[0, 2]], 0.0);
+        assert_f64_near!(m[[1, 0]], 0.0);
+        assert_f64_near!(m[[1, 1]], 1.0);
+        assert_f64_near!(m[[1, 2]], 0.0);
+        assert_f64_near!(m[[2, 0]], 0.0);
+        assert_f64_near!(m[[2, 1]], 0.0);
+        assert_f64_near!(m[[2, 2]], 1.0);
     }
 
     #[test]
diff --git a/vec-utils/src/vec3d.rs b/vec-utils/src/vec3d.rs
index 2420227..de02a43 100644
--- a/vec-utils/src/vec3d.rs
+++ b/vec-utils/src/vec3d.rs
@@ -4,6 +4,8 @@ use core::{f64, fmt};
 use std::vec::Vec;
 
 use crate::angle::AngleRadians;
+#[cfg(feature = "matrix")]
+use crate::matrix::real::Matrix;
 use crate::quat::Quat;
 use crate::{
     impl_dual_op_variants, impl_single_op_comm, impl_single_op_variants,
@@ -122,6 +124,22 @@ impl Vec3d {
         vec![self.x, self.y, self.z]
     }
 
+    /// Convert the Vec3d to a 1x3 matrix
+    #[cfg(feature = "matrix")]
+    pub fn to_hmatrix(&self) -> Matrix<1, 3> {
+        Matrix::<1, 3> {
+            values: [self.x, self.y, self.z]
+        }
+    }
+
+    /// Convert the Vec3d to a 3x1 matrix
+    #[cfg(feature = "matrix")]
+    pub fn to_vmatrix(&self) -> Matrix<3, 1> {
+        Matrix::<3, 1> {
+            values: [self.x, self.y, self.z]
+        }
+    }
+
     /// Calculate the dot product of two Vec3d
     pub fn dot(&self, other: &Vec3d) -> f64 {
         self.x * other.x + self.y * other.y + self.z * other.z
@@ -416,6 +434,18 @@ mod tests {
         assert_eq!(vec, good);
     }
 
+    #[test]
+    #[cfg(feature = "matrix")]
+    fn test_to_matrix() {
+        let vec = Vec3d::new(1.0, 2.0, 3.0);
+        let hmat = vec.to_hmatrix();
+        let vmat = vec.to_vmatrix();
+        let hgood = Matrix::<1, 3>::from_nested_arr([[1.0, 2.0, 3.0]]);
+        let vgood = hgood.transpose();
+        assert_eq!(hmat, hgood);
+        assert_eq!(vmat, vgood);
+    }
+
     #[test]
     fn test_dot() {
         let v1 = Vec3d::new(1.0, 2.0, 3.0);