Frobenius norm for matrix

src/algebra/TinyMatrix.hpp

@@ -504,6 +504,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)
+{
+  static_assert(std::is_arithmetic<T>::value, "norm is not defined for non-arithmetic types");
+
+  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);
+}
+
 template <size_t N, typename T>
 PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
 

tests/test_TinyMatrix.cpp

@@ -209,6 +209,16 @@ TEST_CASE("TinyMatrix", "[algebra]")
     REQUIRE(trace(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 1 + 0 + 1 + 2);
   }
 
+  SECTION("checking for norm calculations")
+  {
+    REQUIRE(frobeniusNorm(TinyMatrix<1, 1, int>(6)) == 6);
+    REQUIRE(frobeniusNorm(TinyMatrix<2, 2>(5, 1, -3, 6)) == Catch::Approx(std::sqrt(25 + 1 + 9 + 36)).epsilon(1E-14));
+    REQUIRE(frobeniusNorm(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
+            Catch::Approx(std::sqrt(1 + 2.3 * 2.3 + 7 * 7 + 6.2 * 6.2 + 9 + 16 + 81 + 1 + 4.1 * 4.1 + 25 + 4 + 9 + 4 +
+                                    27 * 27 + 9 + 17.5 * 17.5))
+              .epsilon(1E-14));
+  }
+
   SECTION("checking for inverse calculations")
   {
     {