Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2019 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_qr_h
16#define dealii_qr_h
17
18#include <deal.II/base/config.h>
19
21
24
25#include <boost/signals2/signal.hpp>
26
27#include <memory>
28
30
42template <typename VectorType>
43class BaseQR
44{
48 using Number = typename VectorType::value_type;
49
50protected:
55
56public:
60 virtual ~BaseQR() = default;
61
68 virtual bool
69 append_column(const VectorType &column) = 0;
70
74 virtual void
75 remove_column(const unsigned int k = 0) = 0;
76
80 unsigned int
81 size() const;
82
87 get_R() const;
88
94 void
96 const Vector<Number> &y,
97 const bool transpose = false) const;
98
103 virtual void
104 multiply_with_Q(VectorType &y, const Vector<Number> &x) const = 0;
105
110 virtual void
111 multiply_with_QT(Vector<Number> &y, const VectorType &x) const = 0;
112
117 virtual void
118 multiply_with_A(VectorType &y, const Vector<Number> &x) const = 0;
119
124 virtual void
125 multiply_with_AT(Vector<Number> &y, const VectorType &x) const = 0;
126
136 boost::signals2::connection
138 const std::function<void(const unsigned int i,
139 const unsigned int j,
140 const std::array<Number, 3> &csr)> &slot);
141
142protected:
147 void
148 multiply_with_cols(VectorType &y, const Vector<Number> &x) const;
149
153 void
154 multiply_with_colsT(Vector<Number> &y, const VectorType &x) const;
155
159 std::vector<std::unique_ptr<VectorType>> columns;
160
165
169 unsigned int current_size;
170
175 boost::signals2::signal<void(const unsigned int i,
176 const unsigned int j,
177 const std::array<Number, 3> &)>
179};
180
181// clang-format off
237// clang-format on
238template <typename VectorType>
239class QR : public BaseQR<VectorType>
240{
241public:
245 using Number = typename VectorType::value_type;
246
250 QR();
251
255 virtual ~QR() = default;
256
262 virtual bool
263 append_column(const VectorType &column);
264
297 virtual void
298 remove_column(const unsigned int k = 0);
299
300 virtual void
301 multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
302
303 virtual void
304 multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
305
306 virtual void
307 multiply_with_A(VectorType &y, const Vector<Number> &x) const;
308
309 virtual void
310 multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
311
312private:
319 void
320 apply_givens_rotation(const unsigned int i, const unsigned int k);
321
325 VectorType tmp;
326};
327
328
329
347template <typename VectorType>
348class ImplicitQR : public BaseQR<VectorType>
349{
350public:
354 using Number = typename VectorType::value_type;
355
360
364 virtual ~ImplicitQR() = default;
365
366 virtual bool
367 append_column(const VectorType &column);
368
381 virtual void
382 remove_column(const unsigned int k = 0);
383
384 virtual void
385 multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
386
387 virtual void
388 multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
389
390 virtual void
391 multiply_with_A(VectorType &y, const Vector<Number> &x) const;
392
393 virtual void
394 multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
395
405 boost::signals2::connection
407 const std::function<bool(const Vector<Number> &u,
408 const Number &rho2,
409 const Number &col_norm_sqr)> &slot);
410
411private:
415 void
416 apply_givens_rotation(const unsigned int i, const unsigned int k);
417
426 boost::signals2::signal<bool(const Vector<Number> &u,
427 const Number &rho,
428 const Number &col_norm_sqr)>
430};
431
432// ------------------- inline and template functions ----------------
433#ifndef DOXYGEN
434
435namespace internal
436{
437 namespace QRImplementation
438 {
439 // We want to avoid including our own LAPACK wrapper header in any external
440 // headers to avoid possible conflicts with other packages that may define
441 // their own such header. At the same time we want to be able to call some
442 // LAPACK functions from the template functions below. To resolve both
443 // problems define some extra wrappers here that can be in the header:
444 template <typename Number>
445 void
446 call_trmv(const char uplo,
447 const char trans,
448 const char diag,
449 const types::blas_int n,
450 const Number *a,
451 const types::blas_int lda,
452 Number *x,
453 const types::blas_int incx);
454
455 template <typename Number>
456 void
457 call_trtrs(const char uplo,
458 const char trans,
459 const char diag,
460 const types::blas_int n,
461 const types::blas_int nrhs,
462 const Number *a,
463 const types::blas_int lda,
464 Number *b,
465 const types::blas_int ldb,
466 types::blas_int *info);
467 } // namespace QRImplementation
468} // namespace internal
469
470
471
472template <typename VectorType>
474 : current_size(0)
475{
477}
478
479
480
481template <typename VectorType>
482unsigned int
484{
485 return current_size;
486}
487
488
489
490template <typename VectorType>
493{
494 return R;
495}
496
497
498
499template <typename VectorType>
500void
502 const Vector<Number> &y,
503 const bool transpose) const
504{
505 Assert(x.size() == this->current_size,
506 ExcDimensionMismatch(x.size(), this->current_size));
507 Assert(y.size() == this->current_size,
508 ExcDimensionMismatch(y.size(), this->current_size));
509
510 // copy if the two vectors are not the same
511 if (&x != &y)
512 x = y;
513
514 const int lda = this->current_size;
515 const int ldb = this->current_size;
516 const int N = this->current_size;
517 const int n_rhs = 1;
518 int info = 0;
520 transpose ? 'T' : 'N',
521 'N',
522 N,
523 n_rhs,
524 &this->R(0, 0),
525 lda,
526 &x(0),
527 ldb,
528 &info);
529}
530
531
532
533template <typename VectorType>
534void
536 const Vector<Number> &x) const
537{
538 Assert(x.size() == this->current_size,
539 ExcDimensionMismatch(x.size(), this->current_size));
540
541 y = 0.;
542 for (unsigned int j = 0; j < this->current_size; ++j)
543 y.add(x[j], *this->columns[j]);
544}
545
546
547
548template <typename VectorType>
549void
551 const VectorType &x) const
552{
553 Assert(y.size() == this->current_size,
554 ExcDimensionMismatch(y.size(), this->current_size));
555
556 for (unsigned int j = 0; j < this->current_size; ++j)
557 y[j] = (*this->columns[j]) * x;
558}
559
560
561
562template <typename VectorType>
563boost::signals2::connection
565 const std::function<void(const unsigned int i,
566 const unsigned int j,
567 const std::array<Number, 3> &)> &slot)
568{
569 return givens_signal.connect(slot);
570}
571
572
573
574template <typename VectorType>
575boost::signals2::connection
577 const std::function<bool(const Vector<Number> &u,
578 const Number &rho,
579 const Number &col_norm_sqr)> &slot)
580{
581 return column_signal.connect(slot);
582}
583
584
585
586template <typename VectorType>
588 : BaseQR<VectorType>()
589{}
590
591
592
593template <typename VectorType>
594bool
595ImplicitQR<VectorType>::append_column(const VectorType &column)
596{
597 if (this->current_size == 0)
598 {
599 this->R.grow_or_shrink(this->current_size + 1);
600 this->columns.push_back(std::make_unique<VectorType>(column));
601 this->R(0, 0) = column.l2_norm();
602 ++this->current_size;
603 }
604 else
605 {
606 // first get scalar products with A^T
607 Vector<Number> u(this->current_size);
608 this->multiply_with_AT(u, column);
609
610 // now solve R^T x = (A^T * column)
611 const int lda = this->current_size;
612 const int ldb = this->current_size;
613 const int N = this->current_size;
614 const int n_rhs = 1;
615 int info = 0;
617 'U', 'T', 'N', N, n_rhs, &this->R(0, 0), lda, &u(0), ldb, &info);
618
619 // finally get the diagonal element:
620 // rho2 = |column|^2 - |u|^2
621 const Number column_norm_sqr = column.norm_sqr();
622 const Number rho2 = column_norm_sqr - u.norm_sqr();
623 const bool linearly_independent =
624 column_signal.empty() ? rho2 > 0 :
625 column_signal(u, rho2, column_norm_sqr).get();
626
627 // bail out if it turns out to be linearly dependent
628 if (!linearly_independent)
629 return false;
630
631 // at this point we update is successful and we can enlarge R
632 // and store the column:
633 this->columns.push_back(std::make_unique<VectorType>(column));
634 this->R.grow_or_shrink(this->current_size + 1);
635 this->R(this->current_size, this->current_size) = std::sqrt(rho2);
636 for (unsigned int i = 0; i < this->current_size; ++i)
637 this->R(i, this->current_size) = u(i);
638
639 this->current_size++;
640 }
641
642 return true;
643}
644
645
646
647template <typename VectorType>
648void
650 const unsigned int k)
651{
652 AssertIndexRange(i, k);
653 AssertIndexRange(k, this->current_size);
654 const std::array<Number, 3> csr =
655 ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
656 this->R(k, k));
657
658 // first, set k'th column:
659 this->R(i, k) = csr[2];
660 this->R(k, k) = 0.;
661 // now do the rest:
662 for (unsigned int j = 0; j < this->R.n(); ++j)
663 if (j != k)
664 {
665 const Number t = this->R(i, j);
666 this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
667 this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
668 }
669
670 if (!this->givens_signal.empty())
671 this->givens_signal(i, k, csr);
672}
673
674
675
676template <typename VectorType>
677void
678ImplicitQR<VectorType>::remove_column(const unsigned int k)
679{
680 // before actually removing a column from Q and resizing R,
681 // apply givens rotations to bring H into upper triangular form:
682 for (unsigned int j = k + 1; j < this->R.n(); ++j)
683 {
684 const unsigned int i = j - 1;
685 apply_givens_rotation(i, j);
686 }
687
688 // remove last row and k-th column
689 --this->current_size;
690 this->R.remove_row_and_column(this->current_size, k);
691
692 // Finally remove the column from A
693 this->columns.erase(this->columns.begin() + k);
694}
695
696
697
698template <typename VectorType>
699void
701 const Vector<Number> &x) const
702{
703 // A = QR
704 // A R^{-1} = Q
705 Vector<Number> x1 = x;
706 BaseQR<VectorType>::solve(x1, x1, false);
707 multiply_with_A(y, x1);
708}
709
710
711
712template <typename VectorType>
713void
715 const VectorType &x) const
716{
717 // A = QR
718 // A^T = R^T Q^T
719 // {R^T}^{-1} A^T = Q^T
720 multiply_with_AT(y, x);
721 BaseQR<VectorType>::solve(y, y, true);
722}
723
724
725
726template <typename VectorType>
727void
729 const Vector<Number> &x) const
730{
732}
733
734
735
736template <typename VectorType>
737void
739 const VectorType &x) const
740{
742}
743
744
745
746template <typename VectorType>
748 : BaseQR<VectorType>()
749{}
750
751
752
753template <typename VectorType>
754bool
755QR<VectorType>::append_column(const VectorType &column)
756{
757 // resize R:
758 this->R.grow_or_shrink(this->current_size + 1);
759 this->columns.push_back(std::make_unique<VectorType>(column));
760
761 // now a Gram-Schmidt part: orthonormalize the new column
762 // against everything we have so far:
763 auto &last_col = *this->columns.back();
764 for (unsigned int i = 0; i < this->current_size; ++i)
765 {
766 const auto &i_col = *this->columns[i];
767 this->R(i, this->current_size) = i_col * last_col;
768 last_col.add(-this->R(i, this->current_size), i_col);
769 }
770
771 this->R(this->current_size, this->current_size) = last_col.l2_norm();
772
773 Assert(this->R(this->current_size, this->current_size) > 0.,
775 last_col *= 1. / this->R(this->current_size, this->current_size);
776
777 ++this->current_size;
778 return true;
779}
780
781
782
783template <typename VectorType>
784void
785QR<VectorType>::apply_givens_rotation(const unsigned int i,
786 const unsigned int k)
787{
788 AssertIndexRange(i, k);
789 AssertIndexRange(k, this->current_size);
790 const std::array<Number, 3> csr =
791 ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
792 this->R(k, k));
793
794 // first, set k'th column:
795 this->R(i, k) = csr[2];
796 this->R(k, k) = 0.;
797 // now do the rest:
798 for (unsigned int j = 0; j < this->R.n(); ++j)
799 if (j != k)
800 {
801 const Number t = this->R(i, j);
802 this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
803 this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
804 }
805
806 // now adjust i,k columns due to multiplication with the
807 // transpose Givens matrix from right:
808 auto &col_i = *this->columns[i];
809 auto &col_k = *this->columns[k];
810 // save column i:
811 tmp = col_i;
812 col_i.sadd(csr[0], csr[1], col_k);
813 col_k.sadd(csr[0], -csr[1], tmp);
814
815 if (!this->givens_signal.empty())
816 this->givens_signal(i, k, csr);
817}
818
819
820
821template <typename VectorType>
822void
823QR<VectorType>::remove_column(const unsigned int k)
824{
825 AssertIndexRange(k, this->current_size);
826 Assert(this->current_size > 0,
827 ExcMessage("Can not remove a column if QR is empty"));
828 // apply a sequence of Givens rotations
829 // see section 6.5 "Updating matrix factorizations" in Golub 2013, Matrix
830 // computations
831
832 // So we want to have QR for \tilde A \in R^{m*(n-1)}
833 // if we remove the column k, we end up with upper Hessenberg matrix
834 // x x x x x
835 // x x x x
836 // H = x x x
837 // x x x
838 // x x
839 // x
840 // where k = 2 (3rd column), m = 7, n = 6
841 //
842 // before actually removing a column from Q and resizing R,
843 // apply givens rotations to bring H into upper triangular form:
844 for (unsigned int j = k + 1; j < this->R.n(); ++j)
845 {
846 const unsigned int i = j - 1;
847 apply_givens_rotation(i, j);
848 }
849
850 // now we can throw away the column from Q and adjust R
851 // since we do thin-QR, after Givens rotations we need to throw
852 // away the last column:
853 const unsigned int size_minus_1 = this->columns.size() - 1;
854 this->columns.erase(this->columns.begin() + size_minus_1);
855
856 // remove last row and k-th column
857 --this->current_size;
858 this->R.remove_row_and_column(this->current_size, k);
859}
860
861
862
863template <typename VectorType>
864void
865QR<VectorType>::multiply_with_Q(VectorType &y, const Vector<Number> &x) const
866{
868}
869
870
871
872template <typename VectorType>
873void
874QR<VectorType>::multiply_with_QT(Vector<Number> &y, const VectorType &x) const
875{
877}
878
879
880
881template <typename VectorType>
882void
883QR<VectorType>::multiply_with_A(VectorType &y, const Vector<Number> &x) const
884{
885 Vector<Number> x1 = x;
886 const int N = this->current_size;
887 const int lda = N;
888 const int incx = 1;
890 'U', 'N', 'N', N, &this->R(0, 0), lda, &x1[0], incx);
891
892 multiply_with_Q(y, x1);
893}
894
895
896
897template <typename VectorType>
898void
899QR<VectorType>::multiply_with_AT(Vector<Number> &y, const VectorType &x) const
900{
901 multiply_with_QT(y, x);
902
903 const int N = this->current_size;
904 const int lda = N;
905 const int incx = 1;
907 'U', 'T', 'N', N, &this->R(0, 0), lda, &y[0], incx);
908}
909
910#endif // no DOXYGEN
911
913
914#endif
Definition qr.h:44
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:159
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:164
typename VectorType::value_type Number
Definition qr.h:48
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:178
void multiply_with_cols(VectorType &y, const Vector< Number > &x) const
unsigned int current_size
Definition qr.h:169
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:354
boost::signals2::signal< bool(const Vector< Number > &u, const Number &rho, const Number &col_norm_sqr)> column_signal
Definition qr.h:429
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)
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:240
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const
typename VectorType::value_type Number
Definition qr.h:245
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:325
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
virtual size_type size() const override
real_type norm_sqr() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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:31
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:45
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)