Make eigenvector computation of diagonal 2x2 matrices more stable

The eigenvector code used to fail for certain seemingly simple
diagonal 2x2 matrices.  The reason was that the code tried
to avoid certain all-zero matrix columns by testing whether
the norm is (close to) zero.  That proved to be surprisingly
unstable.  The fix here is simple: Pick the column with the
larger column norm (we know they are not both zero).

Also add a test case to make sure the bug is gone.

dune/common/fmatrixev.hh

@@ -326,26 +326,17 @@ namespace Dune {
             eigenVectors[1] = {0.0, 1.0};
           }
           else {
-            using std::isnormal;
-            // Try whether the first colum of A - λ_2I is a good eigenvector for λ_1
-            FieldVector<K,2> ev = {matrix[0][0]-eigenValues[1], matrix[1][0]};
-            K norm = std::sqrt(ev[0]*ev[0] + ev[1]*ev[1]);
-            if (!isnormal(norm)) // No, that column is zero, take the other one
-            {
-              ev = {matrix[0][1], matrix[1][1]-eigenValues[1]};
-              norm = std::sqrt(ev[0]*ev[0] + ev[1]*ev[1]);
-            }
-            eigenVectors[0] = ev/norm;
-
-            // Try whether the first colum of A - λ_1I is a good eigenvector for λ_2
-            ev = {matrix[0][0]-eigenValues[0], matrix[1][0]};
-            norm = std::sqrt(ev[0]*ev[0] + ev[1]*ev[1]);
-            if (!isnormal(norm)) // No, that column is zero, take the other one
-            {
-              ev = {matrix[0][1], matrix[1][1]-eigenValues[0]};
-              norm = std::sqrt(ev[0]*ev[0] + ev[1]*ev[1]);
-            }
-            eigenVectors[1] = ev/norm;
+            // The columns of A - λ_2I are eigenvectors for λ_1, or zero.
+            // Take the column with the larger norm to avoid zero columns.
+            FieldVector<K,2> ev0 = {matrix[0][0]-eigenValues[1], matrix[1][0]};
+            FieldVector<K,2> ev1 = {matrix[0][1], matrix[1][1]-eigenValues[1]};
+            eigenVectors[0] = (ev0.two_norm2() >= ev1.two_norm2()) ? ev0/ev0.two_norm() : ev1/ev1.two_norm();
+
+            // The columns of A - λ_1I are eigenvectors for λ_2, or zero.
+            // Take the column with the larger norm to avoid zero columns.
+            ev0 = {matrix[0][0]-eigenValues[0], matrix[1][0]};
+            ev1 = {matrix[0][1], matrix[1][1]-eigenValues[0]};
+            eigenVectors[1] = (ev0.two_norm2() >= ev1.two_norm2()) ? ev0/ev0.two_norm() : ev1/ev1.two_norm();
           }
         }
       }

dune/common/test/eigenvaluestest.cc

@@ -251,6 +251,9 @@ void checkMultiplicity()
   // another singular matrix (triggers a different code path)
   checkMatrixWithReference<FT,2>({{0, 0},{0, 1}}, {{1,0}, {0,1}}, {0, 1});
 
+  // Seemingly simple diagonal matrix -- triggers unstable detection of zero columns
+  checkMatrixWithReference<FT,2>({{1.01, 0},{0, 1}}, {{0,1}, {1,0}}, {1, 1.01});
+
   //--3d--
   //repeated eigenvalue (x3)
   checkMatrixWithReference<FT,3>({{  1,   0,   0},