* \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
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>
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);
}
}
-
}
}
--- /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 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
+}
--- /dev/null
+
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000