From 0e704ded2e874623cc65654706276bf8402ce016 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 18 Dec 2017 13:36:35 +0100 Subject: [PATCH] add Utilities::LinearAlgebra::Givens_rotation() --- include/deal.II/lac/utilities.h | 87 +++++++++++++++++++++++++++++++++ tests/lac/utilities_04.cc | 76 ++++++++++++++++++++++++++++ tests/lac/utilities_04.output | 7 +++ 3 files changed, 170 insertions(+) create mode 100644 tests/lac/utilities_04.cc create mode 100644 tests/lac/utilities_04.output diff --git a/include/deal.II/lac/utilities.h b/include/deal.II/lac/utilities.h index e1530f90c5..9c0bf3da69 100644 --- a/include/deal.II/lac/utilities.h +++ b/include/deal.II/lac/utilities.h @@ -23,6 +23,9 @@ #include #include +#include +#include + DEAL_II_NAMESPACE_OPEN namespace Utilities @@ -33,6 +36,33 @@ namespace Utilities namespace LinearAlgebra { + /** + * Return elements of continuous Givens rotation matrix and norm of the input vector. + * + * That is for a given + * pair @p x and @p y, return $c$ , $s$ and $\sqrt{x^2+y^2}$ such that + * \f[ + * \begin{bmatrix} + * c & s \\ + * -s & c + * \end{bmatrix} + * \begin{bmatrix} + * x \\ + * y + * \end{bmatrix} + * = + * \begin{bmatrix} + * \sqrt{x^2+y^2} \\ + * 0 + * \end{bmatrix} + * \f] + * + * @note The function is implemented for real valued numbers only. + */ + template + std::array Givens_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 @@ -144,6 +174,63 @@ namespace Utilities { namespace LinearAlgebra { + + template + std::array,3> Givens_rotation(const std::complex &f, + const std::complex &g) + { + AssertThrow(false, ExcNotImplemented()); + } + + template + std::array Givens_rotation(const NumberType &f, + const NumberType &g) + { + std::array res; + // naieve calculation for "r" may overflow or underflow: + // c = x / \sqrt{x^2+y^2} + // s = -y / \sqrt{x^2+y^2} + + // See Golub 2013, Matrix computations, Chapter 5.1.8 + // Algorithm 5.1.3 + // and + // Anderson (2000), Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem. LAPACK Working Note 150, University of Tennessee, UT-CS-00-454, December 4, 2000. + // Algorithm 4 + // We implement the latter below: + if (g == NumberType()) + { + res[0] = std::copysign(1., f); + res[1] = NumberType(); + res[2] = std::abs(f); + } + else if (f == NumberType()) + { + res[0] = NumberType(); + res[1] = std::copysign(1., g); + res[2] = std::abs(g); + } + else if (std::abs(f) > std::abs(g)) + { + const NumberType tau = g/f; + const NumberType u = std::copysign(std::sqrt(1.+tau*tau), f); + res[0] = 1./u; // c + res[1] = res[0] * tau; // s + res[2] = f * u; // r + } + else + { + const NumberType tau = f/g; + const NumberType u = std::copysign(std::sqrt(1.+tau*tau), g); + res[1] = 1./u; // s + res[0] = res[1] * tau; // c + res[2] = g * u; // r + } + + return res; + } + + + template double Lanczos_largest_eigenvalue(const OperatorType &H, const VectorType &v0_, diff --git a/tests/lac/utilities_04.cc b/tests/lac/utilities_04.cc new file mode 100644 index 0000000000..bdff2d07fa --- /dev/null +++ b/tests/lac/utilities_04.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 Givens 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::Givens_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 double norm = res.l2_norm(); + deallog << norm << std::endl; + + if (norm > 1e-12) + { + deallog << "x:" << std::endl; + x.print(deallog.get_file_stream()); + deallog << "Givens:" << 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); + + test(1., 0.); + test(0., 1.); + test(1., -2.); + test(-1., 2.); + test(-15.,3.); + test(15.,-3.); +} diff --git a/tests/lac/utilities_04.output b/tests/lac/utilities_04.output new file mode 100644 index 0000000000..687ddbebbc --- /dev/null +++ b/tests/lac/utilities_04.output @@ -0,0 +1,7 @@ + +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 -- 2.39.5