From 701b8fe61bec67a86ba4f0a3644294e27c8ce261 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 18 Dec 2017 17:53:50 +0100 Subject: [PATCH] factor out hyperbolic rotations into Utilities::LinearAlgebra::hyperbolic_rotation() and test --- include/deal.II/lac/utilities.h | 66 +++++++++++++++++++++++++++ source/lac/lapack_full_matrix.cc | 17 ++----- tests/lac/utilities_05.cc | 76 ++++++++++++++++++++++++++++++++ tests/lac/utilities_05.output | 6 +++ 4 files changed, 152 insertions(+), 13 deletions(-) create mode 100644 tests/lac/utilities_05.cc create mode 100644 tests/lac/utilities_05.output diff --git a/include/deal.II/lac/utilities.h b/include/deal.II/lac/utilities.h index 9c0bf3da69..b1c56f5bf4 100644 --- a/include/deal.II/lac/utilities.h +++ b/include/deal.II/lac/utilities.h @@ -58,11 +58,45 @@ namespace Utilities * \f] * * @note The function is implemented for real valued numbers only. + * + * @author Denis Davydov, 2017 */ template std::array Givens_rotation(const NumberType &x, const NumberType &y); + /** + * Return elements of hyperbolic rotation matrix. + * + * That is for a given + * pair @p x and @p y, return $c$ , $s$ and $r$ such that + * \f[ + * \begin{bmatrix} + * c & -s \\ + * -s & c + * \end{bmatrix} + * \begin{bmatrix} + * x \\ + * y + * \end{bmatrix} + * = + * \begin{bmatrix} + * r \\ + * 0 + * \end{bmatrix} + * \f] + * + * Real valued solution only exists if $|x|>|g|$, the function will + * throw an error otherwise. + * + * @note The function is implemented for real valued numbers only. + * + * @author Denis Davydov, 2017 + */ + template + std::array hyperbolic_rotation(const NumberType &x, + const NumberType &y); + /** * Estimate an upper bound for the largest eigenvalue of @p H by a @p k -step * Lanczos process starting from the initial vector @p v0. Typical @@ -175,11 +209,43 @@ namespace Utilities namespace LinearAlgebra { + template + std::array hyperbolic_rotation(const NumberType &f, + const NumberType &g) + { + Assert (f != 0, ExcDivideByZero()); + const NumberType tau = g/f; + AssertThrow (std::abs(tau) < 1., + ExcMessage("real-valued Hyperbolic rotation does not exist for ("+ + std::to_string(f) + + "," + + std::to_string(g)+ + ")")); + const NumberType u = std::copysign(sqrt((1.-tau)*(1.+tau)), f); // <-- more stable than std::sqrt(1.-tau*tau) + std::array csr; + csr[0] = 1./u; // c + csr[1] = csr[0] * tau; // s + csr[2] = f *u; // r + return csr; + } + + template + std::array,3> hyperbolic_rotation(const std::complex &f, + const std::complex &g) + { + AssertThrow(false, ExcNotImplemented()); + std::array res; + return res; + } + + template std::array,3> Givens_rotation(const std::complex &f, const std::complex &g) { AssertThrow(false, ExcNotImplemented()); + std::array res; + return res; } template diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 77d18620ba..db569d11a6 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -281,24 +281,15 @@ namespace z *= std::sqrt(-a); for (unsigned int k = 0; k < N; ++k) { - // we have Cholesky factor of SPD matrix - Assert (A(k,k) != 0, ExcInternalError()); - const number tau = z(k)/A(k,k); - AssertThrow (std::abs(tau) < 1., - ExcMessage("Cholesky downdating led to negative definite matrix")); - const number u = sqrt((1.-tau)*(1.+tau)); // <-- more stable than std::sqrt(1.-tau*tau) - const number c = 1./u; - const number s = c * tau; - const number r = u * std::abs(A(k,k)); - A(k,k) = r; + const std::array csr = Utilities::LinearAlgebra::hyperbolic_rotation(A(k,k),z(k)); + A(k,k) = csr[2]; for (unsigned int i = k+1; i < N; ++i) { const number t = A(i,k); - A(i,k) = c * A(i,k) - s * z(i); - z(i) = -s * t + c * z(i); + A(i,k) = csr[0] * A(i,k) - csr[1] * z(i); + z(i) = -csr[1] * t + csr[0] * z(i); } } - } } diff --git a/tests/lac/utilities_05.cc b/tests/lac/utilities_05.cc new file mode 100644 index 0000000000..92839492af --- /dev/null +++ b/tests/lac/utilities_05.cc @@ -0,0 +1,76 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + + +// Test hyperbolic rotations. + + +#include "../tests.h" +#include +#include +#include + +template +void test (const NumberType a, const NumberType b) +{ + + FullMatrix rotation(2); + Vector x(2), y(2), res(2); + + x[0] = a; + x[1] = b; + y[1] = NumberType(); + + const std::array csr = Utilities::LinearAlgebra::hyperbolic_rotation(a,b); + + rotation(0,0) = csr[0]; // c + rotation(1,1) = csr[0]; // c + rotation(0,1) = -csr[1]; // -s + rotation(1,0) = -csr[1]; // -s + y[0] = csr[2]; // r + + rotation.residual(res, x, y); + + const NumberType norm = res.l2_norm(); + deallog << norm << std::endl; + + if (norm > 1e-12) + { + deallog << "x:" << std::endl; + x.print(deallog.get_file_stream()); + deallog << "Hyperbolic:" << std::endl; + rotation.print(deallog.get_file_stream(), 10, 6); + deallog << "y:" << std::endl; + y.print(deallog.get_file_stream()); + deallog << "res:" << std::endl; + res.print(deallog.get_file_stream()); + AssertThrow(false, ExcInternalError()); + } + +} + +int main() +{ + std::ofstream logfile("output"); + deallog << std::setprecision(6); + deallog.attach(logfile); + + // check all combinations with real solutions: + test( 1., 0.); // g == 0 + test( 2., 1.); // both positive + test( 3., -0.5); // g negative + test(-4., -2.4); // both negative + test(-5., 2); // f negative +} diff --git a/tests/lac/utilities_05.output b/tests/lac/utilities_05.output new file mode 100644 index 0000000000..1f17d1a21b --- /dev/null +++ b/tests/lac/utilities_05.output @@ -0,0 +1,6 @@ + +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 -- 2.39.5