From bdeb4c24ebaf46cbc3d4addd3cf37f3654ae2737 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Fri, 21 Sep 2018 12:49:43 +0200 Subject: [PATCH] add functions to apply Givens rotation --- doc/news/changes/minor/20181230DenisDavydov | 4 + include/deal.II/lac/lapack_full_matrix.h | 20 +++ include/deal.II/lac/vector.h | 12 ++ include/deal.II/lac/vector.templates.h | 14 ++ source/lac/lapack_full_matrix.cc | 34 +++++ tests/lapack/full_matrix_37.cc | 138 ++++++++++++++++++++ tests/lapack/full_matrix_37.output | 13 ++ 7 files changed, 235 insertions(+) create mode 100644 doc/news/changes/minor/20181230DenisDavydov create mode 100644 tests/lapack/full_matrix_37.cc create mode 100644 tests/lapack/full_matrix_37.output diff --git a/doc/news/changes/minor/20181230DenisDavydov b/doc/news/changes/minor/20181230DenisDavydov new file mode 100644 index 0000000000..e5cc0b0af5 --- /dev/null +++ b/doc/news/changes/minor/20181230DenisDavydov @@ -0,0 +1,4 @@ +New: Add LAPACKFullMatrix::apply_givens_rotation() and Vector::apply_givens_rotation() to +apply a Givens rotation matrix. +
+(Denis Davydov, 2018/12/30) diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index f80b549486..802a40c9d7 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -175,6 +175,26 @@ public: void rank1_update(const number a, const Vector &v); + /** + * Apply Givens + * rotation + * @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 true, 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 &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 diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h index a6dfa1cdc9..7547f044c0 100644 --- a/include/deal.II/lac/vector.h +++ b/include/deal.II/lac/vector.h @@ -298,6 +298,18 @@ public: void grow_or_shrink(const size_type N); + /** + * Apply Givens + * rotation + * @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 &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. diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index 14f0d85c28..e67ffd71ae 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -59,6 +59,20 @@ Vector::Vector(const Vector &v) +template +void +Vector::apply_givens_rotation(const std::array &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 Vector::Vector(Vector &&v) noexcept : Subscriptor(std::move(v)) diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 2600f50c9f..8141cdc186 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -300,6 +300,40 @@ LAPACKFullMatrix::grow_or_shrink(const size_type n) +template +void +LAPACKFullMatrix::apply_givens_rotation( + const std::array &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 void LAPACKFullMatrix::remove_row_and_column(const size_type row, diff --git a/tests/lapack/full_matrix_37.cc b/tests/lapack/full_matrix_37.cc new file mode 100644 index 0000000000..d40b70b669 --- /dev/null +++ b/tests/lapack/full_matrix_37.cc @@ -0,0 +1,138 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include +#include +#include +#include + +#include +#include + +#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 +void +test() +{ + const std::array csr = {{0.949833, 0.312758, 13.0767}}; + LAPACKFullMatrix A(4); + Vector 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 A1(A), A2(A); + Vector 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(); +} diff --git a/tests/lapack/full_matrix_37.output b/tests/lapack/full_matrix_37.output new file mode 100644 index 0000000000..fb70c9722e --- /dev/null +++ b/tests/lapack/full_matrix_37.output @@ -0,0 +1,13 @@ + +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 -- 2.39.5