From 22e4f21c8df16ca9ed38b2751f7978d01ac4f71c Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Fri, 21 Sep 2018 14:47:06 +0200 Subject: [PATCH] add QR classes --- include/deal.II/lac/qr.h | 876 +++++++++++++++++++++++++++++++++++++++ tests/lac/qr.cc | 130 ++++++ tests/lac/qr.output | 12 + tests/lac/qr_02.cc | 157 +++++++ tests/lac/qr_02.output | 24 ++ tests/lac/qr_03.cc | 132 ++++++ tests/lac/qr_03.output | 12 + tests/lac/qr_04.cc | 158 +++++++ tests/lac/qr_04.output | 24 ++ tests/lac/qr_05.cc | 113 +++++ tests/lac/qr_05.output | 8 + tests/lac/qr_06.cc | 129 ++++++ tests/lac/qr_06.output | 10 + tests/lac/qr_07.cc | 130 ++++++ tests/lac/qr_07.output | 10 + tests/lac/qr_08.cc | 287 +++++++++++++ tests/lac/qr_08.output | 29 ++ tests/lac/qr_09.cc | 276 ++++++++++++ tests/lac/qr_09.output | 22 + tests/lac/qr_10.cc | 161 +++++++ tests/lac/qr_10.output | 22 + 21 files changed, 2722 insertions(+) create mode 100644 include/deal.II/lac/qr.h create mode 100644 tests/lac/qr.cc create mode 100644 tests/lac/qr.output create mode 100644 tests/lac/qr_02.cc create mode 100644 tests/lac/qr_02.output create mode 100644 tests/lac/qr_03.cc create mode 100644 tests/lac/qr_03.output create mode 100644 tests/lac/qr_04.cc create mode 100644 tests/lac/qr_04.output create mode 100644 tests/lac/qr_05.cc create mode 100644 tests/lac/qr_05.output create mode 100644 tests/lac/qr_06.cc create mode 100644 tests/lac/qr_06.output create mode 100644 tests/lac/qr_07.cc create mode 100644 tests/lac/qr_07.output create mode 100644 tests/lac/qr_08.cc create mode 100644 tests/lac/qr_08.output create mode 100644 tests/lac/qr_09.cc create mode 100644 tests/lac/qr_09.output create mode 100644 tests/lac/qr_10.cc create mode 100644 tests/lac/qr_10.output diff --git a/include/deal.II/lac/qr.h b/include/deal.II/lac/qr.h new file mode 100644 index 0000000000..1898d6fb14 --- /dev/null +++ b/include/deal.II/lac/qr.h @@ -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 +#include + +#include +#include +#include + +#include + +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 +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 true if the result is successful, i.e. + * the columns are linearly independent. Otherwise the @p column + * is rejected and the return value is false. + */ + 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 & + 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 true, $R^Tx=y$ is solved instead. + */ + void + solve(Vector & x, + const Vector &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 &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 &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 &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 &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 &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 &x) const; + + /** + * Multiply with transpose columns stored in the object. + */ + void + multiply_with_colsT(Vector &y, const VectorType &x) const; + + /** + * A vector of unique pointers to store columns. + */ + std::vector> columns; + + /** + * Matrix to store R. + */ + LAPACKFullMatrix 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 &)> + 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 +class QR : public BaseQR +{ +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 true. + */ + 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 &x) const; + + virtual void + multiply_with_QT(Vector &y, const VectorType &x) const; + + virtual void + multiply_with_A(VectorType &y, const Vector &x) const; + + virtual void + multiply_with_AT(Vector &y, const VectorType &x) const; + +private: + /** + * Apply givens rotation in plane @p i @p k to @p Q and @p R so that + * R(k,k) 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 +class ImplicitQR : public BaseQR +{ +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 &x) const; + + virtual void + multiply_with_QT(Vector &y, const VectorType &x) const; + + virtual void + multiply_with_A(VectorType &y, const Vector &x) const; + + virtual void + multiply_with_AT(Vector &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 true if the new column is + * linearly independent. + */ + boost::signals2::connection + connect_append_column_slot( + const std::function &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 true if the new column is + * linearly independent. + */ + boost::signals2::signal &u, + const Number & rho, + const Number & col_norm_sqr)> + column_signal; +}; + +// ------------------- inline and template functions ---------------- +#ifndef DOXYGEN + + + +template +BaseQR::BaseQR() + : current_size(0) +{ + R.set_property(LAPACKSupport::upper_triangular); +} + + + +template +unsigned int +BaseQR::size() const +{ + return current_size; +} + + + +template +const LAPACKFullMatrix::Number> & +BaseQR::get_R() const +{ + return R; +} + + + +template +void +BaseQR::solve(Vector & x, + const Vector &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 +void +BaseQR::multiply_with_cols(VectorType & y, + const Vector &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 +void +BaseQR::multiply_with_colsT(Vector & 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 +boost::signals2::connection +BaseQR::connect_givens_slot( + const std::function &)> &slot) +{ + return givens_signal.connect(slot); +} + + + +template +boost::signals2::connection +ImplicitQR::connect_append_column_slot( + const std::function &u, + const Number & rho, + const Number & col_norm_sqr)> &slot) +{ + return column_signal.connect(slot); +} + + + +template +ImplicitQR::ImplicitQR() + : BaseQR() +{} + + + +template +bool +ImplicitQR::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(column)); + this->R(0, 0) = column.l2_norm(); + ++this->current_size; + } + else + { + // first get scalar products with A^T + Vector 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(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 +void +ImplicitQR::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 csr = + dealii::Utilities::LinearAlgebra::givens_rotation(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 +void +ImplicitQR::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 +void +ImplicitQR::multiply_with_Q(VectorType & y, + const Vector &x) const +{ + // A = QR + // A R^{-1} = Q + Vector x1 = x; + BaseQR::solve(x1, x1, false); + multiply_with_A(y, x1); +} + + + +template +void +ImplicitQR::multiply_with_QT(Vector & 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::solve(y, y, true); +} + + + +template +void +ImplicitQR::multiply_with_A(VectorType & y, + const Vector &x) const +{ + BaseQR::multiply_with_cols(y, x); +} + + + +template +void +ImplicitQR::multiply_with_AT(Vector & y, + const VectorType &x) const +{ + BaseQR::multiply_with_colsT(y, x); +} + + + +template +QR::QR() + : BaseQR() +{} + + + +template +bool +QR::append_column(const VectorType &column) +{ + // resize R: + this->R.grow_or_shrink(this->current_size + 1); + this->columns.push_back(std_cxx14::make_unique(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 +void +QR::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 csr = + dealii::Utilities::LinearAlgebra::givens_rotation(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 +void +QR::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 +void +QR::multiply_with_Q(VectorType &y, const Vector &x) const +{ + BaseQR::multiply_with_cols(y, x); +} + + + +template +void +QR::multiply_with_QT(Vector &y, const VectorType &x) const +{ + BaseQR::multiply_with_colsT(y, x); +} + + + +template +void +QR::multiply_with_A(VectorType &y, const Vector &x) const +{ + Vector 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 +void +QR::multiply_with_AT(Vector &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 index 0000000000..e93e676e50 --- /dev/null +++ b/tests/lac/qr.cc @@ -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 + +#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 +void +test() +{ + typedef Vector VectorType; + QR 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 &R = qr.get_R(); + std::vector Q(qr.size()); + for (unsigned int j = 0; j < qr.size(); ++j) + { + Vector 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(); +} diff --git a/tests/lac/qr.output b/tests/lac/qr.output new file mode 100644 index 0000000000..77e4fa6d72 --- /dev/null +++ b/tests/lac/qr.output @@ -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 index 0000000000..e5bb3cc13d --- /dev/null +++ b/tests/lac/qr_02.cc @@ -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 + +#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 +void +print(const QR &qr, const unsigned int v_size) +{ + const LAPACKFullMatrix &R = qr.get_R(); + std::vector Q(qr.size()); + for (unsigned int j = 0; j < qr.size(); ++j) + { + Vector 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 +void +test() +{ + typedef Vector VectorType; + QR qr; + + auto print_givens = [](const unsigned int i, + const unsigned int j, + const std::array &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(); +} diff --git a/tests/lac/qr_02.output b/tests/lac/qr_02.output new file mode 100644 index 0000000000..86f3738e08 --- /dev/null +++ b/tests/lac/qr_02.output @@ -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 index 0000000000..5777677b92 --- /dev/null +++ b/tests/lac/qr_03.cc @@ -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 + +#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 +void +test() +{ + typedef Vector VectorType; + ImplicitQR 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 &R = qr.get_R(); + std::vector Q(R.n()); + for (unsigned int j = 0; j < R.n(); ++j) + { + Vector 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(); +} diff --git a/tests/lac/qr_03.output b/tests/lac/qr_03.output new file mode 100644 index 0000000000..eec779fe53 --- /dev/null +++ b/tests/lac/qr_03.output @@ -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 index 0000000000..f3e7815d71 --- /dev/null +++ b/tests/lac/qr_04.cc @@ -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 + +#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 +void +print(const ImplicitQR &qr, const unsigned int col_size) +{ + const LAPACKFullMatrix &R = qr.get_R(); + std::vector A(R.n()); + for (unsigned int j = 0; j < R.n(); ++j) + { + Vector 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 +void +test() +{ + typedef Vector VectorType; + ImplicitQR qr; + + auto print_givens = [](const unsigned int i, + const unsigned int j, + const std::array &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(); +} diff --git a/tests/lac/qr_04.output b/tests/lac/qr_04.output new file mode 100644 index 0000000000..ab595f16f5 --- /dev/null +++ b/tests/lac/qr_04.output @@ -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 index 0000000000..5d52ba1ed7 --- /dev/null +++ b/tests/lac/qr_05.cc @@ -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 + +#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 +void +test() +{ + typedef Vector VectorType; + QR 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 &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 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(); +} diff --git a/tests/lac/qr_05.output b/tests/lac/qr_05.output new file mode 100644 index 0000000000..fd177ab6cc --- /dev/null +++ b/tests/lac/qr_05.output @@ -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 index 0000000000..89c95ac92d --- /dev/null +++ b/tests/lac/qr_06.cc @@ -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 + +#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 +void +test() +{ + typedef Vector VectorType; + QR 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 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(); +} diff --git a/tests/lac/qr_06.output b/tests/lac/qr_06.output new file mode 100644 index 0000000000..ee397c8d7f --- /dev/null +++ b/tests/lac/qr_06.output @@ -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 index 0000000000..c9c24bdc3d --- /dev/null +++ b/tests/lac/qr_07.cc @@ -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 + +#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 +void +test() +{ + typedef Vector VectorType; + ImplicitQR 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 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(); +} diff --git a/tests/lac/qr_07.output b/tests/lac/qr_07.output new file mode 100644 index 0000000000..ee397c8d7f --- /dev/null +++ b/tests/lac/qr_07.output @@ -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 index 0000000000..6b88b09a79 --- /dev/null +++ b/tests/lac/qr_08.cc @@ -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 +#include + +#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 +void +print(const QR &qr, const unsigned int col_size) +{ + const LAPACKFullMatrix &R = qr.get_R(); + std::vector Q(R.n()); + for (unsigned int j = 0; j < R.n(); ++j) + { + Vector 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 +void +test() +{ + typedef Vector VectorType; + QR 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 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 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 Q(qr.size()); + for (unsigned int j = 0; j < qr.size(); ++j) + { + Vector x(qr.size()); + x = 0; + x[j] = 1.; + Q[j].reinit(v_size); + qr.multiply_with_Q(Q[j], x); + } + + LAPACKFullMatrix 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 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 &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 V2(3); + LAPACKFullMatrix 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(); +} diff --git a/tests/lac/qr_08.output b/tests/lac/qr_08.output new file mode 100644 index 0000000000..674021249a --- /dev/null +++ b/tests/lac/qr_08.output @@ -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 index 0000000000..3d90c9fa2d --- /dev/null +++ b/tests/lac/qr_09.cc @@ -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 +#include + +#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 +void +print(const ImplicitQR &qr) +{ + const LAPACKFullMatrix &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 +void +test() +{ + typedef Vector VectorType; + QR qr; + ImplicitQR 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 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 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 Q(qr.size()); + for (unsigned int j = 0; j < qr.size(); ++j) + { + Vector x(qr.size()); + x = 0; + x[j] = 1.; + Q[j].reinit(v_size); + qr.multiply_with_Q(Q[j], x); + } + + LAPACKFullMatrix 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 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 &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 V2(3); + LAPACKFullMatrix 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(); +} diff --git a/tests/lac/qr_09.output b/tests/lac/qr_09.output new file mode 100644 index 0000000000..64b93fc783 --- /dev/null +++ b/tests/lac/qr_09.output @@ -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 index 0000000000..aca09f967f --- /dev/null +++ b/tests/lac/qr_10.cc @@ -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 + +#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 +void +test() +{ + typedef Vector VectorType; + ImplicitQR 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 &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 &R = qr.get_R(); + std::vector A(qr.size()); + for (unsigned int j = 0; j < qr.size(); ++j) + { + Vector 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(); +} diff --git a/tests/lac/qr_10.output b/tests/lac/qr_10.output new file mode 100644 index 0000000000..61d47547b3 --- /dev/null +++ b/tests/lac/qr_10.output @@ -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 -- 2.39.5