Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
qr.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2018 - 2021 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
25
26#include <boost/signals2/signal.hpp>
27
28#include <memory>
29
31
43template <typename VectorType>
44class BaseQR
45{
49 using Number = typename VectorType::value_type;
50
51protected:
56
57public:
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
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
143protected:
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
239template <typename VectorType>
240class QR : public BaseQR<VectorType>
241{
242public:
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
313private:
320 void
321 apply_givens_rotation(const unsigned int i, const unsigned int k);
322
326 VectorType tmp;
327};
328
329
330
348template <typename VectorType>
349class ImplicitQR : public BaseQR<VectorType>
350{
351public:
355 using Number = typename VectorType::value_type;
356
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
408 const std::function<bool(const Vector<Number> &u,
409 const Number & rho2,
410 const Number & col_norm_sqr)> &slot);
411
412private:
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)>
431};
432
433// ------------------- inline and template functions ----------------
434#ifndef DOXYGEN
435
436namespace 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
473template <typename VectorType>
475 : current_size(0)
476{
478}
479
480
481
482template <typename VectorType>
483unsigned int
485{
486 return current_size;
487}
488
489
490
491template <typename VectorType>
494{
495 return R;
496}
497
498
499
500template <typename VectorType>
501void
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
534template <typename VectorType>
535void
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
549template <typename VectorType>
550void
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
563template <class VectorType>
564boost::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
575template <class VectorType>
576boost::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
587template <typename VectorType>
589 : BaseQR<VectorType>()
590{}
591
592
593
594template <typename VectorType>
595bool
596ImplicitQR<VectorType>::append_column(const VectorType &column)
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
648template <typename VectorType>
649void
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
677template <typename VectorType>
678void
679ImplicitQR<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
699template <typename VectorType>
700void
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
713template <typename VectorType>
714void
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
727template <typename VectorType>
728void
730 const Vector<Number> &x) const
731{
733}
734
735
736
737template <typename VectorType>
738void
740 const VectorType &x) const
741{
743}
744
745
746
747template <typename VectorType>
749 : BaseQR<VectorType>()
750{}
751
752
753
754template <typename VectorType>
755bool
756QR<VectorType>::append_column(const VectorType &column)
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.,
776 last_col *= 1. / this->R(this->current_size, this->current_size);
777
778 ++this->current_size;
779 return true;
780}
781
782
783
784template <typename VectorType>
785void
786QR<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
822template <typename VectorType>
823void
824QR<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
864template <typename VectorType>
865void
866QR<VectorType>::multiply_with_Q(VectorType &y, const Vector<Number> &x) const
867{
869}
870
871
872
873template <typename VectorType>
874void
875QR<VectorType>::multiply_with_QT(Vector<Number> &y, const VectorType &x) const
876{
878}
879
880
881
882template <typename VectorType>
883void
884QR<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
898template <typename VectorType>
899void
900QR<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
Definition qr.h:45
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)
void multiply_with_colsT(Vector< Number > &y, const VectorType &x) const
unsigned int size() const
virtual void remove_column(const unsigned int k=0)=0
virtual ~BaseQR()=default
virtual void multiply_with_Q(VectorType &y, const Vector< Number > &x) const =0
virtual void multiply_with_A(VectorType &y, const Vector< Number > &x) const =0
std::vector< std::unique_ptr< VectorType > > columns
Definition qr.h:160
virtual void multiply_with_AT(Vector< Number > &y, const VectorType &x) const =0
void solve(Vector< Number > &x, const Vector< Number > &y, const bool transpose=false) const
LAPACKFullMatrix< Number > R
Definition qr.h:165
typename VectorType::value_type Number
Definition qr.h:49
virtual bool append_column(const VectorType &column)=0
boost::signals2::signal< void(const unsigned int i, const unsigned int j, const std::array< Number, 3 > &)> givens_signal
Definition qr.h:179
void multiply_with_cols(VectorType &y, const Vector< Number > &x) const
unsigned int current_size
Definition qr.h:170
const LAPACKFullMatrix< Number > & get_R() const
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const =0
virtual ~ImplicitQR()=default
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const
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)
typename VectorType::value_type Number
Definition qr.h:355
virtual bool append_column(const VectorType &column)
virtual void multiply_with_Q(VectorType &y, const Vector< Number > &x) const
virtual void remove_column(const unsigned int k=0)
boost::signals2::signal< bool(const Vector< Number > &u, const Number &rho, const Number &col_norm_sqr)> column_signal
Definition qr.h:430
void apply_givens_rotation(const unsigned int i, const unsigned int k)
virtual void multiply_with_AT(Vector< Number > &y, const VectorType &x) const
void set_property(const LAPACKSupport::Property property)
Definition qr.h:241
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const
typename VectorType::value_type Number
Definition qr.h:246
virtual bool append_column(const VectorType &column)
void apply_givens_rotation(const unsigned int i, const unsigned int k)
virtual ~QR()=default
virtual void multiply_with_A(VectorType &y, const Vector< Number > &x) const
VectorType tmp
Definition qr.h:326
virtual void multiply_with_AT(Vector< Number > &y, const VectorType &x) const
virtual void remove_column(const unsigned int k=0)
QR()
virtual void multiply_with_Q(VectorType &y, const Vector< Number > &x) const
size_type size() const
real_type norm_sqr() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
@ upper_triangular
Matrix is upper triangular.
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
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
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)