Reference documentation for deal.II version Git b84270a1d4 2020-01-17 16:21:02 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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/config.h>
20 
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/std_cxx14/memory.h>
23 
24 #include <deal.II/lac/lapack_full_matrix.h>
25 #include <deal.II/lac/lapack_templates.h>
26 #include <deal.II/lac/utilities.h>
27 
28 #include <boost/signals2/signal.hpp>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
43 template <typename VectorType>
44 class BaseQR
45 {
49  typedef typename VectorType::value_type Number;
50 
51 protected:
55  BaseQR();
56 
57 public:
61  virtual ~BaseQR() = default;
62 
69  virtual bool
70  append_column(const VectorType &column) = 0;
71 
75  virtual void
76  remove_column(const unsigned int k = 0) = 0;
77 
81  unsigned int
82  size() const;
83 
88  get_R() const;
89 
95  void
96  solve(Vector<Number> & x,
97  const Vector<Number> &y,
98  const bool transpose = false) const;
99 
104  virtual void
105  multiply_with_Q(VectorType &y, const Vector<Number> &x) const = 0;
106 
111  virtual void
112  multiply_with_QT(Vector<Number> &y, const VectorType &x) const = 0;
113 
118  virtual void
119  multiply_with_A(VectorType &y, const Vector<Number> &x) const = 0;
120 
125  virtual void
126  multiply_with_AT(Vector<Number> &y, const VectorType &x) const = 0;
127 
137  boost::signals2::connection
139  const std::function<void(const unsigned int i,
140  const unsigned int j,
141  const std::array<Number, 3> &csr)> &slot);
142 
143 protected:
148  void
149  multiply_with_cols(VectorType &y, const Vector<Number> &x) const;
150 
154  void
155  multiply_with_colsT(Vector<Number> &y, const VectorType &x) const;
156 
160  std::vector<std::unique_ptr<VectorType>> columns;
161 
166 
170  unsigned int current_size;
171 
176  boost::signals2::signal<void(const unsigned int i,
177  const unsigned int j,
178  const std::array<Number, 3> &)>
180 };
181 
182 // clang-format off
238 // clang-format on
239 template <typename VectorType>
240 class QR : public BaseQR<VectorType>
241 {
242 public:
246  typedef typename VectorType::value_type Number;
247 
251  QR();
252 
256  virtual ~QR() = default;
257 
263  virtual bool
264  append_column(const VectorType &column);
265 
298  virtual void
299  remove_column(const unsigned int k = 0);
300 
301  virtual void
302  multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
303 
304  virtual void
305  multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
306 
307  virtual void
308  multiply_with_A(VectorType &y, const Vector<Number> &x) const;
309 
310  virtual void
311  multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
312 
313 private:
320  void
321  apply_givens_rotation(const unsigned int i, const unsigned int k);
322 
326  VectorType tmp;
327 };
328 
329 
330 
348 template <typename VectorType>
349 class ImplicitQR : public BaseQR<VectorType>
350 {
351 public:
355  typedef typename VectorType::value_type Number;
356 
360  ImplicitQR();
361 
365  virtual ~ImplicitQR() = default;
366 
367  virtual bool
368  append_column(const VectorType &column);
369 
382  virtual void
383  remove_column(const unsigned int k = 0);
384 
385  virtual void
386  multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
387 
388  virtual void
389  multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
390 
391  virtual void
392  multiply_with_A(VectorType &y, const Vector<Number> &x) const;
393 
394  virtual void
395  multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
396 
406  boost::signals2::connection
407  connect_append_column_slot(
408  const std::function<bool(const Vector<Number> &u,
409  const Number & rho2,
410  const Number & col_norm_sqr)> &slot);
411 
412 private:
416  void
417  apply_givens_rotation(const unsigned int i, const unsigned int k);
418 
427  boost::signals2::signal<bool(const Vector<Number> &u,
428  const Number & rho,
429  const Number & col_norm_sqr)>
430  column_signal;
431 };
432 
433 // ------------------- inline and template functions ----------------
434 #ifndef DOXYGEN
435 
436 
437 
438 template <typename VectorType>
440  : current_size(0)
441 {
443 }
444 
445 
446 
447 template <typename VectorType>
448 unsigned int
450 {
451  return current_size;
452 }
453 
454 
455 
456 template <typename VectorType>
459 {
460  return R;
461 }
462 
463 
464 
465 template <typename VectorType>
466 void
467 BaseQR<VectorType>::solve(Vector<Number> & x,
468  const Vector<Number> &y,
469  const bool transpose) const
470 {
471  Assert(x.size() == this->current_size,
472  ExcDimensionMismatch(x.size(), this->current_size));
473  Assert(y.size() == this->current_size,
474  ExcDimensionMismatch(y.size(), this->current_size));
475 
476  // copy if the two vectors are not the same
477  if (&x != &y)
478  x = y;
479 
480  const int lda = this->current_size;
481  const int ldb = this->current_size;
482  const int N = this->current_size;
483  const int n_rhs = 1;
484  int info = 0;
485  trtrs("U",
486  transpose ? "T" : "N",
487  "N",
488  &N,
489  &n_rhs,
490  &this->R(0, 0),
491  &lda,
492  &x(0),
493  &ldb,
494  &info);
495 }
496 
497 
498 
499 template <typename VectorType>
500 void
502  const Vector<Number> &x) const
503 {
504  Assert(x.size() == this->current_size,
505  ExcDimensionMismatch(x.size(), this->current_size));
506 
507  y = 0.;
508  for (unsigned int j = 0; j < this->current_size; ++j)
509  y.add(x[j], *this->columns[j]);
510 }
511 
512 
513 
514 template <typename VectorType>
515 void
516 BaseQR<VectorType>::multiply_with_colsT(Vector<Number> & y,
517  const VectorType &x) const
518 {
519  Assert(y.size() == this->current_size,
520  ExcDimensionMismatch(y.size(), this->current_size));
521 
522  for (unsigned int j = 0; j < this->current_size; ++j)
523  y[j] = (*this->columns[j]) * x;
524 }
525 
526 
527 
528 template <class VectorType>
529 boost::signals2::connection
531  const std::function<void(const unsigned int i,
532  const unsigned int j,
533  const std::array<Number, 3> &)> &slot)
534 {
535  return givens_signal.connect(slot);
536 }
537 
538 
539 
540 template <class VectorType>
541 boost::signals2::connection
543  const std::function<bool(const Vector<Number> &u,
544  const Number & rho,
545  const Number & col_norm_sqr)> &slot)
546 {
547  return column_signal.connect(slot);
548 }
549 
550 
551 
552 template <typename VectorType>
555 {}
556 
557 
558 
559 template <typename VectorType>
560 bool
561 ImplicitQR<VectorType>::append_column(const VectorType &column)
562 {
563  if (this->current_size == 0)
564  {
565  this->R.grow_or_shrink(this->current_size + 1);
566  this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
567  this->R(0, 0) = column.l2_norm();
568  ++this->current_size;
569  }
570  else
571  {
572  // first get scalar products with A^T
573  Vector<Number> u(this->current_size);
574  this->multiply_with_AT(u, column);
575 
576  // now solve R^T x = (A^T * column)
577  const int lda = this->current_size;
578  const int ldb = this->current_size;
579  const int N = this->current_size;
580  const int n_rhs = 1;
581  int info = 0;
582  trtrs(
583  "U", "T", "N", &N, &n_rhs, &this->R(0, 0), &lda, &u(0), &ldb, &info);
584 
585  // finally get the diagonal element:
586  // rho2 = |column|^2 - |u|^2
587  const Number column_norm_sqr = column.norm_sqr();
588  const Number rho2 = column_norm_sqr - u.norm_sqr();
589  const bool linearly_independent =
590  column_signal.empty() ? rho2 > 0 :
591  column_signal(u, rho2, column_norm_sqr).get();
592 
593  // bail out if it turns out to be linearly dependent
594  if (!linearly_independent)
595  return false;
596 
597  // at this point we update is successful and we can enlarge R
598  // and store the column:
599  this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
600  this->R.grow_or_shrink(this->current_size + 1);
601  this->R(this->current_size, this->current_size) = std::sqrt(rho2);
602  for (unsigned int i = 0; i < this->current_size; ++i)
603  this->R(i, this->current_size) = u(i);
604 
605  this->current_size++;
606  }
607 
608  return true;
609 }
610 
611 
612 
613 template <typename VectorType>
614 void
616  const unsigned int k)
617 {
618  AssertIndexRange(i, k);
619  AssertIndexRange(k, this->current_size);
620  const std::array<Number, 3> csr =
621  ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
622  this->R(k, k));
623 
624  // first, set k'th column:
625  this->R(i, k) = csr[2];
626  this->R(k, k) = 0.;
627  // now do the rest:
628  for (unsigned int j = 0; j < this->R.n(); ++j)
629  if (j != k)
630  {
631  const Number t = this->R(i, j);
632  this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
633  this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
634  }
635 
636  if (!this->givens_signal.empty())
637  this->givens_signal(i, k, csr);
638 }
639 
640 
641 
642 template <typename VectorType>
643 void
644 ImplicitQR<VectorType>::remove_column(const unsigned int k)
645 {
646  // before actually removing a column from Q and resizing R,
647  // apply givens rotations to bring H into upper triangular form:
648  for (unsigned int j = k + 1; j < this->R.n(); ++j)
649  {
650  const unsigned int i = j - 1;
651  apply_givens_rotation(i, j);
652  }
653 
654  // remove last row and k-th column
655  --this->current_size;
656  this->R.remove_row_and_column(this->current_size, k);
657 
658  // Finally remove the column from A
659  this->columns.erase(this->columns.begin() + k);
660 }
661 
662 
663 
664 template <typename VectorType>
665 void
667  const Vector<Number> &x) const
668 {
669  // A = QR
670  // A R^{-1} = Q
671  Vector<Number> x1 = x;
672  BaseQR<VectorType>::solve(x1, x1, false);
673  multiply_with_A(y, x1);
674 }
675 
676 
677 
678 template <typename VectorType>
679 void
681  const VectorType &x) const
682 {
683  // A = QR
684  // A^T = R^T Q^T
685  // {R^T}^{-1} A^T = Q^T
686  multiply_with_AT(y, x);
687  BaseQR<VectorType>::solve(y, y, true);
688 }
689 
690 
691 
692 template <typename VectorType>
693 void
695  const Vector<Number> &x) const
696 {
698 }
699 
700 
701 
702 template <typename VectorType>
703 void
705  const VectorType &x) const
706 {
708 }
709 
710 
711 
712 template <typename VectorType>
715 {}
716 
717 
718 
719 template <typename VectorType>
720 bool
721 QR<VectorType>::append_column(const VectorType &column)
722 {
723  // resize R:
724  this->R.grow_or_shrink(this->current_size + 1);
725  this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
726 
727  // now a Gram-Schmidt part: orthonormalize the new column
728  // against everything we have so far:
729  auto &last_col = *this->columns.back();
730  for (unsigned int i = 0; i < this->current_size; ++i)
731  {
732  const auto &i_col = *this->columns[i];
733  this->R(i, this->current_size) = i_col * last_col;
734  last_col.add(-this->R(i, this->current_size), i_col);
735  }
736 
737  this->R(this->current_size, this->current_size) = last_col.l2_norm();
738 
739  Assert(this->R(this->current_size, this->current_size) > 0.,
740  ExcDivideByZero());
741  last_col *= 1. / this->R(this->current_size, this->current_size);
742 
743  ++this->current_size;
744  return true;
745 }
746 
747 
748 
749 template <typename VectorType>
750 void
751 QR<VectorType>::apply_givens_rotation(const unsigned int i,
752  const unsigned int k)
753 {
754  AssertIndexRange(i, k);
755  AssertIndexRange(k, this->current_size);
756  const std::array<Number, 3> csr =
757  ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
758  this->R(k, k));
759 
760  // first, set k'th column:
761  this->R(i, k) = csr[2];
762  this->R(k, k) = 0.;
763  // now do the rest:
764  for (unsigned int j = 0; j < this->R.n(); ++j)
765  if (j != k)
766  {
767  const Number t = this->R(i, j);
768  this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
769  this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
770  }
771 
772  // now adjust i,k columns due to multiplication with the
773  // transpose Givens matrix from right:
774  auto &col_i = *this->columns[i];
775  auto &col_k = *this->columns[k];
776  // save column i:
777  tmp = col_i;
778  col_i.sadd(csr[0], csr[1], col_k);
779  col_k.sadd(csr[0], -csr[1], tmp);
780 
781  if (!this->givens_signal.empty())
782  this->givens_signal(i, k, csr);
783 }
784 
785 
786 
787 template <typename VectorType>
788 void
789 QR<VectorType>::remove_column(const unsigned int k)
790 {
791  AssertIndexRange(k, this->current_size);
792  Assert(this->current_size > 0,
793  ExcMessage("Can not remove a column if QR is empty"));
794  // apply a sequence of Givens rotations
795  // see section 6.5 "Updating matrix factorizations" in Golub 2013, Matrix
796  // computations
797 
798  // So we want to have QR for \tilde A \in R^{m*(n-1)}
799  // if we remove the column k, we end up with upper Hessenberg matrix
800  // x x x x x
801  // x x x x
802  // H = x x x
803  // x x x
804  // x x
805  // x
806  // where k = 2 (3rd column), m = 7, n = 6
807  //
808  // before actually removing a column from Q and resizing R,
809  // apply givens rotations to bring H into upper triangular form:
810  for (unsigned int j = k + 1; j < this->R.n(); ++j)
811  {
812  const unsigned int i = j - 1;
813  apply_givens_rotation(i, j);
814  }
815 
816  // now we can throw away the column from Q and adjust R
817  // since we do thin-QR, after Givens rotations we need to throw
818  // away the last column:
819  const unsigned int size_minus_1 = this->columns.size() - 1;
820  this->columns.erase(this->columns.begin() + size_minus_1);
821 
822  // remove last row and k-th column
823  --this->current_size;
824  this->R.remove_row_and_column(this->current_size, k);
825 }
826 
827 
828 
829 template <typename VectorType>
830 void
831 QR<VectorType>::multiply_with_Q(VectorType &y, const Vector<Number> &x) const
832 {
834 }
835 
836 
837 
838 template <typename VectorType>
839 void
840 QR<VectorType>::multiply_with_QT(Vector<Number> &y, const VectorType &x) const
841 {
843 }
844 
845 
846 
847 template <typename VectorType>
848 void
849 QR<VectorType>::multiply_with_A(VectorType &y, const Vector<Number> &x) const
850 {
851  Vector<Number> x1 = x;
852  const int N = this->current_size;
853  const int lda = N;
854  const int incx = 1;
855  trmv("U", "N", "N", &N, &this->R(0, 0), &lda, &x1[0], &incx);
856 
857  multiply_with_Q(y, x1);
858 }
859 
860 
861 
862 template <typename VectorType>
863 void
864 QR<VectorType>::multiply_with_AT(Vector<Number> &y, const VectorType &x) const
865 {
866  multiply_with_QT(y, x);
867 
868  const int N = this->current_size;
869  const int lda = N;
870  const int incx = 1;
871  trmv("U", "T", "N", &N, &this->R(0, 0), &lda, &y[0], &incx);
872 }
873 
874 #endif // no DOXYGEN
875 
876 DEAL_II_NAMESPACE_CLOSE
877 
878 #endif
unsigned int size() const
QR()
VectorType::value_type Number
Definition: qr.h:49
std::vector< std::unique_ptr< VectorType > > columns
Definition: qr.h:160
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)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
virtual ~BaseQR()=default
virtual void multiply_with_A(VectorType &y, const Vector< Number > &x) const =0
unsigned int current_size
Definition: qr.h:170
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)
virtual void multiply_with_Q(VectorType &y, const Vector< Number > &x) const
static ::ExceptionBase & ExcDivideByZero()
LAPACKFullMatrix< Number > R
Definition: qr.h:165
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:430
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:326
#define Assert(cond, exc)
Definition: exceptions.h:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
VectorType::value_type Number
Definition: qr.h:355
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:179
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:349
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:246
Definition: qr.h:240
void apply_givens_rotation(const unsigned int i, const unsigned int k)
Definition: qr.h:44
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