PArpack: extend unit tests to check that we get correct eigenvectors

diff --git a/tests/arpack/step-36_parpack.cc b/tests/arpack/step-36_parpack.cc
index 98387375ebcc1a2033c9b9ca3e6db917d7a79e6e..cb8ec9ab5970083f3a945a21771f542cc670a98b 100644 (file)
--- a/tests/arpack/step-36_parpack.cc
+++ b/tests/arpack/step-36_parpack.cc
@@ -296,6 +296,26 @@ void test ()
     for (unsigned int i=0; i < eigenvalues.size(); i++)
       dealii::deallog << eigenvalues[i] << std::endl;
 
+    // make sure that we have eigenvectors and they are mass-orthonormal:
+    // a) (A*x_i-\lambda*B*x_i).L2() == 0
+    // b) x_i*B*y_i=\delta_{ij}
+    {
+      const double precision = 1e-7;
+      PETScWrappers::MPI::Vector Ax(eigenfunctions[0]), Bx(eigenfunctions[0]);
+      for (unsigned int i=0; i < eigenfunctions.size(); ++i)
+        {
+          mass_matrix.vmult(Bx,eigenfunctions[i]);
+
+          for (unsigned int j=0; j < eigenfunctions.size(); j++)
+            Assert( std::abs( eigenfunctions[j] * Bx - (i==j))< precision,
+                    ExcMessage(std::to_string(eigenfunctions[j] * Bx)));
+
+          stiffness_matrix.vmult(Ax,eigenfunctions[i]);
+          Ax.add(-1.0*std::real(lambda[i]),Bx);
+          Assert (Ax.l2_norm() < precision,
+                  ExcMessage(std::to_string(Ax.l2_norm())));
+        }
+    }
   }
 
 

diff --git a/tests/arpack/step-36_parpack_trilinos.cc b/tests/arpack/step-36_parpack_trilinos.cc
index 17d18057f266877cad2cc3884821ef02a4e66efe..31396a8b28471a25dfb15b24df722b67a5e74227 100644 (file)
--- a/tests/arpack/step-36_parpack_trilinos.cc
+++ b/tests/arpack/step-36_parpack_trilinos.cc
@@ -290,6 +290,26 @@ void test ()
     for (unsigned int i=0; i < eigenvalues.size(); i++)
       deallog << eigenvalues[i] << std::endl;
 
+    // make sure that we have eigenvectors and they are mass-orthonormal:
+    // a) (A*x_i-\lambda*B*x_i).L2() == 0
+    // b) x_i*B*y_i=\delta_{ij}
+    {
+      const double precision = 1e-7;
+      TrilinosWrappers::MPI::Vector Ax(eigenfunctions[0]), Bx(eigenfunctions[0]);
+      for (unsigned int i=0; i < eigenfunctions.size(); ++i)
+        {
+          mass_matrix.vmult(Bx,eigenfunctions[i]);
+
+          for (unsigned int j=0; j < eigenfunctions.size(); j++)
+            Assert( std::abs( eigenfunctions[j] * Bx - (i==j))< precision,
+                    ExcMessage(std::to_string(eigenfunctions[j] * Bx)));
+
+          stiffness_matrix.vmult(Ax,eigenfunctions[i]);
+          Ax.add(-1.0*std::real(lambda[i]),Bx);
+          Assert (Ax.l2_norm() < precision,
+                  ExcMessage(std::to_string(Ax.l2_norm())));
+        }
+    }
   }