--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#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>();
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff