From f9b783d3b3c4468f3bb6ab74ad56c67a67f82b34 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 8 Apr 2018 16:37:50 +0200 Subject: [PATCH] Add another test. --- tests/base/polynomial_jacobi_roots.cc | 67 + ...omial_jacobi_roots.with_lapack=true.output | 1145 +++++++++++++++++ tests/base/polynomial_jacobi_roots_ld.cc | 67 + ...al_jacobi_roots_ld.with_lapack=true.output | 1145 +++++++++++++++++ 4 files changed, 2424 insertions(+) create mode 100644 tests/base/polynomial_jacobi_roots.cc create mode 100644 tests/base/polynomial_jacobi_roots.with_lapack=true.output create mode 100644 tests/base/polynomial_jacobi_roots_ld.cc create mode 100644 tests/base/polynomial_jacobi_roots_ld.with_lapack=true.output 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