Add TinyMatrix's double-dot product

Also add a bunch of non discard qualifiers for TinyMatrix and TinyVector and forbids Frobenius norm for non floating point TinyMatrix

src/algebra/TinyMatrix.hpp

@@ -56,15 +56,13 @@ class [[nodiscard]] TinyMatrix
   }
 
  public:
-  PUGS_INLINE
-  constexpr bool
+  [[nodiscard]] PUGS_INLINE constexpr bool
   isSquare() const noexcept
   {
     return M == N;
   }
 
-  PUGS_INLINE
-  constexpr friend TinyMatrix<N, M, T>
+  [[nodiscard]] PUGS_INLINE constexpr friend TinyMatrix<N, M, T>
   transpose(const TinyMatrix& A)
   {
     TinyMatrix<N, M, T> tA;
@@ -76,36 +74,31 @@ class [[nodiscard]] TinyMatrix
     return tA;
   }
 
-  PUGS_INLINE
-  constexpr size_t
+  [[nodiscard]] PUGS_INLINE constexpr size_t
   dimension() const
   {
     return M * N;
   }
 
-  PUGS_INLINE
-  constexpr size_t
+  [[nodiscard]] PUGS_INLINE constexpr size_t
   numberOfValues() const
   {
     return this->dimension();
   }
 
-  PUGS_INLINE
-  constexpr size_t
+  [[nodiscard]] PUGS_INLINE constexpr size_t
   numberOfRows() const
   {
     return M;
   }
 
-  PUGS_INLINE
-  constexpr size_t
+  [[nodiscard]] PUGS_INLINE constexpr size_t
   numberOfColumns() const
   {
     return N;
   }
 
-  PUGS_INLINE
-  constexpr TinyMatrix
+  PUGS_INLINE constexpr TinyMatrix
   operator-() const
   {
     TinyMatrix opposite;
@@ -140,6 +133,16 @@ class [[nodiscard]] TinyMatrix
     return *this;
   }
 
+  [[nodiscard]] PUGS_INLINE constexpr friend T
+  doubleDot(const TinyMatrix& A, const TinyMatrix& B)
+  {
+    T t = A.m_values[0] * B.m_values[0];
+    for (size_t i = 1; i < M * N; ++i) {
+      t += A.m_values[i] * B.m_values[i];
+    }
+    return t;
+  }
+
   template <size_t P>
   PUGS_INLINE constexpr TinyMatrix<M, P, T>
   operator*(const TinyMatrix<N, P, T>& B) const
@@ -194,8 +197,7 @@ class [[nodiscard]] TinyMatrix
     return os;
   }
 
-  PUGS_INLINE
-  constexpr bool
+  [[nodiscard]] PUGS_INLINE constexpr bool
   operator==(const TinyMatrix& A) const
   {
     for (size_t i = 0; i < M * N; ++i) {
@@ -205,8 +207,7 @@ class [[nodiscard]] TinyMatrix
     return true;
   }
 
-  PUGS_INLINE
-  constexpr bool
+  [[nodiscard]] PUGS_INLINE constexpr bool
   operator!=(const TinyMatrix& A) const
   {
     return not this->operator==(A);
@@ -272,24 +273,21 @@ class [[nodiscard]] TinyMatrix
     return *this;
   }
 
-  PUGS_INLINE
-  constexpr T&
+  [[nodiscard]] PUGS_INLINE constexpr T&
   operator()(size_t i, size_t j) noexcept(NO_ASSERT)
   {
     Assert((i < M) and (j < N));
     return m_values[_index(i, j)];
   }
 
-  PUGS_INLINE
-  constexpr const T&
+  [[nodiscard]] PUGS_INLINE constexpr const T&
   operator()(size_t i, size_t j) const noexcept(NO_ASSERT)
   {
     Assert((i < M) and (j < N));
     return m_values[_index(i, j)];
   }
 
-  PUGS_INLINE
-  constexpr TinyMatrix&
+  PUGS_INLINE constexpr TinyMatrix&
   operator=(ZeroType) noexcept
   {
     static_assert(std::is_arithmetic<T>(), "Cannot assign 'zero' value for non-arithmetic types");
@@ -377,7 +375,7 @@ class [[nodiscard]] TinyMatrix
 };
 
 template <size_t M, size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<M, N, T>
+[[nodiscard]] PUGS_INLINE constexpr TinyMatrix<M, N, T>
 tensorProduct(const TinyVector<M, T>& x, const TinyVector<N, T>& y)
 {
   TinyMatrix<M, N, T> A;
@@ -390,7 +388,7 @@ tensorProduct(const TinyVector<M, T>& x, const TinyVector<N, T>& y)
 }
 
 template <size_t N, typename T>
-PUGS_INLINE constexpr T
+[[nodiscard]] PUGS_INLINE constexpr T
 det(const TinyMatrix<N, N, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "determinant is not defined for non-arithmetic types");
@@ -441,7 +439,7 @@ det(const TinyMatrix<N, N, T>& A)
 }
 
 template <typename T>
-PUGS_INLINE constexpr T
+[[nodiscard]] PUGS_INLINE constexpr T
 det(const TinyMatrix<1, 1, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
@@ -449,7 +447,7 @@ det(const TinyMatrix<1, 1, T>& A)
 }
 
 template <typename T>
-PUGS_INLINE constexpr T
+[[nodiscard]] PUGS_INLINE constexpr T
 det(const TinyMatrix<2, 2, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
@@ -457,7 +455,7 @@ det(const TinyMatrix<2, 2, T>& A)
 }
 
 template <typename T>
-PUGS_INLINE constexpr T
+[[nodiscard]] PUGS_INLINE constexpr T
 det(const TinyMatrix<3, 3, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
@@ -466,7 +464,7 @@ det(const TinyMatrix<3, 3, T>& A)
 }
 
 template <size_t M, size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<M - 1, N - 1, T>
+[[nodiscard]] PUGS_INLINE constexpr TinyMatrix<M - 1, N - 1, T>
 getMinor(const TinyMatrix<M, N, T>& A, size_t I, size_t J)
 {
   static_assert(M >= 2 and N >= 2, "minor calculation requires at least 2x2 matrices");
@@ -492,7 +490,7 @@ getMinor(const TinyMatrix<M, N, T>& A, size_t I, size_t J)
 }
 
 template <size_t N, typename T>
-PUGS_INLINE T
+[[nodiscard]] PUGS_INLINE T
 trace(const TinyMatrix<N, N, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "trace is not defined for non-arithmetic types");
@@ -504,26 +502,21 @@ trace(const TinyMatrix<N, N, T>& A)
   return t;
 }
 
-template <size_t N, typename T>
-PUGS_INLINE T
-frobeniusNorm(const TinyMatrix<N, N, T>& A)
+template <size_t M, size_t N, typename T>
+[[nodiscard]] PUGS_INLINE constexpr decltype(std::sqrt(std::declval<T>()))
+frobeniusNorm(const TinyMatrix<M, N, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "norm is not defined for non-arithmetic types");
+  static_assert(std::is_floating_point<T>::value, "Frobenius norm is defined for floating point types only");
 
-  T t = 0;
-  for (size_t i = 0; i < N; ++i) {
-    for (size_t j = 0; j < N; ++j) {
-      t += A(i, j) * A(i, j);
-    }
-  }
-  return std::sqrt(t);
+  return std::sqrt(doubleDot(A, A));
 }
 
 template <size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
+[[nodiscard]] PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
 
 template <typename T>
-PUGS_INLINE constexpr TinyMatrix<1, 1, T>
+[[nodiscard]] PUGS_INLINE constexpr TinyMatrix<1, 1, T>
 inverse(const TinyMatrix<1, 1, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
@@ -534,7 +527,7 @@ inverse(const TinyMatrix<1, 1, T>& A)
 }
 
 template <size_t N, typename T>
-PUGS_INLINE constexpr T
+[[nodiscard]] PUGS_INLINE constexpr T
 cofactor(const TinyMatrix<N, N, T>& A, size_t i, size_t j)
 {
   static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non-arithmetic types");
@@ -544,7 +537,7 @@ cofactor(const TinyMatrix<N, N, T>& A, size_t i, size_t j)
 }
 
 template <typename T>
-PUGS_INLINE constexpr TinyMatrix<2, 2, T>
+[[nodiscard]] PUGS_INLINE constexpr TinyMatrix<2, 2, T>
 inverse(const TinyMatrix<2, 2, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
@@ -558,7 +551,7 @@ inverse(const TinyMatrix<2, 2, T>& A)
 }
 
 template <typename T>
-PUGS_INLINE constexpr TinyMatrix<3, 3, T>
+[[nodiscard]] PUGS_INLINE constexpr TinyMatrix<3, 3, T>
 inverse(const TinyMatrix<3, 3, T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");