eigenvalues,
eigenvectors,
eigenvalues.size());
+
+ {
+ const double precision = 1e-7;
+ Vector<double> Ax(eigenvectors[0]), Bx(eigenvectors[0]);
+ for (unsigned int i=0; i < eigenvectors.size(); ++i)
+ {
+ B.vmult(Bx,eigenvectors[i]);
+
+ for (unsigned int j=0; j < eigenvectors.size(); j++)
+ Assert( std::abs( eigenvectors[j] * Bx - (i==j))< precision,
+ ExcMessage("Eigenvectors " +
+ Utilities::int_to_string(i) +
+ " and " +
+ Utilities::int_to_string(j) +
+ " are not orthonormal!"));
+
+ A.vmult(Ax,eigenvectors[i]);
+ Ax.add(-1.0*std::real(eigenvalues[i]),Bx);
+ Assert (Ax.l2_norm() < precision,
+ ExcMessage("Returned vector " +
+ Utilities::int_to_string(i) +
+ " is not an eigenvector!"));
+ }
+ }
}
PETScWrappers::PreconditionBlockJacobi preconditioner (A);\r
cg_solver.solve (A, x, b, preconditioner);\r
\r
- // some output\r
- // ?\r
+ PETScWrappers::Vector res(x);\r
+ A.residual(res,x,b);\r
+ AssertThrow(res.l2_norm()<1e-3,\r
+ ExcInternalError());\r
}\r
\r
void LaplaceProblem::run ()\r
\r
void LaplaceEigenspectrumProblem::solve ()\r
{\r
- SolverControl solver_control (1000, 1e-03);\r
+ SolverControl solver_control (1000, 1e-10);\r
SLEPcWrappers::SolverArnoldi eigensolver (solver_control);\r
eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL);\r
eigensolver.solve (A, B, lambda, x, x.size());\r
\r
+ {\r
+ const double precision = 1e-7;\r
+ PETScWrappers::Vector Ax(x[0]), Bx(x[0]);\r
+ for (unsigned int i=0; i < x.size(); ++i)\r
+ {\r
+ B.vmult(Bx,x[i]);\r
+\r
+ for (unsigned int j=0; j < x.size(); j++)\r
+ if (j!=i)\r
+ Assert( std::abs( x[j] * Bx )< precision,\r
+ ExcMessage("Eigenvectors " +\r
+ Utilities::int_to_string(i) +\r
+ " and " +\r
+ Utilities::int_to_string(j) +\r
+ " are not orthogonal!"));\r
+\r
+ A.vmult(Ax,x[i]);\r
+ Ax.add(-1.0*lambda[i],Bx);\r
+ Assert (Ax.l2_norm() < precision,\r
+ ExcMessage("Returned vector " +\r
+ Utilities::int_to_string(i) +\r
+ " is not an eigenvector!"));\r
+ }\r
+ }\r
+\r
+\r
// some output\r
output_table.add_value ("lambda", lambda[0]);\r
output_table.add_value ("error", std::fabs(2.-lambda[0]));\r
preconditioner.initialize (A);\r
cg_solver.solve (A, x, b, preconditioner);\r
\r
- // some output\r
- // ?\r
+ TrilinosWrappers::Vector res(x);\r
+ A.residual(res,x,b);\r
+ AssertThrow(res.l2_norm()<1e-3,\r
+ ExcInternalError());\r
}\r
\r
void LaplaceProblem::run ()\r