const unsigned int n_eigenvecs = eigenvectors.size();
for (size_type i=0; i<n_eigenvecs; ++i)
for (unsigned int j=0; j<n; ++j)
- eigenvectors[i](j) = z[i*n+j];
+ eigenvectors[i](j) = v[i*n+j];
delete[] workd;
eigenfunctions,
eigenvalues.size());
+ // 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}
+ {
+ Vector<double> 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))< 1e-8,
+ ExcMessage(std::to_string(eigenfunctions[j] * Bx)));
+
+ stiffness_matrix.vmult(Ax,eigenfunctions[i]);
+ Ax.add(-1.0*std::real(eigenvalues[i]),Bx);
+ Assert (Ax.l2_norm() < 1e-8,
+ ExcMessage(std::to_string(Ax.l2_norm())));
+ }
+ }
for (unsigned int i=0; i<eigenfunctions.size(); ++i)
eigenfunctions[i] /= eigenfunctions[i].linfty_norm ();