From: Martin Kronbichler Date: Sun, 8 Apr 2018 14:37:50 +0000 (+0200) Subject: Add another test. X-Git-Tag: v9.0.0-rc1~195^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f9b783d3b3c4468f3bb6ab74ad56c67a67f82b34;p=dealii.git Add another test. --- diff --git a/tests/base/polynomial_jacobi_roots.cc b/tests/base/polynomial_jacobi_roots.cc new file mode 100644 index 0000000000..ab74663b3e --- /dev/null +++ b/tests/base/polynomial_jacobi_roots.cc @@ -0,0 +1,67 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// check jacobi_polynomial_roots by comparison to eigenvalue problem with +// LAPACK + +#include "../tests.h" + +#include +#include +#include + +using namespace Polynomials; + + +int main() +{ + initlog(); + + for (int alpha = 0; alpha<4; ++alpha) + for (unsigned int degree=1; degree<40; degree+=3) + { + const int beta = alpha; + deallog << "Jacobi_" << degree << "^(" << alpha << "," << beta << ")" << std::endl; + + std::vector roots = jacobi_polynomial_roots(degree, alpha, beta); + AssertDimension(roots.size(), degree); + Vector roots_reference(degree); + LAPACKFullMatrix eigenvalue_mat(degree, degree); + for (unsigned int k=1; k eigenvectors(roots_reference.size(),roots_reference.size()); + eigenvalue_mat.compute_eigenvalues_symmetric(-1., 1., 1.e-20, + roots_reference, + eigenvectors); + deallog << "Roots (implemented vs reference)" << std::endl; + for (unsigned int i=0; i +#include +#include + +using namespace Polynomials; + + +int main() +{ + initlog(); + + for (int alpha = 0; alpha<4; ++alpha) + for (unsigned int degree=1; degree<40; degree+=3) + { + const int beta = alpha; + deallog << "Jacobi_" << degree << "^(" << alpha << "," << beta << ")" << std::endl; + + std::vector roots = jacobi_polynomial_roots(degree, alpha, beta); + AssertDimension(roots.size(), degree); + Vector roots_reference(degree); + LAPACKFullMatrix eigenvalue_mat(degree, degree); + for (unsigned int k=1; k eigenvectors(roots_reference.size(),roots_reference.size()); + eigenvalue_mat.compute_eigenvalues_symmetric(-1., 1., 1.e-20, + roots_reference, + eigenvectors); + deallog << "Roots (implemented vs reference)" << std::endl; + for (unsigned int i=0; i