]> https://gitweb.dealii.org/ - dealii.git/commitdiff
factor out hyperbolic rotations into Utilities::LinearAlgebra::hyperbolic_rotation...
authorDenis Davydov <davydden@gmail.com>
Mon, 18 Dec 2017 16:53:50 +0000 (17:53 +0100)
committerDenis Davydov <davydden@gmail.com>
Tue, 19 Dec 2017 12:27:14 +0000 (13:27 +0100)
include/deal.II/lac/utilities.h
source/lac/lapack_full_matrix.cc
tests/lac/utilities_05.cc [new file with mode: 0644]
tests/lac/utilities_05.output [new file with mode: 0644]

index 9c0bf3da699ccbb595ce834406eb57aec3bd30da..b1c56f5bf4929173a89846fccd84fcec20651431 100644 (file)
@@ -58,11 +58,45 @@ namespace Utilities
      * \f]
      *
      * @note The function is implemented for real valued numbers only.
+     *
+     * @author Denis Davydov, 2017
      */
     template<typename NumberType>
     std::array<NumberType,3> 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<typename NumberType>
+    std::array<NumberType,3> 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<typename NumberType>
+    std::array<NumberType,3> 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<NumberType,3> csr;
+      csr[0] = 1./u;          // c
+      csr[1] = csr[0] * tau;  // s
+      csr[2] = f *u;          // r
+      return csr;
+    }
+
+    template<typename NumberType>
+    std::array<std::complex<NumberType>,3> hyperbolic_rotation(const std::complex<NumberType> &f,
+                                                               const std::complex<NumberType> &g)
+    {
+      AssertThrow(false, ExcNotImplemented());
+      std::array<NumberType,3> res;
+      return res;
+    }
+
+
     template<typename NumberType>
     std::array<std::complex<NumberType>,3> Givens_rotation(const std::complex<NumberType> &f,
                                                            const std::complex<NumberType> &g)
     {
       AssertThrow(false, ExcNotImplemented());
+      std::array<NumberType,3> res;
+      return res;
     }
 
     template<typename NumberType>
index 77d18620bad6cf2797141be8a80c229d4277d260..db569d11a6e978cdf424e4efd9a35b548c881122 100644 (file)
@@ -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<number,3> 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 (file)
index 0000000..9283949
--- /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 hyperbolic 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::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<double>( 1.,  0.);  // g == 0
+  test<double>( 2.,  1.);  // both positive
+  test<double>( 3., -0.5); // g negative
+  test<double>(-4., -2.4); // both negative
+  test<double>(-5.,  2);   // f negative
+}
diff --git a/tests/lac/utilities_05.output b/tests/lac/utilities_05.output
new file mode 100644 (file)
index 0000000..1f17d1a
--- /dev/null
@@ -0,0 +1,6 @@
+
+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.