#include <deal.II/lac/lapack_support.h>
#include <deal.II/lac/vector_memory.h>
+#include <array>
+#include <complex>
+
DEAL_II_NAMESPACE_OPEN
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<typename NumberType>
+ std::array<NumberType,3> 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
{
namespace LinearAlgebra
{
+
+ template<typename NumberType>
+ std::array<std::complex<NumberType>,3> Givens_rotation(const std::complex<NumberType> &f,
+ const std::complex<NumberType> &g)
+ {
+ AssertThrow(false, ExcNotImplemented());
+ }
+
+ template<typename NumberType>
+ std::array<NumberType,3> Givens_rotation(const NumberType &f,
+ const NumberType &g)
+ {
+ std::array<NumberType,3> 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 <typename OperatorType, typename VectorType>
double Lanczos_largest_eigenvalue(const OperatorType &H,
const VectorType &v0_,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/utilities.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+template <typename NumberType>
+void test (const NumberType a, const NumberType b)
+{
+
+ FullMatrix<NumberType> rotation(2);
+ Vector<NumberType> x(2), y(2), res(2);
+
+ x[0] = a;
+ x[1] = b;
+ y[1] = NumberType();
+
+ const std::array<NumberType,3> 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<double>(1., 0.);
+ test<double>(0., 1.);
+ test<double>(1., -2.);
+ test<double>(-1., 2.);
+ test<double>(-15.,3.);
+ test<double>(15.,-3.);
+}