Correct under- and overflow protection in l2_norm (#184)

Author: julianlitz

include/LinearAlgebra/Vector/vector_operations.h

@@ -112,18 +112,7 @@ T l1_norm(ConstVector<T> x)
     return result;
 }
 
-template <typename T>
-T l2_norm_squared(ConstVector<T> x)
-{
-    T result      = 0.0;
-    std::size_t n = x.size();
-#pragma omp parallel for reduction(+ : result) if (n > 10'000)
-    for (std::size_t i = 0; i < n; ++i) {
-        result += x(i) * x(i);
-    }
-    return result;
-}
-
+// Underflow- and overflow-resistant implementation of the L2 norm
 template <typename T>
 T l2_norm(ConstVector<T> x)
 {
@@ -137,17 +126,17 @@ T l2_norm(ConstVector<T> x)
             scale = abs_val;
         }
     }
-    if (equals(scale, T{0})) {
+    // 2) if the largest absolute value is zero, the norm is zero
+    if (scale == T{0})
         return T{0};
-    }
-    // 2) accumulate sum of squares of scaled entries
+    // 3) accumulate sum of squares of scaled entries
     T sum = 0.0;
 #pragma omp parallel for reduction(+ : sum) if (n > 10'000)
     for (std::size_t i = 0; i < n; ++i) {
         T value = x(i) / scale;
         sum += value * value;
     }
-    // 3) rescale
+    // 4) rescale
     return scale * std::sqrt(sum);
 }
 

tests/LinearAlgebra/Vector/vector_operations.cpp

@@ -163,28 +163,46 @@ TEST(VectorOperations, l1_vector_norm)
     EXPECT_DOUBLE_EQ(l1_norm(ConstVector<double>(v)), 8.0);
 }
 
-/* T l2_norm_squared(ConstVector<T>& x); */
+/* T l2_norm(ConstVector<T>& x); */
 
-TEST(VectorOperations, l2_vector_norm_squared)
+TEST(VectorOperations, l2_vector_norm)
 {
     Vector<double> v("v", 3);
     v(0) = 1;
     v(1) = -5;
     v(2) = 2;
     ConstVector<double> const_v(v);
-    EXPECT_DOUBLE_EQ(l2_norm_squared(const_v), 30.0);
+    EXPECT_DOUBLE_EQ(l2_norm(ConstVector<double>(const_v)), std::sqrt(30.0));
 }
 
-/* T l2_norm(ConstVector<T>& x); */
+TEST(VectorOperations, zero_l2_vector_norm)
+{
+    Vector<double> v("v", 3);
+    v(0) = 0;
+    v(1) = 0;
+    v(2) = 0;
+    ConstVector<double> const_v(v);
+    EXPECT_DOUBLE_EQ(l2_norm(ConstVector<double>(const_v)), 0.0);
+}
 
-TEST(VectorOperations, l2_vector_norm)
+TEST(VectorOperations, underflow_l2_vector_norm)
 {
     Vector<double> v("v", 3);
-    v(0) = 1;
-    v(1) = -5;
-    v(2) = 2;
+    v(0) = 1e-300;
+    v(1) = 1e-300;
+    v(2) = 1e-300;
     ConstVector<double> const_v(v);
-    EXPECT_DOUBLE_EQ(l2_norm(ConstVector<double>(const_v)), std::sqrt(30.0));
+    EXPECT_DOUBLE_EQ(l2_norm(ConstVector<double>(const_v)), std::sqrt(3.0) * 1e-300);
+}
+
+TEST(VectorOperations, overflow_l2_vector_norm)
+{
+    Vector<double> v("v", 3);
+    v(0) = 1e+300;
+    v(1) = 1e+300;
+    v(2) = 1e+300;
+    ConstVector<double> const_v(v);
+    EXPECT_DOUBLE_EQ(l2_norm(ConstVector<double>(const_v)), std::sqrt(3.0) * 1e+300);
 }
 
 /* T infinity_norm(ConstVector<T>& x); */