Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
qr.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_qr_h
17 #define dealii_qr_h
18 
19 #include <deal.II/base/exceptions.h>
20 #include <deal.II/base/std_cxx14/memory.h>
21 
22 #include <deal.II/lac/lapack_full_matrix.h>
23 #include <deal.II/lac/lapack_templates.h>
24 #include <deal.II/lac/utilities.h>
25 
26 #include <boost/signals2/signal.hpp>
27 
28 DEAL_II_NAMESPACE_OPEN
29 
41 template <typename VectorType>
42 class BaseQR
43 {
47  typedef typename VectorType::value_type Number;
48 
49 protected:
53  BaseQR();
54 
55 public:
59  virtual ~BaseQR() = default;
60 
67  virtual bool
68  append_column(const VectorType &column) = 0;
69 
73  virtual void
74  remove_column(const unsigned int k = 0) = 0;
75 
79  unsigned int
80  size() const;
81 
86  get_R() const;
87 
93  void
94  solve(Vector<Number> & x,
95  const Vector<Number> &y,
96  const bool transpose = false) const;
97 
102  virtual void
103  multiply_with_Q(VectorType &y, const Vector<Number> &x) const = 0;
104 
109  virtual void
110  multiply_with_QT(Vector<Number> &y, const VectorType &x) const = 0;
111 
116  virtual void
117  multiply_with_A(VectorType &y, const Vector<Number> &x) const = 0;
118 
123  virtual void
124  multiply_with_AT(Vector<Number> &y, const VectorType &x) const = 0;
125 
135  boost::signals2::connection
137  const std::function<void(const unsigned int i,
138  const unsigned int j,
139  const std::array<Number, 3> &csr)> &slot);
140 
141 protected:
146  void
147  multiply_with_cols(VectorType &y, const Vector<Number> &x) const;
148 
152  void
153  multiply_with_colsT(Vector<Number> &y, const VectorType &x) const;
154 
158  std::vector<std::unique_ptr<VectorType>> columns;
159 
164 
168  unsigned int current_size;
169 
174  boost::signals2::signal<void(const unsigned int i,
175  const unsigned int j,
176  const std::array<Number, 3> &)>
178 };
179 
180 // clang-format off
236 // clang-format on
237 template <typename VectorType>
238 class QR : public BaseQR<VectorType>
239 {
240 public:
244  typedef typename VectorType::value_type Number;
245 
249  QR();
250 
254  virtual ~QR() = default;
255 
261  virtual bool
262  append_column(const VectorType &column);
263 
296  virtual void
297  remove_column(const unsigned int k = 0);
298 
299  virtual void
300  multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
301 
302  virtual void
303  multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
304 
305  virtual void
306  multiply_with_A(VectorType &y, const Vector<Number> &x) const;
307 
308  virtual void
309  multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
310 
311 private:
318  void
319  apply_givens_rotation(const unsigned int i, const unsigned int k);
320 
324  VectorType tmp;
325 };
326 
327 
328 
346 template <typename VectorType>
347 class ImplicitQR : public BaseQR<VectorType>
348 {
349 public:
353  typedef typename VectorType::value_type Number;
354 
358  ImplicitQR();
359 
363  virtual ~ImplicitQR() = default;
364 
365  virtual bool
366  append_column(const VectorType &column);
367 
380  virtual void
381  remove_column(const unsigned int k = 0);
382 
383  virtual void
384  multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
385 
386  virtual void
387  multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
388 
389  virtual void
390  multiply_with_A(VectorType &y, const Vector<Number> &x) const;
391 
392  virtual void
393  multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
394 
404  boost::signals2::connection
405  connect_append_column_slot(
406  const std::function<bool(const Vector<Number> &u,
407  const Number & rho2,
408  const Number & col_norm_sqr)> &slot);
409 
410 private:
414  void
415  apply_givens_rotation(const unsigned int i, const unsigned int k);
416 
425  boost::signals2::signal<bool(const Vector<Number> &u,
426  const Number & rho,
427  const Number & col_norm_sqr)>
428  column_signal;
429 };
430 
431 // ------------------- inline and template functions ----------------
432 #ifndef DOXYGEN
433 
434 
435 
436 template <typename VectorType>
438  : current_size(0)
439 {
441 }
442 
443 
444 
445 template <typename VectorType>
446 unsigned int
448 {
449  return current_size;
450 }
451 
452 
453 
454 template <typename VectorType>
457 {
458  return R;
459 }
460 
461 
462 
463 template <typename VectorType>
464 void
465 BaseQR<VectorType>::solve(Vector<Number> & x,
466  const Vector<Number> &y,
467  const bool transpose) const
468 {
469  Assert(x.size() == this->current_size,
470  ExcDimensionMismatch(x.size(), this->current_size));
471  Assert(y.size() == this->current_size,
472  ExcDimensionMismatch(y.size(), this->current_size));
473 
474  // copy if the two vectors are not the same
475  if (&x != &y)
476  x = y;
477 
478  const int lda = this->current_size;
479  const int ldb = this->current_size;
480  const int N = this->current_size;
481  const int n_rhs = 1;
482  int info = 0;
483  trtrs("U",
484  transpose ? "T" : "N",
485  "N",
486  &N,
487  &n_rhs,
488  &this->R(0, 0),
489  &lda,
490  &x(0),
491  &ldb,
492  &info);
493 }
494 
495 
496 
497 template <typename VectorType>
498 void
500  const Vector<Number> &x) const
501 {
502  Assert(x.size() == this->current_size,
503  ExcDimensionMismatch(x.size(), this->current_size));
504 
505  y = 0.;
506  for (unsigned int j = 0; j < this->current_size; ++j)
507  y.add(x[j], *this->columns[j]);
508 }
509 
510 
511 
512 template <typename VectorType>
513 void
514 BaseQR<VectorType>::multiply_with_colsT(Vector<Number> & y,
515  const VectorType &x) const
516 {
517  Assert(y.size() == this->current_size,
518  ExcDimensionMismatch(y.size(), this->current_size));
519 
520  for (unsigned int j = 0; j < this->current_size; ++j)
521  y[j] = (*this->columns[j]) * x;
522 }
523 
524 
525 
526 template <class VectorType>
527 boost::signals2::connection
529  const std::function<void(const unsigned int i,
530  const unsigned int j,
531  const std::array<Number, 3> &)> &slot)
532 {
533  return givens_signal.connect(slot);
534 }
535 
536 
537 
538 template <class VectorType>
539 boost::signals2::connection
541  const std::function<bool(const Vector<Number> &u,
542  const Number & rho,
543  const Number & col_norm_sqr)> &slot)
544 {
545  return column_signal.connect(slot);
546 }
547 
548 
549 
550 template <typename VectorType>
553 {}
554 
555 
556 
557 template <typename VectorType>
558 bool
559 ImplicitQR<VectorType>::append_column(const VectorType &column)
560 {
561  if (this->current_size == 0)
562  {
563  this->R.grow_or_shrink(this->current_size + 1);
564  this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
565  this->R(0, 0) = column.l2_norm();
566  ++this->current_size;
567  }
568  else
569  {
570  // first get scalar products with A^T
571  Vector<Number> u(this->current_size);
572  this->multiply_with_AT(u, column);
573 
574  // now solve R^T x = (A^T * column)
575  const int lda = this->current_size;
576  const int ldb = this->current_size;
577  const int N = this->current_size;
578  const int n_rhs = 1;
579  int info = 0;
580  trtrs(
581  "U", "T", "N", &N, &n_rhs, &this->R(0, 0), &lda, &u(0), &ldb, &info);
582 
583  // finally get the diagonal element:
584  // rho2 = |column|^2 - |u|^2
585  const Number column_norm_sqr = column.norm_sqr();
586  const Number rho2 = column_norm_sqr - u.norm_sqr();
587  const bool linearly_independent =
588  column_signal.empty() ? rho2 > 0 :
589  column_signal(u, rho2, column_norm_sqr).get();
590 
591  // bail out if it turns out to be linearly dependent
592  if (!linearly_independent)
593  return false;
594 
595  // at this point we update is successful and we can enlarge R
596  // and store the column:
597  this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
598  this->R.grow_or_shrink(this->current_size + 1);
599  this->R(this->current_size, this->current_size) = std::sqrt(rho2);
600  for (unsigned int i = 0; i < this->current_size; ++i)
601  this->R(i, this->current_size) = u(i);
602 
603  this->current_size++;
604  }
605 
606  return true;
607 }
608 
609 
610 
611 template <typename VectorType>
612 void
614  const unsigned int k)
615 {
616  Assert(i < k, ExcIndexRange(i, 0, k));
617  Assert(k < this->current_size, ExcIndexRange(k, 0, this->current_size));
618  const std::array<Number, 3> csr =
619  ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
620  this->R(k, k));
621 
622  // first, set k'th column:
623  this->R(i, k) = csr[2];
624  this->R(k, k) = 0.;
625  // now do the rest:
626  for (unsigned int j = 0; j < this->R.n(); ++j)
627  if (j != k)
628  {
629  const Number t = this->R(i, j);
630  this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
631  this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
632  }
633 
634  if (!this->givens_signal.empty())
635  this->givens_signal(i, k, csr);
636 }
637 
638 
639 
640 template <typename VectorType>
641 void
642 ImplicitQR<VectorType>::remove_column(const unsigned int k)
643 {
644  // before actually removing a column from Q and resizing R,
645  // apply givens rotations to bring H into upper triangular form:
646  for (unsigned int j = k + 1; j < this->R.n(); ++j)
647  {
648  const unsigned int i = j - 1;
649  apply_givens_rotation(i, j);
650  }
651 
652  // remove last row and k-th column
653  --this->current_size;
654  this->R.remove_row_and_column(this->current_size, k);
655 
656  // Finally remove the column from A
657  this->columns.erase(this->columns.begin() + k);
658 }
659 
660 
661 
662 template <typename VectorType>
663 void
665  const Vector<Number> &x) const
666 {
667  // A = QR
668  // A R^{-1} = Q
669  Vector<Number> x1 = x;
670  BaseQR<VectorType>::solve(x1, x1, false);
671  multiply_with_A(y, x1);
672 }
673 
674 
675 
676 template <typename VectorType>
677 void
679  const VectorType &x) const
680 {
681  // A = QR
682  // A^T = R^T Q^T
683  // {R^T}^{-1} A^T = Q^T
684  multiply_with_AT(y, x);
685  BaseQR<VectorType>::solve(y, y, true);
686 }
687 
688 
689 
690 template <typename VectorType>
691 void
693  const Vector<Number> &x) const
694 {
696 }
697 
698 
699 
700 template <typename VectorType>
701 void
703  const VectorType &x) const
704 {
706 }
707 
708 
709 
710 template <typename VectorType>
713 {}
714 
715 
716 
717 template <typename VectorType>
718 bool
719 QR<VectorType>::append_column(const VectorType &column)
720 {
721  // resize R:
722  this->R.grow_or_shrink(this->current_size + 1);
723  this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
724 
725  // now a Gram-Schmidt part: orthonormalize the new column
726  // against everything we have so far:
727  auto &last_col = *this->columns.back();
728  for (unsigned int i = 0; i < this->current_size; ++i)
729  {
730  const auto &i_col = *this->columns[i];
731  this->R(i, this->current_size) = i_col * last_col;
732  last_col.add(-this->R(i, this->current_size), i_col);
733  }
734 
735  this->R(this->current_size, this->current_size) = last_col.l2_norm();
736 
737  Assert(this->R(this->current_size, this->current_size) > 0.,
738  ExcDivideByZero());
739  last_col *= 1. / this->R(this->current_size, this->current_size);
740 
741  ++this->current_size;
742  return true;
743 }
744 
745 
746 
747 template <typename VectorType>
748 void
749 QR<VectorType>::apply_givens_rotation(const unsigned int i,
750  const unsigned int k)
751 {
752  Assert(i < k, ExcIndexRange(i, 0, k));
753  Assert(k < this->current_size, ExcIndexRange(k, 0, this->current_size));
754  const std::array<Number, 3> csr =
755  ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
756  this->R(k, k));
757 
758  // first, set k'th column:
759  this->R(i, k) = csr[2];
760  this->R(k, k) = 0.;
761  // now do the rest:
762  for (unsigned int j = 0; j < this->R.n(); ++j)
763  if (j != k)
764  {
765  const Number t = this->R(i, j);
766  this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
767  this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
768  }
769 
770  // now adjust i,k columns due to multiplication with the
771  // transpose Givens matrix from right:
772  auto &col_i = *this->columns[i];
773  auto &col_k = *this->columns[k];
774  // save column i:
775  tmp = col_i;
776  col_i.sadd(csr[0], csr[1], col_k);
777  col_k.sadd(csr[0], -csr[1], tmp);
778 
779  if (!this->givens_signal.empty())
780  this->givens_signal(i, k, csr);
781 }
782 
783 
784 
785 template <typename VectorType>
786 void
787 QR<VectorType>::remove_column(const unsigned int k)
788 {
789  Assert(k < this->current_size, ExcIndexRange(k, 0, this->current_size));
790  Assert(this->current_size > 0,
791  ExcMessage("Can not remove a column if QR is empty"));
792  // apply a sequence of Givens rotations
793  // see section 6.5 "Updating matrix factorizations" in Golub 2013, Matrix
794  // computations
795 
796  // So we want to have QR for \tilde A \in R^{m*(n-1)}
797  // if we remove the column k, we end up with upper Hessenberg matrix
798  // x x x x x
799  // x x x x
800  // H = x x x
801  // x x x
802  // x x
803  // x
804  // where k = 2 (3rd column), m = 7, n = 6
805  //
806  // before actually removing a column from Q and resizing R,
807  // apply givens rotations to bring H into upper triangular form:
808  for (unsigned int j = k + 1; j < this->R.n(); ++j)
809  {
810  const unsigned int i = j - 1;
811  apply_givens_rotation(i, j);
812  }
813 
814  // now we can throw away the column from Q and adjust R
815  // since we do thin-QR, after Givens rotations we need to throw
816  // away the last column:
817  const unsigned int size_minus_1 = this->columns.size() - 1;
818  this->columns.erase(this->columns.begin() + size_minus_1);
819 
820  // remove last row and k-th column
821  --this->current_size;
822  this->R.remove_row_and_column(this->current_size, k);
823 }
824 
825 
826 
827 template <typename VectorType>
828 void
829 QR<VectorType>::multiply_with_Q(VectorType &y, const Vector<Number> &x) const
830 {
832 }
833 
834 
835 
836 template <typename VectorType>
837 void
838 QR<VectorType>::multiply_with_QT(Vector<Number> &y, const VectorType &x) const
839 {
841 }
842 
843 
844 
845 template <typename VectorType>
846 void
847 QR<VectorType>::multiply_with_A(VectorType &y, const Vector<Number> &x) const
848 {
849  Vector<Number> x1 = x;
850  const int N = this->current_size;
851  const int lda = N;
852  const int incx = 1;
853  trmv("U", "N", "N", &N, &this->R(0, 0), &lda, &x1[0], &incx);
854 
855  multiply_with_Q(y, x1);
856 }
857 
858 
859 
860 template <typename VectorType>
861 void
862 QR<VectorType>::multiply_with_AT(Vector<Number> &y, const VectorType &x) const
863 {
864  multiply_with_QT(y, x);
865 
866  const int N = this->current_size;
867  const int lda = N;
868  const int incx = 1;
869  trmv("U", "T", "N", &N, &this->R(0, 0), &lda, &y[0], &incx);
870 }
871 
872 #endif // no DOXYGEN
873 
874 DEAL_II_NAMESPACE_CLOSE
875 
876 #endif
unsigned int size() const
QR()
VectorType::value_type Number
Definition: qr.h:47
std::vector< std::unique_ptr< VectorType > > columns
Definition: qr.h:158
void remove_row_and_column(const size_type row, const size_type col)
Matrix is upper triangular.
void multiply_with_colsT(Vector< Number > &y, const VectorType &x) const
size_type n() const
virtual bool append_column(const VectorType &column)
virtual ~BaseQR()=default
virtual void multiply_with_A(VectorType &y, const Vector< Number > &x) const =0
unsigned int current_size
Definition: qr.h:168
virtual void multiply_with_AT(Vector< Number > &y, const VectorType &x) const
virtual void multiply_with_Q(VectorType &y, const Vector< Number > &x) const =0
virtual bool append_column(const VectorType &column)=0
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual void multiply_with_Q(VectorType &y, const Vector< Number > &x) const
static ::ExceptionBase & ExcDivideByZero()
LAPACKFullMatrix< Number > R
Definition: qr.h:163
virtual void multiply_with_A(VectorType &y, const Vector< Number > &x) const
boost::signals2::signal< bool(const Vector< Number > &u, const Number &rho, const Number &col_norm_sqr)> column_signal
Definition: qr.h:428
boost::signals2::connection connect_append_column_slot(const std::function< bool(const Vector< Number > &u, const Number &rho2, const Number &col_norm_sqr)> &slot)
static ::ExceptionBase & ExcMessage(std::string arg1)
void solve(Vector< Number > &x, const Vector< Number > &y, const bool transpose=false) const
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)
VectorType tmp
Definition: qr.h:324
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
VectorType::value_type Number
Definition: qr.h:353
virtual void multiply_with_AT(Vector< Number > &y, const VectorType &x) const =0
virtual void remove_column(const unsigned int k=0)
void add(const number a, const LAPACKFullMatrix< number > &B)
boost::signals2::signal< void(const unsigned int i, const unsigned int j, const std::array< Number, 3 > &)> givens_signal
Definition: qr.h:177
virtual void multiply_with_A(VectorType &y, const Vector< Number > &x) const
virtual void multiply_with_AT(Vector< Number > &y, const VectorType &x) const
const LAPACKFullMatrix< Number > & get_R() const
virtual void multiply_with_Q(VectorType &y, const Vector< Number > &x) const
void apply_givens_rotation(const unsigned int i, const unsigned int k)
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const
void grow_or_shrink(const size_type size)
Definition: qr.h:347
virtual bool append_column(const VectorType &column)
virtual void remove_column(const unsigned int k=0)
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const
virtual void remove_column(const unsigned int k=0)=0
VectorType::value_type Number
Definition: qr.h:244
Definition: qr.h:238
void apply_givens_rotation(const unsigned int i, const unsigned int k)
Definition: qr.h:42
void set_property(const LAPACKSupport::Property property)
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const =0
void multiply_with_cols(VectorType &y, const Vector< Number > &x) const