diff --git a/.vscode/settings.json b/.vscode/settings.json
index ad6aaf8..c1cd6c2 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -69,7 +69,8 @@
         "streambuf": "cpp",
         "thread": "cpp",
         "typeinfo": "cpp",
-        "variant": "cpp"
+        "variant": "cpp",
+        "shared_mutex": "cpp"
     },
     "clangd.enable": false
 }
\ No newline at end of file
diff --git a/Matrix.hpp b/Matrix.hpp
index 3fb4848..3a4bea2 100644
--- a/Matrix.hpp
+++ b/Matrix.hpp
@@ -15,6 +15,7 @@ public:
    */
   Matrix(const std::array<float, rows * columns> &array);
 
+  Matrix(const Matrix<rows, columns> &other);
   // TODO: Figure out how to do this
   /**
    * @brief Initialize a matrix directly with any number of arguments
@@ -28,8 +29,8 @@ public:
    * @param result A buffer to store the result into
    * @note there is no problem if result == this
    */
-  void Add(const Matrix<rows, columns> &other,
-           Matrix<rows, columns> &result) const;
+  Matrix<rows, columns> &Add(const Matrix<rows, columns> &other,
+                             Matrix<rows, columns> &result) const;
 
   /**
    * @brief Element-wise subtract matrix
@@ -37,8 +38,8 @@ public:
    * @param result A buffer to store the result into
    * @note there is no problem if result == this
    */
-  void Sub(const Matrix<rows, columns> &other,
-           Matrix<rows, columns> &result) const;
+  Matrix<rows, columns> &Sub(const Matrix<rows, columns> &other,
+                             Matrix<rows, columns> &result) const;
 
   /**
    * @brief Matrix multiply the two matrices
@@ -46,8 +47,8 @@ public:
    * @param result A buffer to store the result into
    */
   template <uint8_t other_columns>
-  void Mult(const Matrix<rows, columns> &other,
-            Matrix<columns, other_columns> &result) const;
+  Matrix<rows, columns> &Mult(const Matrix<columns, other_columns> &other,
+                              Matrix<rows, other_columns> &result) const;
 
   /**
    * @brief Multiply the matrix by a scalar
@@ -55,37 +56,14 @@ public:
    * @param result A buffer to store the result into
    * @note there is no problem if result == this
    */
-  void Mult(float scalar, Matrix<rows, columns> &result) const;
-
-  /**
-   * @brief Invert this matrix
-   * @param result A buffer to store the result into
-   * @warning this is super slow! Only call it if you absolutely have to!!!
-   */
-  void Invert(Matrix<rows, columns> &result) const;
-
-  /**
-   * @brief Transpose this matrix
-   * @param result A buffer to store the result into
-   */
-  void Transpose(Matrix<columns, rows> &result) const;
+  Matrix<rows, columns> &Mult(float scalar,
+                              Matrix<rows, columns> &result) const;
 
   /**
    * @brief Square this matrix
    * @param result A buffer to store the result into
    */
-  void Square(Matrix<rows, columns> &result) const;
-
-  /**
-   * @return Get the determinant of the matrix
-   */
-  float Det() const;
-
-  /**
-   * @brief Calculate the eigenvalues for a square matrix
-   * @param result a buffer to store the result into
-   */
-  void EigenValues(Matrix<rows, 1> &eigenvalues) const;
+  Matrix<rows, columns> &Square(Matrix<rows, rows> &result) const;
 
   /**
    * @brief Element-wise multiply the two matrices
@@ -93,8 +71,8 @@ public:
    * @param result A buffer to store the result into
    * @note there is no problem if result == this
    */
-  void ElementMultiply(const Matrix<rows, columns> &other,
-                       Matrix<rows, columns> &result) const;
+  Matrix<rows, columns> &ElementMultiply(const Matrix<rows, columns> &other,
+                                         Matrix<rows, columns> &result) const;
 
   /**
    * @brief Element-wise divide the two matrices
@@ -102,8 +80,59 @@ public:
    * @param result A buffer to store the result into
    * @note there is no problem if result == this
    */
-  void ElementDivide(const Matrix<rows, columns> &other,
-                     Matrix<rows, columns> &result) const;
+  Matrix<rows, columns> &ElementDivide(const Matrix<rows, columns> &other,
+                                       Matrix<rows, columns> &result) const;
+  /**
+   * @return Get the determinant of the matrix
+   * @note for right now only 2x2 and 3x3 matrices are supported
+   */
+  float Det() const;
+
+  /**
+   * @brief Invert this matrix
+   * @param result A buffer to store the result into
+   * @warning this is super slow! Only call it if you absolutely have to!!!
+   */
+  Matrix<rows, columns> &Invert(Matrix<rows, columns> &result) const;
+
+  /**
+   * @brief Transpose this matrix
+   * @param result A buffer to store the result into
+   */
+  Matrix<columns, rows> &Transpose(Matrix<columns, rows> &result) const;
+
+  /**
+   * @brief reduce the matrix so the sum of its elements equal 1
+   * @param result a buffer to store the result into
+   */
+  Matrix<rows, columns> &Normalize(Matrix<rows, columns> &result) const;
+
+  /**
+   * @brief Get a row from the matrix
+   * @param row_index the row index to get
+   * @param row a buffer to write the row into
+   */
+  Matrix<1, columns> &GetRow(uint8_t row_index, Matrix<1, columns> &row) const;
+
+  /**
+   * @brief Get a row from the matrix
+   * @param column_index the row index to get
+   * @param column a buffer to write the row into
+   */
+  Matrix<rows, 1> &GetColumn(uint8_t column_index,
+                             Matrix<rows, 1> &column) const;
+
+  /**
+   * @brief Get the number of rows in this matrix
+   */
+  constexpr uint8_t GetRowSize() { return rows; }
+
+  /**
+   * @brief Get the number of columns in this matrix
+   */
+  constexpr uint8_t GetColumnSize() { return columns; }
+
+  void ToString(std::string &stringBuffer) const;
 
   /**
    * @brief Get an element from the matrix
@@ -118,6 +147,10 @@ public:
    * internal array
    */
   std::array<float, columns> &operator[](uint8_t row_index) {
+    if (row_index > rows - 1) {
+      return this->matrix[0]; // TODO: We should throw something here instead of
+                              // failing quietly.
+    }
     return this->matrix[row_index];
   }
 
@@ -127,34 +160,10 @@ public:
         this->matrix[row_idx][column_idx] = other.Get(row_idx, column_idx);
       }
     }
+    // return a reference to ourselves so you can chain together these functions
+    return *this;
   }
 
-  /**
-   * @brief Get a row from the matrix
-   * @param row_index the row index to get
-   * @param row a buffer to write the row into
-   */
-  void GetRow(uint8_t row_index, Matrix<1, columns> &row) const;
-
-  /**
-   * @brief Get a row from the matrix
-   * @param column_index the row index to get
-   * @param column a buffer to write the row into
-   */
-  void GetColumn(uint8_t column_index, Matrix<rows, 1> &column) const;
-
-  /**
-   * @brief Get the number of rows in this matrix
-   */
-  constexpr uint8_t GetRowSize() { return rows; }
-
-  /**
-   * @brief Get the number of columns in this matrix
-   */
-  constexpr uint8_t GetColumnSize() { return columns; }
-
-  void ToString(std::string &stringBuffer) const;
-
 private:
   /**
    * @brief take the dot product of the two vectors
@@ -163,25 +172,21 @@ private:
   static float dotProduct(const Matrix<1, vector_size> &vec1,
                           const Matrix<1, vector_size> &vec2);
 
+  template <uint8_t vector_size>
+  static float dotProduct(const Matrix<vector_size, 1> &vec1,
+                          const Matrix<vector_size, 1> &vec2);
   /**
    * @brief Set all elements in this matrix to zero
    */
   void zeroMatrix();
 
-  void matrixOfMinors(Matrix<rows, columns> &result) const;
+  Matrix<rows, columns> &matrixOfMinors(Matrix<rows, columns> &result) const;
 
-  void minorMatrix(Matrix<rows - 1, columns - 1> &result, uint8_t row_idx,
-                   uint8_t column_idx) const;
+  Matrix<rows - 1, columns - 1> &
+      minorMatrix(Matrix<rows - 1, columns - 1> &result, uint8_t row_idx,
+                  uint8_t column_idx) const;
 
-  void adjugate(Matrix<rows, columns> &result) const;
-
-  /**
-   * @brief reduce the matrix so the sum of its elements equal 1
-   * @param result a buffer to store the result into
-   */
-  void normalize(Matrix<rows, columns> &result) const;
-
-  constexpr bool isSquare() { return rows == columns; }
+  Matrix<rows, columns> &adjugate(Matrix<rows, columns> &result) const;
 
   void setMatrixToArray(const std::array<float, rows * columns> &array);
 
@@ -232,31 +237,46 @@ Matrix<rows, columns>::Matrix(const std::array<float, rows * columns> &array) {
 // }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::Add(const Matrix<rows, columns> &other,
-                                Matrix<rows, columns> &result) const {
+Matrix<rows, columns>::Matrix(const Matrix<rows, columns> &other) {
+  for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
+    for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
+      this->matrix[row_idx][column_idx] = other.Get(row_idx, column_idx);
+    }
+  }
+}
+
+template <uint8_t rows, uint8_t columns>
+Matrix<rows, columns> &
+Matrix<rows, columns>::Add(const Matrix<rows, columns> &other,
+                           Matrix<rows, columns> &result) const {
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
     for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
       result[row_idx][column_idx] =
           this->Get(row_idx, column_idx) + other.Get(row_idx, column_idx);
     }
   }
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::Sub(const Matrix<rows, columns> &other,
-                                Matrix<rows, columns> &result) const {
+Matrix<rows, columns> &
+Matrix<rows, columns>::Sub(const Matrix<rows, columns> &other,
+                           Matrix<rows, columns> &result) const {
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
     for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
       result[row_idx][column_idx] =
           this->Get(row_idx, column_idx) - other.Get(row_idx, column_idx);
     }
   }
+
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
 template <uint8_t other_columns>
-void Matrix<rows, columns>::Mult(const Matrix<rows, columns> &other,
-                                 Matrix<columns, other_columns> &result) const {
+Matrix<rows, columns> &
+Matrix<rows, columns>::Mult(const Matrix<columns, other_columns> &other,
+                            Matrix<rows, other_columns> &result) const {
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
     for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
       // get our row
@@ -274,20 +294,25 @@ void Matrix<rows, columns>::Mult(const Matrix<rows, columns> &other,
           Matrix<rows, columns>::dotProduct(this_row, other_column_t);
     }
   }
+
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::Mult(float scalar,
-                                 Matrix<rows, columns> &result) const {
+Matrix<rows, columns> &
+Matrix<rows, columns>::Mult(float scalar, Matrix<rows, columns> &result) const {
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
     for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
       result[row_idx][column_idx] = this->Get(row_idx, column_idx) * scalar;
     }
   }
+
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
+Matrix<rows, columns> &
+Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
   // since all matrix sizes have to be statically specified at compile time we
   // can do this
   static_assert(rows == columns,
@@ -295,10 +320,11 @@ void Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
 
   // unfortunately we can't calculate this at compile time so we'll just reurn
   // zeros
+  float determinant{this->Det()};
   if (this->Det() < 0) {
     // you can't invert a matrix with a negative determinant
     result.zeroMatrix();
-    return;
+    return result;
   }
 
   // TODO: This algorithm is really inneficient because of the matrix of minors.
@@ -311,138 +337,139 @@ void Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
 
   // now adjugate the matrix and save it in our output
   minors.adjugate(result);
-  float determinant = this->Det();
 
   // scale the result by 1/determinant and we have our answer
   result.Mult(1 / determinant);
+
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::Transpose(Matrix<columns, rows> &result) const {
+Matrix<columns, rows> &
+Matrix<rows, columns>::Transpose(Matrix<columns, rows> &result) const {
   for (uint8_t column_idx{0}; column_idx < rows; column_idx++) {
     for (uint8_t row_idx{0}; row_idx < columns; row_idx++) {
       result[row_idx][column_idx] = this->Get(column_idx, row_idx);
     }
   }
+
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::Square(Matrix<rows, columns> &result) const {
-  static_assert(this->isSquare(), "You can't square an non-square matrix.");
+Matrix<rows, columns> &
+Matrix<rows, columns>::Square(Matrix<rows, rows> &result) const {
+  // TODO: Because template requirements are checked before static_assert, this
+  // never throws an error and fails at the Mult call instead.
+  static_assert(rows == columns, "You can't square a non-square matrix.");
 
-  this->Mult(this, result);
+  this->Mult(*this, result);
+
+  return result;
+}
+
+// explicitly define the determinant for a 3x3 matrix because it is definitely
+// the fastest way to calculte a 2x2 matrix determinant
+template <> float Matrix<2, 2>::Det() const {
+  return this->matrix[0][0] * this->matrix[1][1] -
+         this->matrix[0][1] * this->matrix[1][1];
+}
+
+// explicitly define the determinant for a 3x3 matrix because it will probably
+// be faster than the jacobi method for nxn matrices
+template <> float Matrix<3, 3>::Det() const {
+  float a{this->matrix[0][0]};
+  float b{this->matrix[0][1]};
+  float c{this->matrix[0][2]};
+
+  Matrix<2, 2> minors{};
+  this->minorMatrix(minors, 0, 0);
+  float det = a * minors.Det();
+
+  this->minorMatrix(minors, 0, 1);
+  det -= b * minors.Det();
+
+  this->minorMatrix(minors, 0, 2);
+  det += c * minors.Det();
+
+  return det;
 }
 
 template <uint8_t rows, uint8_t columns>
 float Matrix<rows, columns>::Det() const {
-  static_assert(this->isSquare(),
-                "You can't take the determinant of a non-square matrix.");
-  Matrix<1, columns> eigenValues{};
-  this->EigenValues(eigenValues);
-
-  float determinant{1};
-  for (uint8_t i{0}; i < columns; i++) {
-    determinant *= eigenValues.Get(0, i);
-  }
-
-  return determinant;
-}
-
-template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::EigenValues(Matrix<rows, 1> &eigenvalues) const {
   static_assert(rows == columns,
-                "Eigenvalues can only be computed for square matrices.");
-  // I got this code from:
-  // https://www.quora.com/What-is-the-C-code-for-finding-eigenvalues
-  Matrix<rows, 1> v{};
-  Matrix<rows, 1> Av{};
-  Matrix<rows, 1> z{};
+                "You can't take the determinant of a non-square matrix.");
+  // static_assert(
+  //     false,
+  //     "Right now this operation isn't supported for matrices bigger than
+  //     3x3");
+  // Matrix<1, columns> eigenValues{};
+  // this->EigenValues(eigenValues);
 
-  float d = 0.0;
-  float d_old = 0.0;
-  constexpr float convergence_value{1e-6};
-  constexpr uint16_t max_iterations{500};
+  // float determinant{1};
+  // for (uint8_t i{0}; i < columns; i++) {
+  //   determinant *= eigenValues.Get(0, i);
+  // }
 
-  // Initialize v as a random vector
-  for (int i = 0; i < rows; i++) {
-    v[0][i] = rand() / RAND_MAX;
-  }
-
-  // run this loop until the eigenvalues converge or we give up
-  for (uint16_t k{0}; k < max_iterations; k++) {
-    /* Multiply A by v */
-    for (int i = 0; i < rows; i++) {
-      Av[0][i] = 0.0;
-      for (int j = 0; j < rows; j++) {
-        Av[0][i] += this->Get(0, i * rows + j) * v[0][j];
-      }
-    }
-
-    // Calculate the eigenvalue and update v
-    d_old = d;
-    d = dot_product(v, Av, rows);
-    for (int i = 0; i < rows; i++) {
-      z[0][i] = Av[0][i] - d * v[0][i];
-    }
-
-    z.normalize(z);
-
-    for (int i = 0; i < rows; i++) {
-      v[0][i] = z[0][i];
-    }
-
-    /* Check for convergence */
-    if (std::fabs(d - d_old) < convergence_value) {
-      eigenvalues[0][k] = d;
-      k++;
-      d = 0.0;
-      for (int i = 0; i < rows; i++) {
-        v[0][i] = rand() / RAND_MAX;
-      }
-    }
-  }
+  // return determinant;
+  return 0;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::ElementMultiply(
-    const Matrix<rows, columns> &other, Matrix<rows, columns> &result) const {
+Matrix<rows, columns> &
+Matrix<rows, columns>::ElementMultiply(const Matrix<rows, columns> &other,
+                                       Matrix<rows, columns> &result) const {
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
     for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
       result[row_idx][column_idx] =
           this->Get(row_idx, column_idx) * other.Get(row_idx, column_idx);
     }
   }
+
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::ElementDivide(const Matrix<rows, columns> &other,
-                                          Matrix<rows, columns> &result) const {
+Matrix<rows, columns> &
+Matrix<rows, columns>::ElementDivide(const Matrix<rows, columns> &other,
+                                     Matrix<rows, columns> &result) const {
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
     for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
       result[row_idx][column_idx] =
           this->Get(row_idx, column_idx) / other.Get(row_idx, column_idx);
     }
   }
+
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
 float Matrix<rows, columns>::Get(uint8_t row_index,
                                  uint8_t column_index) const {
+  if (row_index > rows - 1 || column_index > columns - 1) {
+    return 0; // TODO: We should throw something here instead of failing quietly
+  }
   return this->matrix[row_index][column_index];
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::GetRow(uint8_t row_index,
-                                   Matrix<1, columns> &row) const {
+Matrix<1, columns> &
+Matrix<rows, columns>::GetRow(uint8_t row_index,
+                              Matrix<1, columns> &row) const {
   row = Matrix<1, columns>(this->matrix[row_index]);
+
+  return row;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::GetColumn(uint8_t column_index,
-                                      Matrix<rows, 1> &column) const {
+Matrix<rows, 1> &
+Matrix<rows, columns>::GetColumn(uint8_t column_index,
+                                 Matrix<rows, 1> &column) const {
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
     column[row_idx][0] = this->Get(row_idx, column_index);
   }
+
+  return column;
 }
 
 template <uint8_t rows, uint8_t columns>
@@ -471,6 +498,18 @@ float Matrix<rows, columns>::dotProduct(const Matrix<1, vector_size> &vec1,
   return sum;
 }
 
+template <uint8_t rows, uint8_t columns>
+template <uint8_t vector_size>
+float Matrix<rows, columns>::dotProduct(const Matrix<vector_size, 1> &vec1,
+                                        const Matrix<vector_size, 1> &vec2) {
+  float sum{0};
+  for (uint8_t i{0}; i < vector_size; i++) {
+    sum += vec1.Get(i, 0) * vec2.Get(i, 0);
+  }
+
+  return sum;
+}
+
 template <uint8_t rows, uint8_t columns>
 void Matrix<rows, columns>::zeroMatrix() {
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
@@ -481,8 +520,8 @@ void Matrix<rows, columns>::zeroMatrix() {
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::matrixOfMinors(
-    Matrix<rows, columns> &result) const {
+Matrix<rows, columns> &
+Matrix<rows, columns>::matrixOfMinors(Matrix<rows, columns> &result) const {
   Matrix<rows - 1, columns - 1> minorMatrix{};
 
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
@@ -491,12 +530,14 @@ void Matrix<rows, columns>::matrixOfMinors(
       result[row_idx][column_idx] = minorMatrix.Det();
     }
   }
+
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::minorMatrix(Matrix<rows - 1, columns - 1> &result,
-                                        uint8_t row_idx,
-                                        uint8_t column_idx) const {
+Matrix<rows - 1, columns - 1> &
+Matrix<rows, columns>::minorMatrix(Matrix<rows - 1, columns - 1> &result,
+                                   uint8_t row_idx, uint8_t column_idx) const {
   std::array<float, (rows - 1) * (columns - 1)> subArray{};
 
   for (uint8_t row_iter{0}; row_iter < rows; row_iter++) {
@@ -511,10 +552,12 @@ void Matrix<rows, columns>::minorMatrix(Matrix<rows - 1, columns - 1> &result,
   }
 
   result = Matrix<rows - 1, columns - 1>{subArray};
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::adjugate(Matrix<rows, columns> &result) const {
+Matrix<rows, columns> &
+Matrix<rows, columns>::adjugate(Matrix<rows, columns> &result) const {
   for (uint8_t row_iter{0}; row_iter < rows; row_iter++) {
     for (uint8_t column_iter{0}; column_iter < columns; column_iter++) {
       float sign = ((row_iter + 1) % 2) == 0 ? -1 : 1;
@@ -522,26 +565,34 @@ void Matrix<rows, columns>::adjugate(Matrix<rows, columns> &result) const {
       result[row_iter][column_iter] = this->Get(row_iter, column_iter) * sign;
     }
   }
+
+  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::normalize(Matrix<rows, columns> &result) const {
+Matrix<rows, columns> &
+Matrix<rows, columns>::Normalize(Matrix<rows, columns> &result) const {
   float sum{0};
-  for (uint8_t column_idx{0}; column_idx < rows; column_idx++) {
-    for (uint8_t row_idx{0}; row_idx < columns; row_idx++) {
-      sum += this->Get(row_idx, column_idx);
+  for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
+    for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
+      float val{this->Get(row_idx, column_idx)};
+      sum += val * val;
     }
   }
 
   if (sum == 0) {
     // this wouldn't do anything anyways
     result.zeroMatrix();
-    return;
+    return result;
   }
 
-  for (uint8_t column_idx{0}; column_idx < rows; column_idx++) {
-    for (uint8_t row_idx{0}; row_idx < columns; row_idx++) {
+  sum = sqrt(sum);
+
+  for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
+    for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
       result[row_idx][column_idx] = this->Get(row_idx, column_idx) / sum;
     }
   }
+
+  return result;
 }
\ No newline at end of file
diff --git a/unit-tests/matrix-tests.cpp b/unit-tests/matrix-tests.cpp
index 25d2d0f..474afee 100644
--- a/unit-tests/matrix-tests.cpp
+++ b/unit-tests/matrix-tests.cpp
@@ -1,14 +1,16 @@
 // include the unit test framework first
 #include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
 
 // include the module you're going to test next
 #include "Matrix.hpp"
 
 // any other libraries
 #include <array>
+#include <cmath>
 #include <iostream>
 
-TEST_CASE("Elementary Matrix Operations", "Matrix::Add") {
+TEST_CASE("Elementary Matrix Operations", "Matrix") {
   std::array<float, 4> arr1{1, 2, 3, 4};
   std::array<float, 4> arr2{5, 6, 7, 8};
   Matrix<2, 2> mat1{arr1};
@@ -58,6 +60,8 @@ TEST_CASE("Elementary Matrix Operations", "Matrix::Add") {
     REQUIRE(mat3.Get(0, 1) == 22);
     REQUIRE(mat3.Get(1, 0) == 43);
     REQUIRE(mat3.Get(1, 1) == 50);
+
+    // try a non-square matrix
   }
 
   SECTION("Scalar Multiplication") {
@@ -68,4 +72,101 @@ TEST_CASE("Elementary Matrix Operations", "Matrix::Add") {
     REQUIRE(mat3.Get(1, 0) == 6);
     REQUIRE(mat3.Get(1, 1) == 8);
   }
+
+  SECTION("Squaring") {
+    mat1.Square(mat3);
+
+    REQUIRE(mat3.Get(0, 0) == 7);
+    REQUIRE(mat3.Get(0, 1) == 10);
+    REQUIRE(mat3.Get(1, 0) == 15);
+    REQUIRE(mat3.Get(1, 1) == 22);
+  }
+
+  SECTION("Element Multiply") {
+    mat1.ElementMultiply(mat2, mat3);
+
+    REQUIRE(mat3.Get(0, 0) == 5);
+    REQUIRE(mat3.Get(0, 1) == 12);
+    REQUIRE(mat3.Get(1, 0) == 21);
+    REQUIRE(mat3.Get(1, 1) == 32);
+  }
+
+  SECTION("Element Divide") {
+    mat1.ElementDivide(mat2, mat3);
+
+    REQUIRE(mat3.Get(0, 0) == 1 / 5);
+    REQUIRE(mat3.Get(0, 1) == 2 / 6);
+    REQUIRE(mat3.Get(1, 0) == 3 / 7);
+    REQUIRE(mat3.Get(1, 1) == 4 / 8);
+  }
+
+  SECTION("Determinant") {
+    float det1 = mat1.Det();
+
+    REQUIRE_THAT(det1, Catch::Matchers::WithinRel(-2.0F, 1e-6f));
+
+    std::array<float, 9> arr4{1, 2, 3, 4, 5, 6, 7, 8, 9};
+    Matrix<3, 3> mat4{arr4};
+
+    float det4 = mat4.Det();
+
+    REQUIRE_THAT(det4, Catch::Matchers::WithinRel(0.0F, 1e-6f));
+
+    std::array<float, 9> arr5{1, 0, 0, 0, 2, 0, 0, 0, 3};
+    Matrix<3, 3> mat5{arr5};
+
+    float det5 = mat5.Det();
+    REQUIRE_THAT(det5, Catch::Matchers::WithinRel(6.0F, 1e-6f));
+  }
+
+  SECTION("Invert"){};
+
+  SECTION("Transpose") {
+    // transpose a square matrix
+    mat1.Transpose(mat3);
+
+    REQUIRE(mat3.Get(0, 0) == 1);
+    REQUIRE(mat3.Get(0, 1) == 3);
+    REQUIRE(mat3.Get(1, 0) == 2);
+    REQUIRE(mat3.Get(1, 1) == 4);
+
+    // transpose a non-square matrix
+    std::array<float, 6> arr4{1, 2, 3, 4, 5, 6};
+    Matrix<2, 3> mat4{arr4};
+    Matrix<3, 2> mat5{};
+
+    mat4.Transpose(mat5);
+
+    REQUIRE(mat5.Get(0, 0) == 1);
+    REQUIRE(mat5.Get(0, 1) == 4);
+    REQUIRE(mat5.Get(1, 0) == 2);
+    REQUIRE(mat5.Get(1, 1) == 5);
+    REQUIRE(mat5.Get(2, 0) == 3);
+    REQUIRE(mat5.Get(2, 1) == 6);
+  }
+
+  SECTION("Normalize") {
+    mat1.Normalize(mat3);
+
+    float sqrt_30{sqrt(30)};
+
+    REQUIRE(mat3.Get(0, 0) == 1 / sqrt_30);
+    REQUIRE(mat3.Get(0, 1) == 2 / sqrt_30);
+    REQUIRE(mat3.Get(1, 0) == 3 / sqrt_30);
+    REQUIRE(mat3.Get(1, 1) == 4 / sqrt_30);
+
+    std::array<float, 2> arr4{-0.878877044, 2.92092276};
+    Matrix<2, 1> mat4{arr4};
+    Matrix<2, 1> mat5{};
+    mat4.Normalize(mat5);
+
+    REQUIRE_THAT(mat5.Get(0, 0),
+                 Catch::Matchers::WithinRel(-0.288129855179f, 1e-6f));
+    REQUIRE_THAT(mat5.Get(1, 0),
+                 Catch::Matchers::WithinRel(0.957591346325f, 1e-6f));
+  }
+
+  SECTION("GET ROW") {}
+
+  SECTION("GET COLUMN") {}
 }
\ No newline at end of file