Reference documentation for deal.II version Git 409ee4b167 2020-08-14 09:46:12 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
qr.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2020 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 
22 
24 #include <deal.II/lac/utilities.h>
25 
26 #include <boost/signals2/signal.hpp>
27 
28 #include <memory>
29 
31 
43 template <typename VectorType>
44 class BaseQR
45 {
49  using Number = typename VectorType::value_type;
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  using Number = typename VectorType::value_type;
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 
327 };
328 
329 
330 
348 template <typename VectorType>
349 class ImplicitQR : public BaseQR<VectorType>
350 {
351 public:
355  using Number = typename VectorType::value_type;
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 namespace internal
437 {
438  namespace QRImplementation
439  {
440  // We want to avoid including our own LAPACK wrapper header in any external
441  // headers to avoid possible conflicts with other packages that may define
442  // their own such header. At the same time we want to be able to call some
443  // LAPACK functions from the template functions below. To resolve both
444  // problems define some extra wrappers here that can be in the header:
445  template <typename Number>
446  void
447  call_trmv(const char uplo,
448  const char trans,
449  const char diag,
450  const types::blas_int n,
451  const Number * a,
452  const types::blas_int lda,
453  Number * x,
454  const types::blas_int incx);
455 
456  template <typename Number>
457  void
458  call_trtrs(const char uplo,
459  const char trans,
460  const char diag,
461  const types::blas_int n,
462  const types::blas_int nrhs,
463  const Number * a,
464  const types::blas_int lda,
465  Number * b,
466  const types::blas_int ldb,
467  types::blas_int * info);
468  } // namespace QRImplementation
469 } // namespace internal
470 
471 
472 
473 template <typename VectorType>
475  : current_size(0)
476 {
478 }
479 
480 
481 
482 template <typename VectorType>
483 unsigned int
485 {
486  return current_size;
487 }
488 
489 
490 
491 template <typename VectorType>
494 {
495  return R;
496 }
497 
498 
499 
500 template <typename VectorType>
501 void
502 BaseQR<VectorType>::solve(Vector<Number> & x,
503  const Vector<Number> &y,
504  const bool transpose) const
505 {
506  Assert(x.size() == this->current_size,
507  ExcDimensionMismatch(x.size(), this->current_size));
508  Assert(y.size() == this->current_size,
509  ExcDimensionMismatch(y.size(), this->current_size));
510 
511  // copy if the two vectors are not the same
512  if (&x != &y)
513  x = y;
514 
515  const int lda = this->current_size;
516  const int ldb = this->current_size;
517  const int N = this->current_size;
518  const int n_rhs = 1;
519  int info = 0;
521  transpose ? 'T' : 'N',
522  'N',
523  N,
524  n_rhs,
525  &this->R(0, 0),
526  lda,
527  &x(0),
528  ldb,
529  &info);
530 }
531 
532 
533 
534 template <typename VectorType>
535 void
537  const Vector<Number> &x) const
538 {
539  Assert(x.size() == this->current_size,
540  ExcDimensionMismatch(x.size(), this->current_size));
541 
542  y = 0.;
543  for (unsigned int j = 0; j < this->current_size; ++j)
544  y.add(x[j], *this->columns[j]);
545 }
546 
547 
548 
549 template <typename VectorType>
550 void
551 BaseQR<VectorType>::multiply_with_colsT(Vector<Number> & y,
552  const VectorType &x) const
553 {
554  Assert(y.size() == this->current_size,
555  ExcDimensionMismatch(y.size(), this->current_size));
556 
557  for (unsigned int j = 0; j < this->current_size; ++j)
558  y[j] = (*this->columns[j]) * x;
559 }
560 
561 
562 
563 template <class VectorType>
564 boost::signals2::connection
566  const std::function<void(const unsigned int i,
567  const unsigned int j,
568  const std::array<Number, 3> &)> &slot)
569 {
570  return givens_signal.connect(slot);
571 }
572 
573 
574 
575 template <class VectorType>
576 boost::signals2::connection
578  const std::function<bool(const Vector<Number> &u,
579  const Number & rho,
580  const Number & col_norm_sqr)> &slot)
581 {
582  return column_signal.connect(slot);
583 }
584 
585 
586 
587 template <typename VectorType>
590 {}
591 
592 
593 
594 template <typename VectorType>
595 bool
597 {
598  if (this->current_size == 0)
599  {
600  this->R.grow_or_shrink(this->current_size + 1);
601  this->columns.push_back(std::make_unique<VectorType>(column));
602  this->R(0, 0) = column.l2_norm();
603  ++this->current_size;
604  }
605  else
606  {
607  // first get scalar products with A^T
608  Vector<Number> u(this->current_size);
609  this->multiply_with_AT(u, column);
610 
611  // now solve R^T x = (A^T * column)
612  const int lda = this->current_size;
613  const int ldb = this->current_size;
614  const int N = this->current_size;
615  const int n_rhs = 1;
616  int info = 0;
618  'U', 'T', 'N', N, n_rhs, &this->R(0, 0), lda, &u(0), ldb, &info);
619 
620  // finally get the diagonal element:
621  // rho2 = |column|^2 - |u|^2
622  const Number column_norm_sqr = column.norm_sqr();
623  const Number rho2 = column_norm_sqr - u.norm_sqr();
624  const bool linearly_independent =
625  column_signal.empty() ? rho2 > 0 :
626  column_signal(u, rho2, column_norm_sqr).get();
627 
628  // bail out if it turns out to be linearly dependent
629  if (!linearly_independent)
630  return false;
631 
632  // at this point we update is successful and we can enlarge R
633  // and store the column:
634  this->columns.push_back(std::make_unique<VectorType>(column));
635  this->R.grow_or_shrink(this->current_size + 1);
636  this->R(this->current_size, this->current_size) = std::sqrt(rho2);
637  for (unsigned int i = 0; i < this->current_size; ++i)
638  this->R(i, this->current_size) = u(i);
639 
640  this->current_size++;
641  }
642 
643  return true;
644 }
645 
646 
647 
648 template <typename VectorType>
649 void
651  const unsigned int k)
652 {
653  AssertIndexRange(i, k);
654  AssertIndexRange(k, this->current_size);
655  const std::array<Number, 3> csr =
656  ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
657  this->R(k, k));
658 
659  // first, set k'th column:
660  this->R(i, k) = csr[2];
661  this->R(k, k) = 0.;
662  // now do the rest:
663  for (unsigned int j = 0; j < this->R.n(); ++j)
664  if (j != k)
665  {
666  const Number t = this->R(i, j);
667  this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
668  this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
669  }
670 
671  if (!this->givens_signal.empty())
672  this->givens_signal(i, k, csr);
673 }
674 
675 
676 
677 template <typename VectorType>
678 void
679 ImplicitQR<VectorType>::remove_column(const unsigned int k)
680 {
681  // before actually removing a column from Q and resizing R,
682  // apply givens rotations to bring H into upper triangular form:
683  for (unsigned int j = k + 1; j < this->R.n(); ++j)
684  {
685  const unsigned int i = j - 1;
686  apply_givens_rotation(i, j);
687  }
688 
689  // remove last row and k-th column
690  --this->current_size;
691  this->R.remove_row_and_column(this->current_size, k);
692 
693  // Finally remove the column from A
694  this->columns.erase(this->columns.begin() + k);
695 }
696 
697 
698 
699 template <typename VectorType>
700 void
702  const Vector<Number> &x) const
703 {
704  // A = QR
705  // A R^{-1} = Q
706  Vector<Number> x1 = x;
707  BaseQR<VectorType>::solve(x1, x1, false);
708  multiply_with_A(y, x1);
709 }
710 
711 
712 
713 template <typename VectorType>
714 void
716  const VectorType &x) const
717 {
718  // A = QR
719  // A^T = R^T Q^T
720  // {R^T}^{-1} A^T = Q^T
721  multiply_with_AT(y, x);
722  BaseQR<VectorType>::solve(y, y, true);
723 }
724 
725 
726 
727 template <typename VectorType>
728 void
730  const Vector<Number> &x) const
731 {
733 }
734 
735 
736 
737 template <typename VectorType>
738 void
740  const VectorType &x) const
741 {
743 }
744 
745 
746 
747 template <typename VectorType>
750 {}
751 
752 
753 
754 template <typename VectorType>
755 bool
757 {
758  // resize R:
759  this->R.grow_or_shrink(this->current_size + 1);
760  this->columns.push_back(std::make_unique<VectorType>(column));
761 
762  // now a Gram-Schmidt part: orthonormalize the new column
763  // against everything we have so far:
764  auto &last_col = *this->columns.back();
765  for (unsigned int i = 0; i < this->current_size; ++i)
766  {
767  const auto &i_col = *this->columns[i];
768  this->R(i, this->current_size) = i_col * last_col;
769  last_col.add(-this->R(i, this->current_size), i_col);
770  }
771 
772  this->R(this->current_size, this->current_size) = last_col.l2_norm();
773 
774  Assert(this->R(this->current_size, this->current_size) > 0.,
775  ExcDivideByZero());
776  last_col *= 1. / this->R(this->current_size, this->current_size);
777 
778  ++this->current_size;
779  return true;
780 }
781 
782 
783 
784 template <typename VectorType>
785 void
786 QR<VectorType>::apply_givens_rotation(const unsigned int i,
787  const unsigned int k)
788 {
789  AssertIndexRange(i, k);
790  AssertIndexRange(k, this->current_size);
791  const std::array<Number, 3> csr =
792  ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
793  this->R(k, k));
794 
795  // first, set k'th column:
796  this->R(i, k) = csr[2];
797  this->R(k, k) = 0.;
798  // now do the rest:
799  for (unsigned int j = 0; j < this->R.n(); ++j)
800  if (j != k)
801  {
802  const Number t = this->R(i, j);
803  this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
804  this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
805  }
806 
807  // now adjust i,k columns due to multiplication with the
808  // transpose Givens matrix from right:
809  auto &col_i = *this->columns[i];
810  auto &col_k = *this->columns[k];
811  // save column i:
812  tmp = col_i;
813  col_i.sadd(csr[0], csr[1], col_k);
814  col_k.sadd(csr[0], -csr[1], tmp);
815 
816  if (!this->givens_signal.empty())
817  this->givens_signal(i, k, csr);
818 }
819 
820 
821 
822 template <typename VectorType>
823 void
824 QR<VectorType>::remove_column(const unsigned int k)
825 {
826  AssertIndexRange(k, this->current_size);
827  Assert(this->current_size > 0,
828  ExcMessage("Can not remove a column if QR is empty"));
829  // apply a sequence of Givens rotations
830  // see section 6.5 "Updating matrix factorizations" in Golub 2013, Matrix
831  // computations
832 
833  // So we want to have QR for \tilde A \in R^{m*(n-1)}
834  // if we remove the column k, we end up with upper Hessenberg matrix
835  // x x x x x
836  // x x x x
837  // H = x x x
838  // x x x
839  // x x
840  // x
841  // where k = 2 (3rd column), m = 7, n = 6
842  //
843  // before actually removing a column from Q and resizing R,
844  // apply givens rotations to bring H into upper triangular form:
845  for (unsigned int j = k + 1; j < this->R.n(); ++j)
846  {
847  const unsigned int i = j - 1;
848  apply_givens_rotation(i, j);
849  }
850 
851  // now we can throw away the column from Q and adjust R
852  // since we do thin-QR, after Givens rotations we need to throw
853  // away the last column:
854  const unsigned int size_minus_1 = this->columns.size() - 1;
855  this->columns.erase(this->columns.begin() + size_minus_1);
856 
857  // remove last row and k-th column
858  --this->current_size;
859  this->R.remove_row_and_column(this->current_size, k);
860 }
861 
862 
863 
864 template <typename VectorType>
865 void
866 QR<VectorType>::multiply_with_Q(VectorType &y, const Vector<Number> &x) const
867 {
869 }
870 
871 
872 
873 template <typename VectorType>
874 void
875 QR<VectorType>::multiply_with_QT(Vector<Number> &y, const VectorType &x) const
876 {
878 }
879 
880 
881 
882 template <typename VectorType>
883 void
884 QR<VectorType>::multiply_with_A(VectorType &y, const Vector<Number> &x) const
885 {
886  Vector<Number> x1 = x;
887  const int N = this->current_size;
888  const int lda = N;
889  const int incx = 1;
891  'U', 'N', 'N', N, &this->R(0, 0), lda, &x1[0], incx);
892 
893  multiply_with_Q(y, x1);
894 }
895 
896 
897 
898 template <typename VectorType>
899 void
900 QR<VectorType>::multiply_with_AT(Vector<Number> &y, const VectorType &x) const
901 {
902  multiply_with_QT(y, x);
903 
904  const int N = this->current_size;
905  const int lda = N;
906  const int incx = 1;
908  'U', 'T', 'N', N, &this->R(0, 0), lda, &y[0], incx);
909 }
910 
911 #endif // no DOXYGEN
912 
914 
915 #endif
unsigned int size() const
QR()
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:1636
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
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::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)
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
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
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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
Definition: qr.h:240
void apply_givens_rotation(const unsigned int i, const unsigned int k)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
Definition: qr.h:44
void set_property(const LAPACKSupport::Property property)
static const char N
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const =0
void multiply_with_cols(VectorType &y, const Vector< Number > &x) const
typename VectorType::value_type Number
Definition: qr.h:49
void call_trtrs(const char uplo, const char trans, const char diag, const types::blas_int n, const types::blas_int nrhs, const Number *a, const types::blas_int lda, Number *b, const types::blas_int ldb, types::blas_int *info)
Definition: qr.cc:46
void call_trmv(const char uplo, const char trans, const char diag, const types::blas_int n, const Number *a, const types::blas_int lda, Number *x, const types::blas_int incx)
Definition: qr.cc:32