From: Denis Davydov <davydden@gmail.com>
Date: Thu, 17 Sep 2015 20:25:37 +0000 (+0200)
Subject: PArpack: extend unit tests to check that we get correct eigenvectors
X-Git-Tag: v8.4.0-rc2~396^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d21c2dd440e9672771bcf8ec9d79b064e6d681ce;p=dealii.git

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 98387375eb..cb8ec9ab59 100644
--- 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 17d18057f2..31396a8b28 100644
--- 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())));
+        }
+    }
   }