Unitarity check

matrix.hpp

@@ -123,7 +123,7 @@ public:
     }
     return m;
   }
-  const bool is_diagonal() {
+  const bool is_diagonal() const {
     bool result = true;
     if (this->num_entries != this->entry_dimension)
       return !result;
@@ -135,8 +135,12 @@ public:
     }
     return result;
   }
-  // Fixme: This is very dumb
+  const bool is_unitary() const{
-  const vector get_eigenvalues() {
+	  matrix m = this->hermitian_conjugate();
+	  return (*this * m == this->I());
+  }
+  // FIXME: This is very dumb
+  const vector get_eigenvalues() const {
     if (this->is_diagonal())
       return this->get_diagonal();
     return vector(0);