]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add functions to apply Givens rotation 7559/head
authorDenis Davydov <davydden@gmail.com>
Fri, 21 Sep 2018 10:49:43 +0000 (12:49 +0200)
committerDenis Davydov <davydden@gmail.com>
Sun, 30 Dec 2018 18:32:28 +0000 (19:32 +0100)
doc/news/changes/minor/20181230DenisDavydov [new file with mode: 0644]
include/deal.II/lac/lapack_full_matrix.h
include/deal.II/lac/vector.h
include/deal.II/lac/vector.templates.h
source/lac/lapack_full_matrix.cc
tests/lapack/full_matrix_37.cc [new file with mode: 0644]
tests/lapack/full_matrix_37.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20181230DenisDavydov b/doc/news/changes/minor/20181230DenisDavydov
new file mode 100644 (file)
index 0000000..e5cc0b0
--- /dev/null
@@ -0,0 +1,4 @@
+New: Add LAPACKFullMatrix::apply_givens_rotation() and Vector::apply_givens_rotation() to
+apply a Givens rotation matrix.
+<br>
+(Denis Davydov, 2018/12/30)
index f80b549486b4dbc92943b305efb5942a4b5afdcb..802a40c9d7e63fc3906e5224b454f89f7420ab72 100644 (file)
@@ -175,6 +175,26 @@ public:
   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
index a6dfa1cdc9e8584fe143dcb7f43a2a0709560572..7547f044c0f270c6963d87fca27da6ddc83d8a1d 100644 (file)
@@ -298,6 +298,18 @@ public:
   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.
index 14f0d85c28dbba7777926d14e9ff1afacbc229f4..e67ffd71ae183f84c068a4b1ff7700d5061aecca 100644 (file)
@@ -59,6 +59,20 @@ Vector<Number>::Vector(const Vector<Number> &v)
 
 
 
+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))
index 2600f50c9f407feae7904a099a00a9f7bfe45cb5..8141cdc186b4b16bbc1eedb8680a8c325103038e 100644 (file)
@@ -300,6 +300,40 @@ LAPACKFullMatrix<number>::grow_or_shrink(const size_type n)
 
 
 
+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,
diff --git a/tests/lapack/full_matrix_37.cc b/tests/lapack/full_matrix_37.cc
new file mode 100644 (file)
index 0000000..d40b70b
--- /dev/null
@@ -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 <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>();
+}
diff --git a/tests/lapack/full_matrix_37.output b/tests/lapack/full_matrix_37.output
new file mode 100644 (file)
index 0000000..fb70c97
--- /dev/null
@@ -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 

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.