]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add QR classes 7222/head
authorDenis Davydov <davydden@gmail.com>
Fri, 21 Sep 2018 12:47:06 +0000 (14:47 +0200)
committerDenis Davydov <davydden@gmail.com>
Tue, 8 Jan 2019 18:20:29 +0000 (19:20 +0100)
21 files changed:
include/deal.II/lac/qr.h [new file with mode: 0644]
tests/lac/qr.cc [new file with mode: 0644]
tests/lac/qr.output [new file with mode: 0644]
tests/lac/qr_02.cc [new file with mode: 0644]
tests/lac/qr_02.output [new file with mode: 0644]
tests/lac/qr_03.cc [new file with mode: 0644]
tests/lac/qr_03.output [new file with mode: 0644]
tests/lac/qr_04.cc [new file with mode: 0644]
tests/lac/qr_04.output [new file with mode: 0644]
tests/lac/qr_05.cc [new file with mode: 0644]
tests/lac/qr_05.output [new file with mode: 0644]
tests/lac/qr_06.cc [new file with mode: 0644]
tests/lac/qr_06.output [new file with mode: 0644]
tests/lac/qr_07.cc [new file with mode: 0644]
tests/lac/qr_07.output [new file with mode: 0644]
tests/lac/qr_08.cc [new file with mode: 0644]
tests/lac/qr_08.output [new file with mode: 0644]
tests/lac/qr_09.cc [new file with mode: 0644]
tests/lac/qr_09.output [new file with mode: 0644]
tests/lac/qr_10.cc [new file with mode: 0644]
tests/lac/qr_10.output [new file with mode: 0644]

diff --git a/include/deal.II/lac/qr.h b/include/deal.II/lac/qr.h
new file mode 100644 (file)
index 0000000..1898d6f
--- /dev/null
@@ -0,0 +1,876 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_qr_h
+#define dealii_qr_h
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/std_cxx14/memory.h>
+
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/lapack_templates.h>
+#include <deal.II/lac/utilities.h>
+
+#include <boost/signals2/signal.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * A base class for thin QR implementations.
+ *
+ * This class and classes derived from it are meant to build $Q$ and $R$
+ * matrices one row/column at a time, i.e., by growing $R$ matrix from an empty
+ * $0\times0$ matrix to $N\timesN$, where $N$ is the number of added column
+ * vectors.
+ *
+ * As a consequence, matrices which have the same number of rows as each vector
+ * (i.e. $Q$ matrix) is stored as a collection of vectors of `VectorType`.
+ */
+template <typename VectorType>
+class BaseQR
+{
+  /**
+   * Number type for R matrix.
+   */
+  typedef typename VectorType::value_type Number;
+
+protected:
+  /**
+   * Default private constructor.
+   */
+  BaseQR();
+
+public:
+  /**
+   * Destructor.
+   */
+  virtual ~BaseQR() = default;
+
+  /**
+   * Append @p column to the QR factorization.
+   * Returns <code>true</code> if the result is successful, i.e.
+   * the columns are linearly independent. Otherwise the @p column
+   * is rejected and the return value is <code>false</code>.
+   */
+  virtual bool
+  append_column(const VectorType &column) = 0;
+
+  /**
+   * Remove a column @p k and update QR factorization.
+   */
+  virtual void
+  remove_column(const unsigned int k = 0) = 0;
+
+  /**
+   * Return size of the subspace.
+   */
+  unsigned int
+  size() const;
+
+  /**
+   * Return the current upper triangular matrix R.
+   */
+  const LAPACKFullMatrix<Number> &
+  get_R() const;
+
+  /**
+   * Solve $Rx=y$. Vectors @p x and @p y should be consistent
+   * with the current size of the subspace.
+   * If @p transpose is <code>true</code>, $R^Tx=y$ is solved instead.
+   */
+  void
+  solve(Vector<Number> &      x,
+        const Vector<Number> &y,
+        const bool            transpose = false) const;
+
+  /**
+   * Set $y = Qx$. The size of $x$ should be consistent with the
+   * size of the R matrix.
+   */
+  virtual void
+  multiply_with_Q(VectorType &y, const Vector<Number> &x) const = 0;
+
+  /**
+   * Set $y = Q^Tx$. The size of $x$ should be consistent with the
+   * size of column vectors.
+   */
+  virtual void
+  multiply_with_QT(Vector<Number> &y, const VectorType &x) const = 0;
+
+  /**
+   * Set $y = QRx$. The size of $x$ should be consistent with the
+   * size of the R matrix.
+   */
+  virtual void
+  multiply_with_A(VectorType &y, const Vector<Number> &x) const = 0;
+
+  /**
+   * Set $y = R^T Q^Tx$. The size of $x$ should be consistent with the
+   * size of column vectors.
+   */
+  virtual void
+  multiply_with_AT(Vector<Number> &y, const VectorType &x) const = 0;
+
+  /**
+   * Connect a slot to retrieve a notification when the Givens rotations
+   * are performed.
+   *
+   * The function takes @p i'th and @p j'th indices of the plane of rotation
+   * and a triplet of numbers @p csr (cosine, sine and radius, see
+   * Utilities::LinearAlgebra::givens_rotation()) which represents the rotation
+   * matrix.
+   */
+  boost::signals2::connection
+  connect_givens_slot(
+    const std::function<void(const unsigned int           i,
+                             const unsigned int           j,
+                             const std::array<Number, 3> &csr)> &slot);
+
+protected:
+  /**
+   * Compute $y=Hx$ where $H$ is the matrix formed by the column vectors stored
+   * by this object.
+   */
+  void
+  multiply_with_cols(VectorType &y, const Vector<Number> &x) const;
+
+  /**
+   * Multiply with transpose columns stored in the object.
+   */
+  void
+  multiply_with_colsT(Vector<Number> &y, const VectorType &x) const;
+
+  /**
+   * A vector of unique pointers to store columns.
+   */
+  std::vector<std::unique_ptr<VectorType>> columns;
+
+  /**
+   * Matrix to store R.
+   */
+  LAPACKFullMatrix<Number> R;
+
+  /**
+   * Current size (number of columns in Q).
+   */
+  unsigned int current_size;
+
+  /**
+   * Signal used to retrieve a notification
+   * when Givens rotations are performed in plane @p i and @p j.
+   */
+  boost::signals2::signal<void(const unsigned int i,
+                               const unsigned int j,
+                               const std::array<Number, 3> &)>
+    givens_signal;
+};
+
+// clang-format off
+/**
+ * A class to compute and store the QR factorization of a matrix represented by a set of column vectors.
+ *
+ * The class is design to update a given (possibly empty) QR factorization
+ * of a matrix $A$ (constructed incrementally by providing its columns)
+ * due to the addition of a new column vector to $A$. This is equivalent to
+ * constructing an orthonormal basis by the Gram-Schmidt procedure.
+ * The class also provides update functionality when the first column
+ * is removed.
+ *
+ * The `VectorType` template argument may either be a parallel and serial vector, and only need
+ * to have basic operations such as additions, scalar product, etc.
+ * It also needs to have a copy-constructor.
+ *
+ * See sections 6.5.2-6.5.3 on pp. 335-338 in
+ * @code{.bib}
+ * @Book{Golub2013,
+ *   title     = {Matrix computations},
+ *   publisher = {Johns Hopkins University Press},
+ *   year      = {2013},
+ *   author    = {Golub, Gene H and Van Loan, Charles F},
+ *   edition   = {4},
+ *  }
+ * @endcode
+ * as well as
+ * @code{.bib}
+ * @article{Daniel1976,
+ *   author   = {Daniel, James W and Gragg, Walter Bill and Kaufman, Linda and Stewart, Gilbert W},
+ *   title    = {{Reorthogonalization and stable algorithms for updating the Gram-Schmidt QR factorization}},
+ *   journal  = {Mathematics of Computation},
+ *   year     = {1976},
+ *   volume   = {30},
+ *   number   = {136},
+ *   pages    = {772--795},
+ * }
+ * @Article{Reichel1990,
+ *   author     = {Reichel, L. and Gragg, W. B.},
+ *   title      = {{Algorithm 686: FORTRAN Subroutines for Updating the QR Decomposition}},
+ *   journal    = {ACM Trans. Math. Softw.},
+ *   year       = {1990},
+ *   volume     = {16},
+ *   number     = {4},
+ *   pages      = {369--377},
+ *   month      = dec,
+ *   issn       = {0098-3500},
+ *   acmid      = {98291},
+ *   address    = {New York, NY, USA},
+ *   doi        = {10.1145/98267.98291},
+ *   issue_date = {Dec. 1990},
+ *   numpages   = {9},
+ *   publisher  = {ACM},
+ *   url        = {http://doi.acm.org/10.1145/98267.98291},
+ *  }
+ * @endcode
+ */
+// clang-format on
+template <typename VectorType>
+class QR : public BaseQR<VectorType>
+{
+public:
+  /**
+   * Number type for R matrix.
+   */
+  typedef typename VectorType::value_type Number;
+
+  /**
+   * Default constructor.
+   */
+  QR();
+
+  /**
+   * Destructor.
+   */
+  virtual ~QR() = default;
+
+  /**
+   * @copyfrom BaseQR::append_column
+   *
+   * @note Currently this function always returns <code>true</code>.
+   */
+  virtual bool
+  append_column(const VectorType &column);
+
+  /**
+   * Remove first column and update QR factorization.
+   *
+   * Starting from the given QR decomposition
+   * $QR= A = [a_1\,\dots a_n], \quad a_i \in {\mathbb R}^m$
+   * we aim at computing factorization of
+   * $\tilde Q \tilde R= \tilde A = [a_2\,\dots a_n], \quad a_i \in {\mathbb
+   * R}^m$.
+   *
+   * The standard approach is to partition $R$ as
+   * \f[
+   * R =
+   * \begin{bmatrix}
+   * r_{11} & w^T \\
+   * 0      & R_{33}
+   * \end{bmatrix}
+   * \f]
+   * It then follows that
+   * \f[
+   * Q^T \tilde A =
+   * \begin{bmatrix}
+   * 0 & w^T \\
+   * 0 & R_{33}
+   * \end{bmatrix}
+   * \f]
+   * is upper Hessenberg where unwanted sub-diagonal elements can be
+   * zeroed by a sequence of Givens rotations.
+   *
+   * Note that $\tilde R^T \tilde R = \tilde A^T \tilde A$,
+   * where the RHS is included in $A^T A = R^T R$. Therefore
+   * $\tilde R$ can be obtained by Cholesky decomposition.
+   */
+  virtual void
+  remove_column(const unsigned int k = 0);
+
+  virtual void
+  multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
+
+  virtual void
+  multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
+
+  virtual void
+  multiply_with_A(VectorType &y, const Vector<Number> &x) const;
+
+  virtual void
+  multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
+
+private:
+  /**
+   * Apply givens rotation in plane @p i @p k to @p Q and @p R so that
+   * <code>R(k,k)</code> is zeroed.
+   *
+   * See Chapter 5.1.9 of Golub 2013, Matrix computations.
+   */
+  void
+  apply_givens_rotation(const unsigned int i, const unsigned int k);
+
+  /**
+   * Temporary vector needed to do Givens rotation of Q.
+   */
+  VectorType tmp;
+};
+
+
+
+/**
+ * A class to obtain the triangular $R$ matrix of the $A=QR$ factorization
+ * together with the matrix $A$ itself. The orthonormal matrix $Q$ is not stored
+ * explicitly, the name of the class.
+ * The multiplication with $Q$ can be represented as $Q=A R^{-1}$, whereas the
+ * multiplication with $Q^T$ is given by $Q^T=R^{-T}A^T$.
+ *
+ * The class is designed to update a given (possibly empty) QR factorization
+ * due to the addition of a new column vector. This is equivalent to
+ * constructing an orthonormal basis by the Gram-Schmidt procedure.
+ * The class also provides update functionality when the column
+ * is removed.
+ *
+ * The `VectorType` template argument may either be a parallel and serial
+ * vector, and only need to have basic operations such as additions, scalar
+ * product, etc. It also needs to have a copy-constructor.
+ */
+template <typename VectorType>
+class ImplicitQR : public BaseQR<VectorType>
+{
+public:
+  /**
+   * Number type for R matrix.
+   */
+  typedef typename VectorType::value_type Number;
+
+  /**
+   * Default constructor.
+   */
+  ImplicitQR();
+
+  /**
+   * Destructor.
+   */
+  virtual ~ImplicitQR() = default;
+
+  virtual bool
+  append_column(const VectorType &column);
+
+  /**
+   * Remove column and update QR factorization.
+   *
+   * Starting from the given QR decomposition
+   * $QR= A = [a_1\,\dots a_n], \quad a_i \in R^m$
+   * we aim at computing factorization of
+   * $\tilde Q \tilde R= \tilde A = [a_2\,\dots a_n], \quad a_i \in R^m$.
+   *
+   * Note that $\tilde R^T \tilde R = \tilde A^T \tilde A$,
+   * where the RHS is included in $A^T A = R^T R$. Therefore
+   * $\tilde R$ can be obtained by Cholesky decomposition.
+   */
+  virtual void
+  remove_column(const unsigned int k = 0);
+
+  virtual void
+  multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
+
+  virtual void
+  multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
+
+  virtual void
+  multiply_with_A(VectorType &y, const Vector<Number> &x) const;
+
+  virtual void
+  multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
+
+  /**
+   * Connect a slot to implement a custom check of linear dependency
+   * during addition of a column.
+   *
+   * Here, @p u is the last column of the to-be R matrix , @p rho
+   * is its diagonal and @p col_norm_sqr is the square of the $l2$ norm of the column.
+   * The function should return <code>true</code> if the new column is
+   * linearly independent.
+   */
+  boost::signals2::connection
+  connect_append_column_slot(
+    const std::function<bool(const Vector<Number> &u,
+                             const Number &        rho2,
+                             const Number &        col_norm_sqr)> &slot);
+
+private:
+  /**
+   * Apply givens rotation in plane @i @p k to zero out $R(k,k)$.
+   */
+  void
+  apply_givens_rotation(const unsigned int i, const unsigned int k);
+
+  /**
+   * Signal used to decide if the new column is linear dependent.
+   *
+   * Here, @p u is the last column of the to-be R matrix , @p rho
+   * is its diagonal and @p col_norm_sqr is the square of the $l2$ norm of the column.
+   * The function should return <code>true</code> if the new column is
+   * linearly independent.
+   */
+  boost::signals2::signal<bool(const Vector<Number> &u,
+                               const Number &        rho,
+                               const Number &        col_norm_sqr)>
+    column_signal;
+};
+
+// -------------------  inline and template functions ----------------
+#ifndef DOXYGEN
+
+
+
+template <typename VectorType>
+BaseQR<VectorType>::BaseQR()
+  : current_size(0)
+{
+  R.set_property(LAPACKSupport::upper_triangular);
+}
+
+
+
+template <typename VectorType>
+unsigned int
+BaseQR<VectorType>::size() const
+{
+  return current_size;
+}
+
+
+
+template <typename VectorType>
+const LAPACKFullMatrix<typename BaseQR<VectorType>::Number> &
+BaseQR<VectorType>::get_R() const
+{
+  return R;
+}
+
+
+
+template <typename VectorType>
+void
+BaseQR<VectorType>::solve(Vector<Number> &      x,
+                          const Vector<Number> &y,
+                          const bool            transpose) const
+{
+  Assert(x.size() == this->current_size,
+         ExcDimensionMismatch(x.size(), this->current_size));
+  Assert(y.size() == this->current_size,
+         ExcDimensionMismatch(y.size(), this->current_size));
+
+  // copy if the two vectors are not the same
+  if (&x != &y)
+    x = y;
+
+  const int lda   = this->current_size;
+  const int ldb   = this->current_size;
+  const int N     = this->current_size;
+  const int n_rhs = 1;
+  int       info  = 0;
+  trtrs("U",
+        transpose ? "T" : "N",
+        "N",
+        &N,
+        &n_rhs,
+        &this->R(0, 0),
+        &lda,
+        &x(0),
+        &ldb,
+        &info);
+}
+
+
+
+template <typename VectorType>
+void
+BaseQR<VectorType>::multiply_with_cols(VectorType &          y,
+                                       const Vector<Number> &x) const
+{
+  Assert(x.size() == this->current_size,
+         ExcDimensionMismatch(x.size(), this->current_size));
+
+  y = 0.;
+  for (unsigned int j = 0; j < this->current_size; ++j)
+    y.add(x[j], *this->columns[j]);
+}
+
+
+
+template <typename VectorType>
+void
+BaseQR<VectorType>::multiply_with_colsT(Vector<Number> &  y,
+                                        const VectorType &x) const
+{
+  Assert(y.size() == this->current_size,
+         ExcDimensionMismatch(y.size(), this->current_size));
+
+  for (unsigned int j = 0; j < this->current_size; ++j)
+    y[j] = (*this->columns[j]) * x;
+}
+
+
+
+template <class VectorType>
+boost::signals2::connection
+BaseQR<VectorType>::connect_givens_slot(
+  const std::function<void(const unsigned int i,
+                           const unsigned int j,
+                           const std::array<Number, 3> &)> &slot)
+{
+  return givens_signal.connect(slot);
+}
+
+
+
+template <class VectorType>
+boost::signals2::connection
+ImplicitQR<VectorType>::connect_append_column_slot(
+  const std::function<bool(const Vector<Number> &u,
+                           const Number &        rho,
+                           const Number &        col_norm_sqr)> &slot)
+{
+  return column_signal.connect(slot);
+}
+
+
+
+template <typename VectorType>
+ImplicitQR<VectorType>::ImplicitQR()
+  : BaseQR<VectorType>()
+{}
+
+
+
+template <typename VectorType>
+bool
+ImplicitQR<VectorType>::append_column(const VectorType &column)
+{
+  if (this->current_size == 0)
+    {
+      this->R.grow_or_shrink(this->current_size + 1);
+      this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
+      this->R(0, 0) = column.l2_norm();
+      ++this->current_size;
+    }
+  else
+    {
+      // first get scalar products with A^T
+      Vector<Number> u(this->current_size);
+      this->multiply_with_AT(u, column);
+
+      // now solve R^T x = (A^T * column)
+      const int lda   = this->current_size;
+      const int ldb   = this->current_size;
+      const int N     = this->current_size;
+      const int n_rhs = 1;
+      int       info  = 0;
+      trtrs(
+        "U", "T", "N", &N, &n_rhs, &this->R(0, 0), &lda, &u(0), &ldb, &info);
+
+      // finally get the diagonal element:
+      // rho2 = |column|^2 - |u|^2
+      const Number column_norm_sqr = column.norm_sqr();
+      const Number rho2            = column_norm_sqr - u.norm_sqr();
+      const bool   linearly_independent =
+        column_signal.empty() ? rho2 > 0 :
+                                column_signal(u, rho2, column_norm_sqr).get();
+
+      // bail out if it turns out to be linearly dependent
+      if (!linearly_independent)
+        return false;
+
+      // at this point we update is successful and we can enlarge R
+      // and store the column:
+      this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
+      this->R.grow_or_shrink(this->current_size + 1);
+      this->R(this->current_size, this->current_size) = std::sqrt(rho2);
+      for (unsigned int i = 0; i < this->current_size; ++i)
+        this->R(i, this->current_size) = u(i);
+
+      this->current_size++;
+    }
+
+  return true;
+}
+
+
+
+template <typename VectorType>
+void
+ImplicitQR<VectorType>::apply_givens_rotation(const unsigned int i,
+                                              const unsigned int k)
+{
+  Assert(i < k, ExcIndexRange(i, 0, k));
+  Assert(k < this->current_size, ExcIndexRange(k, 0, this->current_size));
+  const std::array<Number, 3> csr =
+    dealii::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
+                                                              this->R(k, k));
+
+  // first, set k'th column:
+  this->R(i, k) = csr[2];
+  this->R(k, k) = 0.;
+  // now do the rest:
+  for (unsigned int j = 0; j < this->R.n(); ++j)
+    if (j != k)
+      {
+        const Number t = this->R(i, j);
+        this->R(i, j)  = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
+        this->R(k, j)  = -csr[1] * t + csr[0] * this->R(k, j);
+      }
+
+  if (!this->givens_signal.empty())
+    this->givens_signal(i, k, csr);
+}
+
+
+
+template <typename VectorType>
+void
+ImplicitQR<VectorType>::remove_column(const unsigned int k)
+{
+  // before actually removing a column from Q and resizing R,
+  // apply givens rotations to bring H into upper triangular form:
+  for (unsigned int j = k + 1; j < this->R.n(); ++j)
+    {
+      const unsigned int i = j - 1;
+      apply_givens_rotation(i, j);
+    }
+
+  // remove last row and k-th column
+  --this->current_size;
+  this->R.remove_row_and_column(this->current_size, k);
+
+  // Finally remove the column from A
+  this->columns.erase(this->columns.begin() + k);
+}
+
+
+
+template <typename VectorType>
+void
+ImplicitQR<VectorType>::multiply_with_Q(VectorType &          y,
+                                        const Vector<Number> &x) const
+{
+  // A = QR
+  // A R^{-1} = Q
+  Vector<Number> x1 = x;
+  BaseQR<VectorType>::solve(x1, x1, false);
+  multiply_with_A(y, x1);
+}
+
+
+
+template <typename VectorType>
+void
+ImplicitQR<VectorType>::multiply_with_QT(Vector<Number> &  y,
+                                         const VectorType &x) const
+{
+  // A = QR
+  // A^T = R^T Q^T
+  // {R^T}^{-1} A^T = Q^T
+  multiply_with_AT(y, x);
+  BaseQR<VectorType>::solve(y, y, true);
+}
+
+
+
+template <typename VectorType>
+void
+ImplicitQR<VectorType>::multiply_with_A(VectorType &          y,
+                                        const Vector<Number> &x) const
+{
+  BaseQR<VectorType>::multiply_with_cols(y, x);
+}
+
+
+
+template <typename VectorType>
+void
+ImplicitQR<VectorType>::multiply_with_AT(Vector<Number> &  y,
+                                         const VectorType &x) const
+{
+  BaseQR<VectorType>::multiply_with_colsT(y, x);
+}
+
+
+
+template <typename VectorType>
+QR<VectorType>::QR()
+  : BaseQR<VectorType>()
+{}
+
+
+
+template <typename VectorType>
+bool
+QR<VectorType>::append_column(const VectorType &column)
+{
+  // resize R:
+  this->R.grow_or_shrink(this->current_size + 1);
+  this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
+
+  // now a Gram-Schmidt part: orthonormalize the new column
+  // against everything we have so far:
+  auto &last_col = *this->columns.back();
+  for (unsigned int i = 0; i < this->current_size; ++i)
+    {
+      const auto &i_col              = *this->columns[i];
+      this->R(i, this->current_size) = i_col * last_col;
+      last_col.add(-this->R(i, this->current_size), i_col);
+    }
+
+  this->R(this->current_size, this->current_size) = last_col.l2_norm();
+
+  Assert(this->R(this->current_size, this->current_size) > 0.,
+         ExcDivideByZero());
+  last_col *= 1. / this->R(this->current_size, this->current_size);
+
+  ++this->current_size;
+  return true;
+}
+
+
+
+template <typename VectorType>
+void
+QR<VectorType>::apply_givens_rotation(const unsigned int i,
+                                      const unsigned int k)
+{
+  Assert(i < k, ExcIndexRange(i, 0, k));
+  Assert(k < this->current_size, ExcIndexRange(k, 0, this->current_size));
+  const std::array<Number, 3> csr =
+    dealii::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
+                                                              this->R(k, k));
+
+  // first, set k'th column:
+  this->R(i, k) = csr[2];
+  this->R(k, k) = 0.;
+  // now do the rest:
+  for (unsigned int j = 0; j < this->R.n(); ++j)
+    if (j != k)
+      {
+        const Number t = this->R(i, j);
+        this->R(i, j)  = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
+        this->R(k, j)  = -csr[1] * t + csr[0] * this->R(k, j);
+      }
+
+  // now adjust i,k columns due to multiplication with the
+  // transpose Givens matrix from right:
+  auto &col_i = *this->columns[i];
+  auto &col_k = *this->columns[k];
+  // save column i:
+  tmp = col_i;
+  col_i.sadd(csr[0], csr[1], col_k);
+  col_k.sadd(csr[0], -csr[1], tmp);
+
+  if (!this->givens_signal.empty())
+    this->givens_signal(i, k, csr);
+}
+
+
+
+template <typename VectorType>
+void
+QR<VectorType>::remove_column(const unsigned int k)
+{
+  Assert(k < this->current_size, ExcIndexRange(k, 0, this->current_size));
+  Assert(this->current_size > 0,
+         ExcMessage("Can not remove a column if QR is empty"));
+  // apply a sequence of Givens rotations
+  // see section 6.5 "Updating matrix factorizations" in Golub 2013, Matrix
+  // computations
+
+  // So we want to have QR for \tilde A \in R^{m*(n-1)}
+  // if we remove the column k, we end up with upper Hessenberg matrix
+  //      x x x x x
+  //        x x x x
+  // H =      x x x
+  //          x x x
+  //            x x
+  //              x
+  // where k = 2 (3rd column), m = 7, n = 6
+  //
+  // before actually removing a column from Q and resizing R,
+  // apply givens rotations to bring H into upper triangular form:
+  for (unsigned int j = k + 1; j < this->R.n(); ++j)
+    {
+      const unsigned int i = j - 1;
+      apply_givens_rotation(i, j);
+    }
+
+  // now we can throw away the column from Q and adjust R
+  // since we do thin-QR, after Givens rotations we need to throw
+  // away the last column:
+  const unsigned int size_minus_1 = this->columns.size() - 1;
+  this->columns.erase(this->columns.begin() + size_minus_1);
+
+  // remove last row and k-th column
+  --this->current_size;
+  this->R.remove_row_and_column(this->current_size, k);
+}
+
+
+
+template <typename VectorType>
+void
+QR<VectorType>::multiply_with_Q(VectorType &y, const Vector<Number> &x) const
+{
+  BaseQR<VectorType>::multiply_with_cols(y, x);
+}
+
+
+
+template <typename VectorType>
+void
+QR<VectorType>::multiply_with_QT(Vector<Number> &y, const VectorType &x) const
+{
+  BaseQR<VectorType>::multiply_with_colsT(y, x);
+}
+
+
+
+template <typename VectorType>
+void
+QR<VectorType>::multiply_with_A(VectorType &y, const Vector<Number> &x) const
+{
+  Vector<Number> x1   = x;
+  const int      N    = this->current_size;
+  const int      lda  = N;
+  const int      incx = 1;
+  trmv("U", "N", "N", &N, &this->R(0, 0), &lda, &x1[0], &incx);
+
+  multiply_with_Q(y, x1);
+}
+
+
+
+template <typename VectorType>
+void
+QR<VectorType>::multiply_with_AT(Vector<Number> &y, const VectorType &x) const
+{
+  multiply_with_QT(y, x);
+
+  const int N    = this->current_size;
+  const int lda  = N;
+  const int incx = 1;
+  trmv("U", "T", "N", &N, &this->R(0, 0), &lda, &y[0], &incx);
+}
+
+#endif // no DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/tests/lac/qr.cc b/tests/lac/qr.cc
new file mode 100644 (file)
index 0000000..e93e676
--- /dev/null
@@ -0,0 +1,130 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+/*
+ * Test QR::append_column()
+ */
+
+/*
+ * MWE in Python:
+
+import numpy as np
+from scipy import linalg
+
+A = np.array([[0, 1, 2], [3, 4, 5], [6, 3, 4.5], [7,9, 8], [10,8,0]])
+Q, R = linalg.qr(A)
+# thin QR
+ncol = A.shape[1]
+R = -R[:ncol,:]
+Q = -Q[:, :ncol]
+np.set_printoptions(precision=6)
+
+print Q
+print R
+
+>>> print A
+[[  0.    1.    2. ]
+ [  3.    4.    5. ]
+ [  6.    3.    4.5]
+ [  7.    9.    8. ]
+ [ 10.    8.    0. ]]
+
+>>> print R
+[[ 13.928388  12.420676   7.03599 ]
+ [ -0.         4.089842   4.916632]
+ [ -0.        -0.         6.290594]]
+>>> print Q
+[[-0.        0.244508  0.126831]
+ [ 0.215387  0.32391   0.300765]
+ [ 0.430775 -0.57472   0.682727]
+ [ 0.502571  0.674288  0.182604]
+ [ 0.717958 -0.224343 -0.627689]]
+
+*/
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  QR<VectorType>         qr;
+
+  const unsigned int v_size = 5;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+
+  qr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+
+  qr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+
+  qr.append_column(v);
+
+
+  const LAPACKFullMatrix<double> &R = qr.get_R();
+  std::vector<VectorType>         Q(qr.size());
+  for (unsigned int j = 0; j < qr.size(); ++j)
+    {
+      Vector<double> x(qr.size());
+      x    = 0;
+      x[j] = 1.;
+      Q[j].reinit(v_size);
+      qr.multiply_with_Q(Q[j], x);
+    }
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+  deallog << "R:" << std::endl;
+  R.print_formatted(deallog.get_file_stream(), 6, false, 10);
+  deallog << "Q:" << std::endl;
+
+  for (unsigned int i = 0; i < v_size; ++i)
+    for (unsigned int j = 0; j < size; ++j)
+      {
+        deallog.get_file_stream() << std::setw(9) << Q[j](i);
+        if (j < size - 1)
+          deallog.get_file_stream() << " ";
+        else
+          deallog.get_file_stream() << std::endl;
+      }
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr.output b/tests/lac/qr.output
new file mode 100644 (file)
index 0000000..77e4fa6
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::3
+DEAL::R:
+13.928388  12.420676  7.035990   
+           4.089842   4.916632   
+                      6.290594   
+DEAL::Q:
+0.00000   0.244508  0.126831 
+0.215387  0.323910  0.300765 
+0.430775  -0.574720 0.682727 
+0.502571  0.674288  0.182604 
+0.717958  -0.224343 -0.627689
diff --git a/tests/lac/qr_02.cc b/tests/lac/qr_02.cc
new file mode 100644 (file)
index 0000000..e5bb3cc
--- /dev/null
@@ -0,0 +1,157 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+/*
+ * Same qr.cc but test QR::remove_first_column()
+ */
+
+/*
+ * MWE in Python:
+
+import numpy as np
+from scipy import linalg
+from math import sqrt
+
+A = np.array([[0, 1, 2], [3, 4, 5], [6, 3, 4.5], [7,9, 8], [10,8,0]])
+Q, R = linalg.qr(A)
+Q = -Q
+R = -R
+# first element on diagonal should be
+# sqrt(R[0,1]**2+R[1,1]**2) = 13.076696830622021
+Q1, R1 = linalg.qr_delete(Q, R, 0, 1, 'col', False)
+# thin QR
+ncol = A.shape[1]-1
+R1 = R1[:ncol,:]
+Q1 = Q1[:, :ncol]
+np.set_printoptions(precision=6)
+
+>>> print A
+[[  0.    1.    2. ]
+ [  3.    4.    5. ]
+ [  6.    3.    4.5]
+ [  7.    9.    8. ]
+ [ 10.    8.    0. ]]
+
+>>> print R1
+[[ 13.076697   8.22073 ]
+ [  0.         6.757928]]
+
+>>> print Q1
+[[ 0.076472  0.202924]
+ [ 0.305888  0.367773]
+ [ 0.229416  0.38681 ]
+ [ 0.688247  0.346572]
+ [ 0.611775 -0.744198]]
+
+# first Givens values for plane 0,1 are 0.949833 0.312758 13.0767
+# so we expect the first column to be
+# 0.949833*Q[:,0]+0.312758*Q[:,1]
+*/
+
+
+template <typename VectorType>
+void
+print(const QR<VectorType> &qr, const unsigned int v_size)
+{
+  const LAPACKFullMatrix<double> &R = qr.get_R();
+  std::vector<VectorType>         Q(qr.size());
+  for (unsigned int j = 0; j < qr.size(); ++j)
+    {
+      Vector<double> x(qr.size());
+      x    = 0;
+      x[j] = 1.;
+      Q[j].reinit(v_size);
+      qr.multiply_with_Q(Q[j], x);
+    }
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+  deallog << "R:" << std::endl;
+  R.print_formatted(deallog.get_file_stream(), 6, false, 10);
+  deallog << "Q:" << std::endl;
+
+  for (unsigned int i = 0; i < v_size; ++i)
+    for (unsigned int j = 0; j < size; ++j)
+      {
+        deallog.get_file_stream() << std::setw(9) << Q[j](i);
+        if (j < size - 1)
+          deallog.get_file_stream() << " ";
+        else
+          deallog.get_file_stream() << std::endl;
+      }
+}
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  QR<VectorType>         qr;
+
+  auto print_givens = [](const unsigned int           i,
+                         const unsigned int           j,
+                         const std::array<number, 3> &csr) {
+    deallog.get_file_stream() << "Givens " << i << " " << j << ": " << csr[0]
+                              << " " << csr[1] << " " << csr[2] << std::endl;
+  };
+  qr.connect_givens_slot(print_givens);
+
+  const unsigned int v_size = 5;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+
+  qr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+
+  qr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+
+  qr.append_column(v);
+
+  // up to here ^^^^ exactly the same as in qr_01 test
+  print(qr, v_size);
+
+  // remove first column
+  qr.remove_column();
+
+  print(qr, v_size);
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr_02.output b/tests/lac/qr_02.output
new file mode 100644 (file)
index 0000000..86f3738
--- /dev/null
@@ -0,0 +1,24 @@
+
+DEAL::3
+DEAL::R:
+13.928388  12.420676  7.035990
+           4.089842   4.916632
+                      6.290594
+DEAL::Q:
+0.00000   0.244508  0.126831
+0.215387  0.323910  0.300765
+0.430775  -0.574720 0.682727
+0.502571  0.674288  0.182604
+0.717958  -0.224343 -0.627689
+Givens 0 1: 0.949833 0.312758 13.0767
+Givens 1 2: 0.365410 0.930847 6.75793
+DEAL::2
+DEAL::R:
+13.076697  8.220730
+           6.757928
+DEAL::Q:
+0.076472 0.202924
+0.305888 0.367773
+0.229416 0.38681
+0.688247 0.346572
+0.611775  -0.744198
diff --git a/tests/lac/qr_03.cc b/tests/lac/qr_03.cc
new file mode 100644 (file)
index 0000000..5777677
--- /dev/null
@@ -0,0 +1,132 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+
+/*
+ * Test ImplicitQR::append_column()
+ */
+
+/*
+ * MWE in Python for standard QR:
+
+import numpy as np
+from scipy import linalg
+
+A = np.array([[0, 1, 2], [3, 4, 5], [6, 3, 4.5], [7,9, 8], [10,8,0]])
+Q, R = linalg.qr(A)
+# thin QR
+ncol = A.shape[1]
+R = -R[:ncol,:]
+Q = -Q[:, :ncol]
+np.set_printoptions(precision=6)
+
+print Q
+print R
+
+>>> print A
+[[  0.    1.    2. ]
+ [  3.    4.    5. ]
+ [  6.    3.    4.5]
+ [  7.    9.    8. ]
+ [ 10.    8.    0. ]]
+
+>>> print R
+[[ 13.928388  12.420676   7.03599 ]
+ [ -0.         4.089842   4.916632]
+ [ -0.        -0.         6.290594]]
+>>> print Q
+[[-0.        0.244508  0.126831]
+ [ 0.215387  0.32391   0.300765]
+ [ 0.430775 -0.57472   0.682727]
+ [ 0.502571  0.674288  0.182604]
+ [ 0.717958 -0.224343 -0.627689]]
+
+*/
+
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  ImplicitQR<VectorType> qr;
+
+  const unsigned int v_size = 5;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+
+  qr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+
+  qr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+
+  qr.append_column(v);
+
+
+  const LAPACKFullMatrix<double> &R = qr.get_R();
+  std::vector<VectorType>         Q(R.n());
+  for (unsigned int j = 0; j < R.n(); ++j)
+    {
+      Vector<double> x(R.n());
+      x    = 0;
+      x[j] = 1.;
+      Q[j].reinit(v);
+      qr.multiply_with_A(Q[j], x);
+    }
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+  deallog << "R:" << std::endl;
+  R.print_formatted(deallog.get_file_stream(), 6, false, 10);
+  deallog << "A:" << std::endl;
+
+  for (unsigned int i = 0; i < v_size; ++i)
+    for (unsigned int j = 0; j < size; ++j)
+      {
+        deallog.get_file_stream() << std::setw(9) << Q[j](i);
+        if (j < size - 1)
+          deallog.get_file_stream() << " ";
+        else
+          deallog.get_file_stream() << std::endl;
+      }
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr_03.output b/tests/lac/qr_03.output
new file mode 100644 (file)
index 0000000..eec779f
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::3
+DEAL::R:
+13.928388  12.420676  7.035990   
+           4.089842   4.916632   
+                      6.290594   
+DEAL::A:
+0.00000   1.00000   2.00000  
+3.00000   4.00000   5.00000  
+6.00000   3.00000   4.50000  
+7.00000   9.00000   8.00000  
+10.0000   8.00000   0.00000  
diff --git a/tests/lac/qr_04.cc b/tests/lac/qr_04.cc
new file mode 100644 (file)
index 0000000..f3e7815
--- /dev/null
@@ -0,0 +1,158 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+/*
+ * Same qr_02.cc but test ImplicitQR::remove_column()
+ */
+
+/*
+ * MWE in Python:
+
+import numpy as np
+from scipy import linalg
+from math import sqrt
+
+A = np.array([[0, 1, 2], [3, 4, 5], [6, 3, 4.5], [7,9, 8], [10,8,0]])
+Q, R = linalg.qr(A)
+Q = -Q
+R = -R
+# first element on diagonal should be
+# sqrt(R[0,1]**2+R[1,1]**2) = 13.076696830622021
+Q1, R1 = linalg.qr_delete(Q, R, 0, 1, 'col', False)
+# thin QR
+ncol = A.shape[1]-1
+R1 = R1[:ncol,:]
+Q1 = Q1[:, :ncol]
+np.set_printoptions(precision=6)
+
+>>> print A
+[[  0.    1.    2. ]
+ [  3.    4.    5. ]
+ [  6.    3.    4.5]
+ [  7.    9.    8. ]
+ [ 10.    8.    0. ]]
+
+>>> print R1
+[[ 13.076697   8.22073 ]
+ [  0.         6.757928]]
+
+>>> print Q1
+[[ 0.076472  0.202924]
+ [ 0.305888  0.367773]
+ [ 0.229416  0.38681 ]
+ [ 0.688247  0.346572]
+ [ 0.611775 -0.744198]]
+
+# first Givens values for plane 0,1 are 0.949833 0.312758 13.0767
+# so we expect the first column to be
+# 0.949833*Q[:,0]+0.312758*Q[:,1]
+*/
+
+
+template <typename VectorType>
+void
+print(const ImplicitQR<VectorType> &qr, const unsigned int col_size)
+{
+  const LAPACKFullMatrix<double> &R = qr.get_R();
+  std::vector<VectorType>         A(R.n());
+  for (unsigned int j = 0; j < R.n(); ++j)
+    {
+      Vector<double> x(R.n());
+      x    = 0;
+      x[j] = 1.;
+      A[j].reinit(col_size);
+      qr.multiply_with_A(A[j], x);
+    }
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+  deallog << "R:" << std::endl;
+  R.print_formatted(deallog.get_file_stream(), 6, false, 10);
+  deallog << "A:" << std::endl;
+
+  for (unsigned int i = 0; i < col_size; ++i)
+    for (unsigned int j = 0; j < size; ++j)
+      {
+        deallog.get_file_stream() << std::setw(9) << A[j](i);
+        if (j < size - 1)
+          deallog.get_file_stream() << " ";
+        else
+          deallog.get_file_stream() << std::endl;
+      }
+}
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  ImplicitQR<VectorType> qr;
+
+  auto print_givens = [](const unsigned int           i,
+                         const unsigned int           j,
+                         const std::array<number, 3> &csr) {
+    deallog.get_file_stream() << "Givens " << i << " " << j << ": " << csr[0]
+                              << " " << csr[1] << " " << csr[2] << std::endl;
+  };
+  qr.connect_givens_slot(print_givens);
+
+
+  const unsigned int v_size = 5;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+
+  qr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+
+  qr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+
+  qr.append_column(v);
+
+  // up to here ^^^^ exactly the same as in qr_03 test
+  print(qr, v_size);
+
+  // remove first column
+  qr.remove_column();
+
+  print(qr, v_size);
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr_04.output b/tests/lac/qr_04.output
new file mode 100644 (file)
index 0000000..ab595f1
--- /dev/null
@@ -0,0 +1,24 @@
+
+DEAL::3
+DEAL::R:
+13.928388  12.420676  7.035990
+           4.089842   4.916632
+                      6.290594
+DEAL::A:
+0.00000   1.00000   2.00000
+3.00000   4.00000   5.00000
+6.00000   3.00000   4.50000
+7.00000   9.00000   8.00000
+10.0000   8.00000   0.00000
+Givens 0 1: 0.949833 0.312758 13.0767
+Givens 1 2: 0.365410 0.930847 6.75793
+DEAL::2
+DEAL::R:
+13.076697  8.220730
+           6.757928
+DEAL::A:
+1.00000   2.00000
+4.00000   5.00000
+3.00000   4.50000
+9.00000   8.00000
+8.00000   0.00000
diff --git a/tests/lac/qr_05.cc b/tests/lac/qr_05.cc
new file mode 100644 (file)
index 0000000..5d52ba1
--- /dev/null
@@ -0,0 +1,113 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+/*
+ * Test BaseQR::solve()
+ * QR is the same as in qr.cc
+ */
+
+/*
+ * MWE in Python:
+
+import numpy as np
+from scipy import linalg
+
+A = np.array([[0, 1, 2], [3, 4, 5], [6, 3, 4.5], [7,9, 8], [10,8,0]])
+Q, R = linalg.qr(A)
+# thin QR
+ncol = A.shape[1]
+R = -R[:ncol,:]
+np.set_printoptions(precision=6)
+d = np.array([1,3,5])
+y = linalg.solve_triangular(R, d)
+y2 = linalg.solve_triangular(R, d, trans='T')
+
+>>> print R
+[[ 13.928388  12.420676   7.03599 ]
+ [ -0.         4.089842   4.916632]
+ [ -0.        -0.         6.290594]]
+
+>>> print y
+[-0.131756 -0.221995  0.794838]
+
+>>> print y2
+[ 0.071796  0.515484  0.31164 ]
+*/
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  QR<VectorType>         qr;
+
+  const unsigned int v_size = 5;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+
+  qr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+
+  qr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+
+  qr.append_column(v);
+
+
+  const LAPACKFullMatrix<double> &R = qr.get_R();
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+  deallog << "R:" << std::endl;
+  R.print_formatted(deallog.get_file_stream(), 6, false, 10);
+
+  Vector<double> x(size), y(size);
+  y[0] = 1.;
+  y[1] = 3.;
+  y[2] = 5.;
+
+  qr.solve(x, y, false);
+  x.print(deallog.get_file_stream(), 6, false);
+
+  qr.solve(x, y, true);
+  x.print(deallog.get_file_stream(), 6, false);
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr_05.output b/tests/lac/qr_05.output
new file mode 100644 (file)
index 0000000..fd177ab
--- /dev/null
@@ -0,0 +1,8 @@
+
+DEAL::3
+DEAL::R:
+13.928388  12.420676  7.035990   
+           4.089842   4.916632   
+                      6.290594   
+-0.131756 -0.221995 0.794838 
+0.071796 0.515484 0.311640 
diff --git a/tests/lac/qr_06.cc b/tests/lac/qr_06.cc
new file mode 100644 (file)
index 0000000..89c95ac
--- /dev/null
@@ -0,0 +1,129 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+/*
+ * Test QR::multiply_with_XY()
+ */
+
+/*
+ * MWE in Python:
+
+import numpy as np
+from scipy import linalg
+
+A = np.array([[0, 1, 2], [3, 4, 5], [6, 3, 4.5], [7,9, 8], [10,8,0]])
+Q, R = linalg.qr(A)
+# thin QR
+ncol = A.shape[1]
+R = -R[:ncol,:]
+Q = -Q[:, :ncol]
+
+d1 = np.array([1,3,5,7,9])
+d2 = np.array([10,9,9])
+
+np.set_printoptions(precision=6)
+
+>>> A.dot(d2)
+array([  27. ,  111. ,  127.5,  223. ,  172. ])
+
+>>> A.T.dot(d1)
+array([ 178. ,  163. ,   95.5])
+
+>>> Q.dot(d2)
+array([  3.342054,   7.775949,   5.279812,  12.737741,  -0.488702])
+
+>>> Q.T.dot(d1)
+array([ 12.779655,   1.043571,   0.071793])
+
+*/
+
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  QR<VectorType>         qr;
+
+  const unsigned int v_size = 5;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+
+  qr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+
+  qr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+
+  qr.append_column(v);
+
+  Vector<number> d1(5), d2(3), res1(3), res2(5);
+  d1[0] = 1;
+  d1[1] = 3;
+  d1[2] = 5;
+  d1[3] = 7;
+  d1[4] = 9;
+
+  d2[0] = 10;
+  d2[1] = 9;
+  d2[2] = 9;
+
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+
+  deallog << "A:" << std::endl;
+  qr.multiply_with_A(res2, d2);
+  res2.print(deallog.get_file_stream(), 6, false);
+
+  deallog << "AT:" << std::endl;
+  qr.multiply_with_AT(res1, d1);
+  res1.print(deallog.get_file_stream(), 6, false);
+
+  deallog << "Q:" << std::endl;
+  qr.multiply_with_Q(res2, d2);
+  res2.print(deallog.get_file_stream(), 6, false);
+
+  deallog << "QT:" << std::endl;
+  qr.multiply_with_QT(res1, d1);
+  res1.print(deallog.get_file_stream(), 6, false);
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr_06.output b/tests/lac/qr_06.output
new file mode 100644 (file)
index 0000000..ee397c8
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::3
+DEAL::A:
+27.000000 111.000000 127.500000 223.000000 172.000000 
+DEAL::AT:
+178.000000 163.000000 95.500000 
+DEAL::Q:
+3.342054 7.775949 5.279812 12.737741 -0.488702 
+DEAL::QT:
+12.779655 1.043571 0.071793 
diff --git a/tests/lac/qr_07.cc b/tests/lac/qr_07.cc
new file mode 100644 (file)
index 0000000..c9c24bd
--- /dev/null
@@ -0,0 +1,130 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+/*
+ * Test ImplicitQR::multiply_with_XY()
+ */
+
+/*
+ * MWE in Python:
+
+import numpy as np
+from scipy import linalg
+
+A = np.array([[0, 1, 2], [3, 4, 5], [6, 3, 4.5], [7,9, 8], [10,8,0]])
+Q, R = linalg.qr(A)
+# thin QR
+ncol = A.shape[1]
+R = -R[:ncol,:]
+Q = -Q[:, :ncol]
+
+d1 = np.array([1,3,5,7,9])
+d2 = np.array([10,9,9])
+
+np.set_printoptions(precision=6)
+
+>>> A.dot(d2)
+array([  27. ,  111. ,  127.5,  223. ,  172. ])
+
+>>> A.T.dot(d1)
+array([ 178. ,  163. ,   95.5])
+
+>>> Q.dot(d2)
+array([  3.342054,   7.775949,   5.279812,  12.737741,  -0.488702])
+
+>>> Q.T.dot(d1)
+array([ 12.779655,   1.043571,   0.071793])
+
+*/
+
+
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  ImplicitQR<VectorType> qr;
+
+  const unsigned int v_size = 5;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+
+  qr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+
+  qr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+
+  qr.append_column(v);
+
+  Vector<number> d1(5), d2(3), res1(3), res2(5);
+  d1[0] = 1;
+  d1[1] = 3;
+  d1[2] = 5;
+  d1[3] = 7;
+  d1[4] = 9;
+
+  d2[0] = 10;
+  d2[1] = 9;
+  d2[2] = 9;
+
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+
+  deallog << "A:" << std::endl;
+  qr.multiply_with_A(res2, d2);
+  res2.print(deallog.get_file_stream(), 6, false);
+
+  deallog << "AT:" << std::endl;
+  qr.multiply_with_AT(res1, d1);
+  res1.print(deallog.get_file_stream(), 6, false);
+
+  deallog << "Q:" << std::endl;
+  qr.multiply_with_Q(res2, d2);
+  res2.print(deallog.get_file_stream(), 6, false);
+
+  deallog << "QT:" << std::endl;
+  qr.multiply_with_QT(res1, d1);
+  res1.print(deallog.get_file_stream(), 6, false);
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr_07.output b/tests/lac/qr_07.output
new file mode 100644 (file)
index 0000000..ee397c8
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::3
+DEAL::A:
+27.000000 111.000000 127.500000 223.000000 172.000000 
+DEAL::AT:
+178.000000 163.000000 95.500000 
+DEAL::Q:
+3.342054 7.775949 5.279812 12.737741 -0.488702 
+DEAL::QT:
+12.779655 1.043571 0.071793 
diff --git a/tests/lac/qr_08.cc b/tests/lac/qr_08.cc
new file mode 100644 (file)
index 0000000..6b88b09
--- /dev/null
@@ -0,0 +1,287 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+/*
+ * Use signal for givens rotation to update factorizations Q^T H Q and Q^T u
+ * when the first column is removed from the subspace
+ */
+
+/*
+ * MWE in Python:
+
+import numpy as np
+from scipy import linalg
+from math import sqrt
+
+# vector subspace of R^6:
+A = np.array([[0, 1, 2, 7], [3, 4, 5, 1], [6, 3, 4.5, 2.2], [7,9, 8, 0],
+[10,8,0, 1.1], [7, 1, 9, 5]]) Q, R = linalg.qr(A) Q = -Q R = -R ncol =
+A.shape[1] # thin QR: R1 = R[:ncol,:] Q1 = Q[:, :ncol]
+
+# 1D Laplace matrix and some vector to be represented in the subspace
+M =
+np.array([[2,-1,0,0,0,0],[-1,2,-1,0,0,0],[0,-1,2,-1,0,0],[0,0,-1,2,-1,0],[0,0,0,-1,2,-1],[0,0,0,0,-1,2]])
+V = np.array([1,2,2,1,5,9])
+
+# projection
+H1 = Q1.T.dot(M.dot(Q1))
+V1 = Q1.T.dot(V)
+
+# now delete the first column
+Q2, R2 = linalg.qr_delete(Q, R, 0, 1, 'col', False)
+R2 = R2[:ncol-1,:]
+Q2 = Q2[:, :ncol-1]
+
+# updated projections:
+H2 = Q2.T.dot(M.dot(Q2))
+V2 = Q2.T.dot(V)
+
+np.set_printoptions(precision=6)
+
+>>> print A
+[[  0.    1.    2.    7. ]
+ [  3.    4.    5.    1. ]
+ [  6.    3.    4.5   2.2]
+ [  7.    9.    8.    0. ]
+ [ 10.    8.    0.    1.1]
+ [  7.    1.    9.    5. ]]
+
+>>> print Q1
+[[-0.        0.160817  0.221587 -0.955862]
+ [ 0.19245   0.285897  0.335821  0.033159]
+ [ 0.3849   -0.232291  0.045103  0.005283]
+ [ 0.44905   0.613487  0.388792  0.239359]
+ [ 0.6415    0.095299 -0.70425  -0.16568 ]
+ [ 0.44905  -0.673048  0.434697 -0.021413]]
+
+>>> print R1
+[[ 15.588457  11.547005  10.328155   3.990132]
+ [ -0.         6.218253  -0.443735  -2.359839]
+ [ -0.        -0.         9.347851   3.384966]
+ [ -0.        -0.        -0.        -6.935562]]
+
+>>> print H1
+[[ 0.353909 -0.275486 -0.246262  0.084659]
+ [-0.275486  2.337236 -0.110024  0.295859]
+ [-0.246262 -0.110024  2.945693  0.587455]
+ [ 0.084659  0.295859  0.587455  2.132729]]
+
+>>> print V1
+[ 8.852704 -4.699427  1.76325  -1.660732]
+
+>>> print Q2
+[[ 0.076249  0.123157  0.867562]
+ [ 0.304997  0.213292 -0.108274]
+ [ 0.228748  0.229803  0.073457]
+ [ 0.686244  0.177292 -0.3507  ]
+ [ 0.609994 -0.504539  0.294887]
+ [ 0.076249  0.774943  0.142366]]
+
+>>> print H2
+[[ 0.569767 -0.567113 -0.321939]
+ [-0.567113  2.728827 -0.480842]
+ [-0.321939 -0.480842  2.378169]]
+
+>>> print V2
+[ 5.566198  5.638437  3.202953]
+*/
+
+template <typename VectorType>
+void
+print(const QR<VectorType> &qr, const unsigned int col_size)
+{
+  const LAPACKFullMatrix<double> &R = qr.get_R();
+  std::vector<VectorType>         Q(R.n());
+  for (unsigned int j = 0; j < R.n(); ++j)
+    {
+      Vector<double> x(R.n());
+      x    = 0;
+      x[j] = 1.;
+      Q[j].reinit(col_size);
+      qr.multiply_with_Q(Q[j], x);
+    }
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+  deallog << "R:" << std::endl;
+  R.print_formatted(deallog.get_file_stream(), 6, false, 10);
+  deallog << "Q:" << std::endl;
+
+  for (unsigned int i = 0; i < col_size; ++i)
+    for (unsigned int j = 0; j < size; ++j)
+      {
+        deallog.get_file_stream() << std::setw(9) << Q[j](i);
+        if (j < size - 1)
+          deallog.get_file_stream() << " ";
+        else
+          deallog.get_file_stream() << std::endl;
+      }
+}
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  QR<VectorType>         qr;
+
+  const unsigned int v_size = 6;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+  v[5] = 7.;
+
+  qr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+  v[5] = 1.;
+
+  qr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+  v[5] = 9.;
+
+  qr.append_column(v);
+
+  v[0] = 7;
+  v[1] = 1;
+  v[2] = 2.2;
+  v[3] = 0;
+  v[4] = 1.1;
+  v[5] = 5.;
+
+  qr.append_column(v);
+
+  print(qr, v_size);
+
+  // matrix to be projected:
+  FullMatrix<number> M(v_size);
+  M       = 0.;
+  M(0, 0) = 2.;
+  M(0, 1) = -1;
+  M(1, 0) = -1;
+  M(1, 1) = 2;
+  M(1, 2) = -1;
+  M(2, 1) = -1;
+  M(2, 2) = 2;
+  M(2, 3) = -1;
+  M(3, 2) = -1;
+  M(3, 3) = 2;
+  M(3, 4) = -1;
+  M(4, 3) = -1;
+  M(4, 4) = 2;
+  M(4, 5) = -1;
+  M(5, 4) = -1;
+  M(5, 5) = 2;
+
+  // vector to be projected:
+  Vector<number> V(v_size);
+  V[0] = 1;
+  V[1] = 2;
+  V[2] = 2;
+  V[3] = 1;
+  V[4] = 5;
+  V[5] = 9;
+
+  // reduced matrix
+  // H_{ij} = Q(k,i) M(k,l) Q(l,j)
+  std::vector<VectorType> Q(qr.size());
+  for (unsigned int j = 0; j < qr.size(); ++j)
+    {
+      Vector<double> x(qr.size());
+      x    = 0;
+      x[j] = 1.;
+      Q[j].reinit(v_size);
+      qr.multiply_with_Q(Q[j], x);
+    }
+
+  LAPACKFullMatrix<number> H1(4);
+  for (unsigned int i = 0; i < 4; ++i)
+    for (unsigned int j = 0; j < 4; ++j)
+      {
+        const VectorType &qi  = Q[i];
+        const VectorType &qj  = Q[j];
+        number            res = 0.;
+        for (unsigned int k = 0; k < v_size; ++k)
+          for (unsigned int l = 0; l < v_size; ++l)
+            res += qi[k] * M(k, l) * qj[l];
+
+        H1(i, j) = res;
+      }
+
+  deallog.get_file_stream() << std::endl;
+  deallog << "H1:" << std::endl;
+  H1.print_formatted(deallog.get_file_stream(), 6, false, 10);
+
+  Vector<number> V1(4);
+  qr.multiply_with_QT(V1, V);
+  deallog << "V1:" << std::endl;
+  V1.print(deallog.get_file_stream(), 6, false);
+
+  // remove first column
+  auto update_factorization = [&](const unsigned int           i,
+                                  const unsigned int           k,
+                                  const std::array<number, 3> &csr) {
+    H1.apply_givens_rotation(csr, i, k, true);
+    H1.apply_givens_rotation(csr, i, k, false);
+    V1.apply_givens_rotation(csr, i, k);
+  };
+
+  qr.connect_givens_slot(update_factorization);
+  qr.remove_column();
+
+  // remove one row (and column)
+  Vector<number>           V2(3);
+  LAPACKFullMatrix<number> H2(3);
+  for (unsigned int i = 0; i < 3; ++i)
+    {
+      V2(i) = V1(i);
+      for (unsigned int j = 0; j < 3; ++j)
+        H2(i, j) = H1(i, j);
+    }
+
+
+  deallog.get_file_stream() << std::endl;
+  deallog << "H2:" << std::endl;
+  H2.print_formatted(deallog.get_file_stream(), 6, false, 10);
+
+  deallog << "V2:" << std::endl;
+  V2.print(deallog.get_file_stream(), 6, false);
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr_08.output b/tests/lac/qr_08.output
new file mode 100644 (file)
index 0000000..6740212
--- /dev/null
@@ -0,0 +1,29 @@
+
+DEAL::4
+DEAL::R:
+15.588457  11.547005  10.328155  3.990132   
+           6.218253   -0.443735  -2.359839  
+                      9.347851   3.384966   
+                                 6.935562   
+DEAL::Q:
+0.00000   0.160817  0.221587  0.955862 
+0.192450  0.285897  0.335821  -0.0331588
+0.384900  -0.232291 0.0451029 -0.00528347
+0.449050  0.613487  0.388792  -0.239359
+0.641500  0.0952989 -0.704250 0.165680 
+0.449050  -0.673048 0.434697  0.0214128
+
+DEAL::H1:
+0.353909   -0.275486  -0.246262  -0.084659  
+-0.275486  2.337236   -0.110024  -0.295859  
+-0.246262  -0.110024  2.945693   -0.587455  
+-0.084659  -0.295859  -0.587455  2.132729   
+DEAL::V1:
+8.852704 -4.699427 1.763250 1.660732 
+
+DEAL::H2:
+0.569767   -0.567113  -0.321939  
+-0.567113  2.728827   -0.480842  
+-0.321939  -0.480842  2.378169   
+DEAL::V2:
+5.566198 5.638437 3.202953 
diff --git a/tests/lac/qr_09.cc b/tests/lac/qr_09.cc
new file mode 100644 (file)
index 0000000..3d90c9f
--- /dev/null
@@ -0,0 +1,276 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+/*
+ * Same as qr_08, but test ImplicitQR class, namely
+ * use signal for givens rotation to update factorizations Q^T H Q and Q^T u
+ * when the first column is removed from the subspace.
+ * We use QR class to build Q and project matrix and vector on the subspace,
+ * but attach signal to ImplicitQR.
+ */
+
+/*
+ * MWE in Python:
+
+import numpy as np
+from scipy import linalg
+from math import sqrt
+
+# vector subspace of R^6:
+A = np.array([[0, 1, 2, 7], [3, 4, 5, 1], [6, 3, 4.5, 2.2], [7,9, 8, 0],
+[10,8,0, 1.1], [7, 1, 9, 5]]) Q, R = linalg.qr(A) Q = -Q R = -R ncol =
+A.shape[1] # thin QR: R1 = R[:ncol,:] Q1 = Q[:, :ncol]
+
+# 1D Laplace matrix and some vector to be represented in the subspace
+M =
+np.array([[2,-1,0,0,0,0],[-1,2,-1,0,0,0],[0,-1,2,-1,0,0],[0,0,-1,2,-1,0],[0,0,0,-1,2,-1],[0,0,0,0,-1,2]])
+V = np.array([1,2,2,1,5,9])
+
+# projection
+H1 = Q1.T.dot(M.dot(Q1))
+V1 = Q1.T.dot(V)
+
+# now delete the first column
+Q2, R2 = linalg.qr_delete(Q, R, 0, 1, 'col', False)
+R2 = R2[:ncol-1,:]
+Q2 = Q2[:, :ncol-1]
+
+# updated projections:
+H2 = Q2.T.dot(M.dot(Q2))
+V2 = Q2.T.dot(V)
+
+np.set_printoptions(precision=6)
+
+>>> print A
+[[  0.    1.    2.    7. ]
+ [  3.    4.    5.    1. ]
+ [  6.    3.    4.5   2.2]
+ [  7.    9.    8.    0. ]
+ [ 10.    8.    0.    1.1]
+ [  7.    1.    9.    5. ]]
+
+>>> print Q1
+[[-0.        0.160817  0.221587 -0.955862]
+ [ 0.19245   0.285897  0.335821  0.033159]
+ [ 0.3849   -0.232291  0.045103  0.005283]
+ [ 0.44905   0.613487  0.388792  0.239359]
+ [ 0.6415    0.095299 -0.70425  -0.16568 ]
+ [ 0.44905  -0.673048  0.434697 -0.021413]]
+
+>>> print R1
+[[ 15.588457  11.547005  10.328155   3.990132]
+ [ -0.         6.218253  -0.443735  -2.359839]
+ [ -0.        -0.         9.347851   3.384966]
+ [ -0.        -0.        -0.        -6.935562]]
+
+>>> print H1
+[[ 0.353909 -0.275486 -0.246262  0.084659]
+ [-0.275486  2.337236 -0.110024  0.295859]
+ [-0.246262 -0.110024  2.945693  0.587455]
+ [ 0.084659  0.295859  0.587455  2.132729]]
+
+>>> print V1
+[ 8.852704 -4.699427  1.76325  -1.660732]
+
+>>> print Q2
+[[ 0.076249  0.123157  0.867562]
+ [ 0.304997  0.213292 -0.108274]
+ [ 0.228748  0.229803  0.073457]
+ [ 0.686244  0.177292 -0.3507  ]
+ [ 0.609994 -0.504539  0.294887]
+ [ 0.076249  0.774943  0.142366]]
+
+>>> print H2
+[[ 0.569767 -0.567113 -0.321939]
+ [-0.567113  2.728827 -0.480842]
+ [-0.321939 -0.480842  2.378169]]
+
+>>> print V2
+[ 5.566198  5.638437  3.202953]
+*/
+
+
+template <typename VectorType>
+void
+print(const ImplicitQR<VectorType> &qr)
+{
+  const LAPACKFullMatrix<double> &R = qr.get_R();
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+  deallog << "R:" << std::endl;
+  R.print_formatted(deallog.get_file_stream(), 6, false, 10);
+}
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  QR<VectorType>         qr;
+  ImplicitQR<VectorType> pqr;
+
+  const unsigned int v_size = 6;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+  v[5] = 7.;
+
+  qr.append_column(v);
+  pqr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+  v[5] = 1.;
+
+  qr.append_column(v);
+  pqr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+  v[5] = 9.;
+
+  qr.append_column(v);
+  pqr.append_column(v);
+
+  v[0] = 7;
+  v[1] = 1;
+  v[2] = 2.2;
+  v[3] = 0;
+  v[4] = 1.1;
+  v[5] = 5.;
+
+  qr.append_column(v);
+  pqr.append_column(v);
+
+  print(pqr);
+
+  // matrix to be projected:
+  FullMatrix<number> M(v_size);
+  M       = 0.;
+  M(0, 0) = 2.;
+  M(0, 1) = -1;
+  M(1, 0) = -1;
+  M(1, 1) = 2;
+  M(1, 2) = -1;
+  M(2, 1) = -1;
+  M(2, 2) = 2;
+  M(2, 3) = -1;
+  M(3, 2) = -1;
+  M(3, 3) = 2;
+  M(3, 4) = -1;
+  M(4, 3) = -1;
+  M(4, 4) = 2;
+  M(4, 5) = -1;
+  M(5, 4) = -1;
+  M(5, 5) = 2;
+
+  // vector to be projected:
+  Vector<number> V(v_size);
+  V[0] = 1;
+  V[1] = 2;
+  V[2] = 2;
+  V[3] = 1;
+  V[4] = 5;
+  V[5] = 9;
+
+  // reduced matrix
+  // H_{ij} = Q(k,i) M(k,l) Q(l,j)
+  std::vector<VectorType> Q(qr.size());
+  for (unsigned int j = 0; j < qr.size(); ++j)
+    {
+      Vector<double> x(qr.size());
+      x    = 0;
+      x[j] = 1.;
+      Q[j].reinit(v_size);
+      qr.multiply_with_Q(Q[j], x);
+    }
+
+  LAPACKFullMatrix<number> H1(4);
+  for (unsigned int i = 0; i < 4; ++i)
+    for (unsigned int j = 0; j < 4; ++j)
+      {
+        const VectorType &qi  = Q[i];
+        const VectorType &qj  = Q[j];
+        number            res = 0.;
+        for (unsigned int k = 0; k < v_size; ++k)
+          for (unsigned int l = 0; l < v_size; ++l)
+            res += qi[k] * M(k, l) * qj[l];
+
+        H1(i, j) = res;
+      }
+
+  deallog.get_file_stream() << std::endl;
+  deallog << "H1:" << std::endl;
+  H1.print_formatted(deallog.get_file_stream(), 6, false, 10);
+
+  Vector<number> V1(4);
+  qr.multiply_with_QT(V1, V);
+  deallog << "V1:" << std::endl;
+  V1.print(deallog.get_file_stream(), 6, false);
+
+  // remove first column
+  auto update_factorization = [&](const unsigned int           i,
+                                  const unsigned int           k,
+                                  const std::array<number, 3> &csr) {
+    H1.apply_givens_rotation(csr, i, k, true);
+    H1.apply_givens_rotation(csr, i, k, false);
+    V1.apply_givens_rotation(csr, i, k);
+  };
+
+  pqr.connect_givens_slot(update_factorization);
+  pqr.remove_column();
+
+  // remove one row (and column)
+  Vector<number>           V2(3);
+  LAPACKFullMatrix<number> H2(3);
+  for (unsigned int i = 0; i < 3; ++i)
+    {
+      V2(i) = V1(i);
+      for (unsigned int j = 0; j < 3; ++j)
+        H2(i, j) = H1(i, j);
+    }
+
+
+  deallog.get_file_stream() << std::endl;
+  deallog << "H2:" << std::endl;
+  H2.print_formatted(deallog.get_file_stream(), 6, false, 10);
+
+  deallog << "V2:" << std::endl;
+  V2.print(deallog.get_file_stream(), 6, false);
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr_09.output b/tests/lac/qr_09.output
new file mode 100644 (file)
index 0000000..64b93fc
--- /dev/null
@@ -0,0 +1,22 @@
+
+DEAL::4
+DEAL::R:
+15.588457  11.547005  10.328155  3.990132   
+           6.218253   -0.443735  -2.359839  
+                      9.347851   3.384966   
+                                 6.935562   
+
+DEAL::H1:
+0.353909   -0.275486  -0.246262  -0.084659  
+-0.275486  2.337236   -0.110024  -0.295859  
+-0.246262  -0.110024  2.945693   -0.587455  
+-0.084659  -0.295859  -0.587455  2.132729   
+DEAL::V1:
+8.852704 -4.699427 1.763250 1.660732 
+
+DEAL::H2:
+0.569767   -0.567113  -0.321939  
+-0.567113  2.728827   -0.480842  
+-0.321939  -0.480842  2.378169   
+DEAL::V2:
+5.566198 5.638437 3.202953 
diff --git a/tests/lac/qr_10.cc b/tests/lac/qr_10.cc
new file mode 100644 (file)
index 0000000..aca09f9
--- /dev/null
@@ -0,0 +1,161 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/qr.h>
+
+#include "../tests.h"
+
+/*
+ * Test ImplicitQR::append_column() with lineary dependent columns.
+ * The output should be the same as qr_03 test plus extra output due to
+ * the rejected column.
+ */
+
+/*
+ * MWE in Python for standard QR with 3 linearly independet columns. We will add
+the fourth one:
+
+import numpy as np
+from scipy import linalg
+
+A = np.array([[0, 1, 2], [3, 4, 5], [6, 3, 4.5], [7,9, 8], [10,8,0]])
+Q, R = linalg.qr(A)
+# thin QR
+ncol = A.shape[1]
+R = -R[:ncol,:]
+Q = -Q[:, :ncol]
+np.set_printoptions(precision=6)
+
+print Q
+print R
+
+>>> print A
+[[  0.    1.    2. ]
+ [  3.    4.    5. ]
+ [  6.    3.    4.5]
+ [  7.    9.    8. ]
+ [ 10.    8.    0. ]]
+
+>>> print R
+[[ 13.928388  12.420676   7.03599 ]
+ [ -0.         4.089842   4.916632]
+ [ -0.        -0.         6.290594]]
+>>> print Q
+[[-0.        0.244508  0.126831]
+ [ 0.215387  0.32391   0.300765]
+ [ 0.430775 -0.57472   0.682727]
+ [ 0.502571  0.674288  0.182604]
+ [ 0.717958 -0.224343 -0.627689]]
+
+*/
+
+
+template <typename number>
+void
+test()
+{
+  typedef Vector<number> VectorType;
+  ImplicitQR<VectorType> qr;
+
+  const unsigned int v_size = 5;
+  VectorType         v(v_size);
+  v[0] = 0;
+  v[1] = 3;
+  v[2] = 6;
+  v[3] = 7;
+  v[4] = 10.;
+
+  qr.append_column(v);
+
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+
+  qr.append_column(v);
+
+  v[0] = 2;
+  v[1] = 5;
+  v[2] = 4.5;
+  v[3] = 8;
+  v[4] = 0;
+
+  qr.append_column(v);
+
+  // Now try adding the second column again
+  v[0] = 1;
+  v[1] = 4;
+  v[2] = 3;
+  v[3] = 9;
+  v[4] = 8.;
+
+  const auto fun = [&](const Vector<number> &u,
+                       const number &        rho2,
+                       const number &        col_l2_norm2) {
+    deallog << "Linearly dependent column" << std::endl;
+    deallog << "u:" << std::endl;
+    u.print(deallog.get_file_stream());
+    deallog << "rho2:" << std::endl;
+    deallog.get_file_stream() << rho2 << std::endl;
+    deallog << "col2:" << std::endl;
+    deallog.get_file_stream() << col_l2_norm2 << std::endl;
+    const bool reject =
+      rho2 < dealii::Utilities::fixed_power<2>(1e-4) * col_l2_norm2;
+    deallog << "reject:" << std::endl;
+    deallog.get_file_stream() << reject << std::endl;
+    return !reject;
+  };
+
+  qr.connect_append_column_slot(fun);
+  const bool res = qr.append_column(v);
+  deallog << "Success: " << res << std::endl;
+
+  const LAPACKFullMatrix<double> &R = qr.get_R();
+  std::vector<VectorType>         A(qr.size());
+  for (unsigned int j = 0; j < qr.size(); ++j)
+    {
+      Vector<double> x(qr.size());
+      x    = 0;
+      x[j] = 1.;
+      A[j].reinit(v_size);
+      qr.multiply_with_A(A[j], x);
+    }
+
+  const unsigned int size = qr.size();
+  deallog << size << std::endl;
+  deallog << "R:" << std::endl;
+  R.print_formatted(deallog.get_file_stream(), 6, false, 10);
+  deallog << "A:" << std::endl;
+
+  for (unsigned int i = 0; i < v_size; ++i)
+    for (unsigned int j = 0; j < size; ++j)
+      {
+        deallog.get_file_stream() << std::setw(9) << A[j](i);
+        if (j < size - 1)
+          deallog.get_file_stream() << " ";
+        else
+          deallog.get_file_stream() << std::endl;
+      }
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(6);
+
+  test<double>();
+}
diff --git a/tests/lac/qr_10.output b/tests/lac/qr_10.output
new file mode 100644 (file)
index 0000000..61d4754
--- /dev/null
@@ -0,0 +1,22 @@
+
+DEAL::Linearly dependent column
+DEAL::u:
+1.242e+01 4.090e+00 0.000e+00 
+DEAL::rho2:
+0.00000
+DEAL::col2:
+171.000
+DEAL::reject:
+1
+DEAL::Success: 0
+DEAL::3
+DEAL::R:
+13.928388  12.420676  7.035990   
+           4.089842   4.916632   
+                      6.290594   
+DEAL::A:
+0.00000   1.00000   2.00000  
+3.00000   4.00000   5.00000  
+6.00000   3.00000   4.50000  
+7.00000   9.00000   8.00000  
+10.0000   8.00000   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.