--- /dev/null
+New: Add LAPACKFullMatrix::apply_givens_rotation() and Vector::apply_givens_rotation() to
+apply a Givens rotation matrix.
+<br>
+(Denis Davydov, 2018/12/30)
void
rank1_update(const number a, const Vector<number> &v);
+ /**
+ * Apply <a href="https://en.wikipedia.org/wiki/Givens_rotation">Givens
+ * rotation</a>
+ * @p csr (a triplet of cosine, sine and radius, see
+ * Utilities::LinearAlgebra::givens_rotation() for the definition of the
+ * rotation matrix $\mathbf G$)
+ * to this matrix in the plane spanned by the @p i'th and @p k'th unit vectors.
+ * If @p left is <code>true</code>, the rotation is applied from left
+ * $\mathbf A \leftarrow \mathbf G \mathbf A$
+ * and only rows @p i and @p k are affected.
+ * Otherwise, transpose of the rotation matrix is applied from right
+ * $\mathbf A \leftarrow \mathbf A \mathbf G^T$
+ * and only columns @p i and @p k are affected.
+ */
+ void
+ apply_givens_rotation(const std::array<number, 3> &csr,
+ const size_type i,
+ const size_type k,
+ const bool left = true);
+
/**
* Assignment from different matrix classes, performing the usual conversion
* to the transposed format expected by LAPACK. This assignment operator
void
grow_or_shrink(const size_type N);
+ /**
+ * Apply <a href="https://en.wikipedia.org/wiki/Givens_rotation">Givens
+ * rotation</a>
+ * @p csr (a triplet of cosine, sine and radius, see
+ * Utilities::LinearAlgebra::givens_rotation())
+ * to the vector in the plane spanned by the @p i'th and @p k'th unit vectors.
+ */
+ void
+ apply_givens_rotation(const std::array<Number, 3> &csr,
+ const size_type i,
+ const size_type k);
+
/**
* Change the dimension to that of the vector @p V. The same applies as for
* the other @p reinit function.
+template <typename Number>
+void
+Vector<Number>::apply_givens_rotation(const std::array<Number, 3> &csr,
+ const size_type i,
+ const size_type k)
+{
+ auto & V = *this;
+ const Number t = V(i);
+ V(i) = csr[0] * V(i) + csr[1] * V(k);
+ V(k) = -csr[1] * t + csr[0] * V(k);
+}
+
+
+
template <typename Number>
Vector<Number>::Vector(Vector<Number> &&v) noexcept
: Subscriptor(std::move(v))
+template <typename number>
+void
+LAPACKFullMatrix<number>::apply_givens_rotation(
+ const std::array<number, 3> &csr,
+ const size_type i,
+ const size_type k,
+ const bool left)
+{
+ auto &A = *this;
+ // see Golub 2013 "Matrix computations", p241 5.1.9 Applying Givens
+ // Rotations but note the difference in notation, namely the sign of s: we
+ // have G * A, where G[1,1] = s
+ if (left)
+ {
+ for (size_type j = 0; j < A.n(); ++j)
+ {
+ const number t = A(i, j);
+ A(i, j) = csr[0] * A(i, j) + csr[1] * A(k, j);
+ A(k, j) = -csr[1] * t + csr[0] * A(k, j);
+ }
+ }
+ else
+ {
+ for (size_type j = 0; j < A.m(); ++j)
+ {
+ const number t = A(j, i);
+ A(j, i) = csr[0] * A(j, i) + csr[1] * A(j, k);
+ A(j, k) = -csr[1] * t + csr[0] * A(j, k);
+ }
+ }
+}
+
+
+
template <typename number>
void
LAPACKFullMatrix<number>::remove_row_and_column(const size_type row,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test application of Givens rotations to matrix and vector
+ */
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+/*
+ * MWE in Python:
+
+import numpy as np
+from scipy import linalg
+
+A = np.array([[1.1, 1, 2, 7], [3, 4, 5, 1], [6, 3, 4.5, 2.2], [7,9, 8, 2.3]])
+V = np.array([3,5,2,1])
+c = 0.949833
+s = 0.312758
+
+i = 1
+k = 3
+G = np.array([[1,0,0,0], [0,c,0,s], [0,0,1,0], [0,-s,0,c]])
+
+# from left
+A1 = G.dot(A)
+
+# from right
+A2 = A.dot(G.T)
+
+# vector
+V1 = G.dot(V)
+
+np.set_printoptions(precision=6)
+
+>>> print A
+[[ 1.1 1. 2. 7. ]
+ [ 3. 4. 5. 1. ]
+ [ 6. 3. 4.5 2.2]
+ [ 7. 9. 8. 2.3]]
+
+>>> print A1
+[[ 1.1 1. 2. 7. ]
+ [ 5.038805 6.614154 7.251229 1.669176]
+ [ 6. 3. 4.5 2.2 ]
+ [ 5.710557 7.297465 6.034874 1.871858]]
+
+ >>> print A2
+[[ 1.1 3.139139 2. 6.336073]
+ [ 3. 4.11209 5. -0.301199]
+ [ 6. 3.537567 4.5 1.151359]
+ [ 7. 9.26784 8. -0.630206]]
+
+>>> print V1
+[ 3. 5.061923 2. -0.613957]
+
+*/
+
+
+template <typename NumberType>
+void
+test()
+{
+ const std::array<NumberType, 3> csr = {{0.949833, 0.312758, 13.0767}};
+ LAPACKFullMatrix<NumberType> A(4);
+ Vector<NumberType> V(4);
+ const unsigned int i = 1;
+ const unsigned int k = 3;
+
+ A(0, 0) = 1.1;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(0, 3) = 7;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+ A(1, 3) = 1;
+ A(2, 0) = 6;
+ A(2, 1) = 3;
+ A(2, 2) = 4.5;
+ A(2, 3) = 2.2;
+ A(3, 0) = 7;
+ A(3, 1) = 9;
+ A(3, 2) = 8;
+ A(3, 3) = 2.3;
+
+ V(0) = 3;
+ V(1) = 5;
+ V(2) = 2;
+ V(3) = 1;
+
+ LAPACKFullMatrix<NumberType> A1(A), A2(A);
+ Vector<NumberType> V1(V);
+
+ A1.apply_givens_rotation(csr, i, k, true);
+ A2.apply_givens_rotation(csr, i, k, false);
+ V1.apply_givens_rotation(csr, i, k);
+
+ deallog << "A1:" << std::endl;
+ A1.print_formatted(deallog.get_file_stream(), 6, false, 10);
+
+ deallog << "A2:" << std::endl;
+ A2.print_formatted(deallog.get_file_stream(), 6, false, 10);
+
+ deallog << "V1:" << std::endl;
+ V1.print(deallog.get_file_stream(), 6, false);
+}
+
+int
+main()
+{
+ initlog();
+ deallog << std::setprecision(6);
+
+ test<double>();
+}
--- /dev/null
+
+DEAL::A1:
+1.100000 1.000000 2.000000 7.000000
+5.038805 6.614154 7.251229 1.669176
+6.000000 3.000000 4.500000 2.200000
+5.710557 7.297465 6.034874 1.871858
+DEAL::A2:
+1.100000 3.139139 2.000000 6.336073
+3.000000 4.112090 5.000000 -0.301199
+6.000000 3.537567 4.500000 1.151359
+7.000000 9.267840 8.000000 -0.630206
+DEAL::V1:
+3.000000 5.061923 2.000000 -0.613957
projects / dealii.git / commitdiff