]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add Utilities::LinearAlgebra::Givens_rotation()
authorDenis Davydov <davydden@gmail.com>
Mon, 18 Dec 2017 12:36:35 +0000 (13:36 +0100)
committerDenis Davydov <davydden@gmail.com>
Tue, 19 Dec 2017 12:27:13 +0000 (13:27 +0100)
include/deal.II/lac/utilities.h
tests/lac/utilities_04.cc [new file with mode: 0644]
tests/lac/utilities_04.output [new file with mode: 0644]

index e1530f90c5d4254a24dd5a888e2d1345f0568b70..9c0bf3da699ccbb595ce834406eb57aec3bd30da 100644 (file)
@@ -23,6 +23,9 @@
 #include <deal.II/lac/lapack_support.h>
 #include <deal.II/lac/vector_memory.h>
 
+#include <array>
+#include <complex>
+
 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<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
@@ -144,6 +174,63 @@ namespace Utilities
 {
   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_,
diff --git a/tests/lac/utilities_04.cc b/tests/lac/utilities_04.cc
new file mode 100644 (file)
index 0000000..bdff2d0
--- /dev/null
@@ -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 <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.);
+}
diff --git a/tests/lac/utilities_04.output b/tests/lac/utilities_04.output
new file mode 100644 (file)
index 0000000..687ddbe
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.