deal.II version GIT relicensing-1929-g70a5450eef 2024-10-03 18:00:00+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 - 2024 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
199// clang-format on
200template <typename VectorType>
201class QR : public BaseQR<VectorType>
202{
203public:
207 using Number = typename VectorType::value_type;
208
212 QR();
213
217 virtual ~QR() = default;
218
224 virtual bool
225 append_column(const VectorType &column);
226
259 virtual void
260 remove_column(const unsigned int k = 0);
261
262 virtual void
263 multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
264
265 virtual void
266 multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
267
268 virtual void
269 multiply_with_A(VectorType &y, const Vector<Number> &x) const;
270
271 virtual void
272 multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
273
274private:
281 void
282 apply_givens_rotation(const unsigned int i, const unsigned int k);
283
287 VectorType tmp;
288};
289
290
291
309template <typename VectorType>
310class ImplicitQR : public BaseQR<VectorType>
311{
312public:
316 using Number = typename VectorType::value_type;
317
322
326 virtual ~ImplicitQR() = default;
327
328 virtual bool
329 append_column(const VectorType &column);
330
343 virtual void
344 remove_column(const unsigned int k = 0);
345
346 virtual void
347 multiply_with_Q(VectorType &y, const Vector<Number> &x) const;
348
349 virtual void
350 multiply_with_QT(Vector<Number> &y, const VectorType &x) const;
351
352 virtual void
353 multiply_with_A(VectorType &y, const Vector<Number> &x) const;
354
355 virtual void
356 multiply_with_AT(Vector<Number> &y, const VectorType &x) const;
357
367 boost::signals2::connection
369 const std::function<bool(const Vector<Number> &u,
370 const Number &rho2,
371 const Number &col_norm_sqr)> &slot);
372
373private:
377 void
378 apply_givens_rotation(const unsigned int i, const unsigned int k);
379
388 boost::signals2::signal<bool(const Vector<Number> &u,
389 const Number &rho,
390 const Number &col_norm_sqr)>
392};
393
394// ------------------- inline and template functions ----------------
395#ifndef DOXYGEN
396
397namespace internal
398{
399 namespace QRImplementation
400 {
401 // We want to avoid including our own LAPACK wrapper header in any external
402 // headers to avoid possible conflicts with other packages that may define
403 // their own such header. At the same time we want to be able to call some
404 // LAPACK functions from the template functions below. To resolve both
405 // problems define some extra wrappers here that can be in the header:
406 template <typename Number>
407 void
408 call_trmv(const char uplo,
409 const char trans,
410 const char diag,
411 const types::blas_int n,
412 const Number *a,
413 const types::blas_int lda,
414 Number *x,
415 const types::blas_int incx);
416
417 template <typename Number>
418 void
419 call_trtrs(const char uplo,
420 const char trans,
421 const char diag,
422 const types::blas_int n,
423 const types::blas_int nrhs,
424 const Number *a,
425 const types::blas_int lda,
426 Number *b,
427 const types::blas_int ldb,
428 types::blas_int *info);
429 } // namespace QRImplementation
430} // namespace internal
431
432
433
434template <typename VectorType>
436 : current_size(0)
437{
439}
440
441
442
443template <typename VectorType>
444unsigned int
446{
447 return current_size;
448}
449
450
451
452template <typename VectorType>
455{
456 return R;
457}
458
459
460
461template <typename VectorType>
462void
464 const Vector<Number> &y,
465 const bool transpose) const
466{
467 Assert(x.size() == this->current_size,
468 ExcDimensionMismatch(x.size(), this->current_size));
469 Assert(y.size() == this->current_size,
470 ExcDimensionMismatch(y.size(), this->current_size));
471
472 // copy if the two vectors are not the same
473 if (&x != &y)
474 x = y;
475
476 const int lda = this->current_size;
477 const int ldb = this->current_size;
478 const int N = this->current_size;
479 const int n_rhs = 1;
480 int info = 0;
482 transpose ? 'T' : 'N',
483 'N',
484 N,
485 n_rhs,
486 &this->R(0, 0),
487 lda,
488 &x(0),
489 ldb,
490 &info);
491}
492
493
494
495template <typename VectorType>
496void
498 const Vector<Number> &x) const
499{
500 Assert(x.size() == this->current_size,
501 ExcDimensionMismatch(x.size(), this->current_size));
502
503 y = 0.;
504 for (unsigned int j = 0; j < this->current_size; ++j)
505 y.add(x[j], *this->columns[j]);
506}
507
508
509
510template <typename VectorType>
511void
513 const VectorType &x) const
514{
515 Assert(y.size() == this->current_size,
516 ExcDimensionMismatch(y.size(), this->current_size));
517
518 for (unsigned int j = 0; j < this->current_size; ++j)
519 y[j] = (*this->columns[j]) * x;
520}
521
522
523
524template <typename VectorType>
525boost::signals2::connection
527 const std::function<void(const unsigned int i,
528 const unsigned int j,
529 const std::array<Number, 3> &)> &slot)
530{
531 return givens_signal.connect(slot);
532}
533
534
535
536template <typename VectorType>
537boost::signals2::connection
539 const std::function<bool(const Vector<Number> &u,
540 const Number &rho,
541 const Number &col_norm_sqr)> &slot)
542{
543 return column_signal.connect(slot);
544}
545
546
547
548template <typename VectorType>
550 : BaseQR<VectorType>()
551{}
552
553
554
555template <typename VectorType>
556bool
557ImplicitQR<VectorType>::append_column(const VectorType &column)
558{
559 if (this->current_size == 0)
560 {
561 this->R.grow_or_shrink(this->current_size + 1);
562 this->columns.push_back(std::make_unique<VectorType>(column));
563 this->R(0, 0) = column.l2_norm();
564 ++this->current_size;
565 }
566 else
567 {
568 // first get scalar products with A^T
569 Vector<Number> u(this->current_size);
570 this->multiply_with_AT(u, column);
571
572 // now solve R^T x = (A^T * column)
573 const int lda = this->current_size;
574 const int ldb = this->current_size;
575 const int N = this->current_size;
576 const int n_rhs = 1;
577 int info = 0;
579 'U', 'T', 'N', N, n_rhs, &this->R(0, 0), lda, &u(0), ldb, &info);
580
581 // finally get the diagonal element:
582 // rho2 = |column|^2 - |u|^2
583 const Number column_norm_sqr = column.norm_sqr();
584 const Number rho2 = column_norm_sqr - u.norm_sqr();
585 const bool linearly_independent =
586 column_signal.empty() ? rho2 > 0 :
587 column_signal(u, rho2, column_norm_sqr).get();
588
589 // bail out if it turns out to be linearly dependent
590 if (!linearly_independent)
591 return false;
592
593 // at this point we update is successful and we can enlarge R
594 // and store the column:
595 this->columns.push_back(std::make_unique<VectorType>(column));
596 this->R.grow_or_shrink(this->current_size + 1);
597 this->R(this->current_size, this->current_size) = std::sqrt(rho2);
598 for (unsigned int i = 0; i < this->current_size; ++i)
599 this->R(i, this->current_size) = u(i);
600
601 this->current_size++;
602 }
603
604 return true;
605}
606
607
608
609template <typename VectorType>
610void
612 const unsigned int k)
613{
614 AssertIndexRange(i, k);
615 AssertIndexRange(k, this->current_size);
616 const std::array<Number, 3> csr =
617 ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
618 this->R(k, k));
619
620 // first, set k'th column:
621 this->R(i, k) = csr[2];
622 this->R(k, k) = 0.;
623 // now do the rest:
624 for (unsigned int j = 0; j < this->R.n(); ++j)
625 if (j != k)
626 {
627 const Number t = this->R(i, j);
628 this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
629 this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
630 }
631
632 if (!this->givens_signal.empty())
633 this->givens_signal(i, k, csr);
634}
635
636
637
638template <typename VectorType>
639void
640ImplicitQR<VectorType>::remove_column(const unsigned int k)
641{
642 // before actually removing a column from Q and resizing R,
643 // apply givens rotations to bring H into upper triangular form:
644 for (unsigned int j = k + 1; j < this->R.n(); ++j)
645 {
646 const unsigned int i = j - 1;
647 apply_givens_rotation(i, j);
648 }
649
650 // remove last row and k-th column
651 --this->current_size;
652 this->R.remove_row_and_column(this->current_size, k);
653
654 // Finally remove the column from A
655 this->columns.erase(this->columns.begin() + k);
656}
657
658
659
660template <typename VectorType>
661void
663 const Vector<Number> &x) const
664{
665 // A = QR
666 // A R^{-1} = Q
667 Vector<Number> x1 = x;
668 BaseQR<VectorType>::solve(x1, x1, false);
669 multiply_with_A(y, x1);
670}
671
672
673
674template <typename VectorType>
675void
677 const VectorType &x) const
678{
679 // A = QR
680 // A^T = R^T Q^T
681 // {R^T}^{-1} A^T = Q^T
682 multiply_with_AT(y, x);
683 BaseQR<VectorType>::solve(y, y, true);
684}
685
686
687
688template <typename VectorType>
689void
691 const Vector<Number> &x) const
692{
694}
695
696
697
698template <typename VectorType>
699void
701 const VectorType &x) const
702{
704}
705
706
707
708template <typename VectorType>
710 : BaseQR<VectorType>()
711{}
712
713
714
715template <typename VectorType>
716bool
717QR<VectorType>::append_column(const VectorType &column)
718{
719 // resize R:
720 this->R.grow_or_shrink(this->current_size + 1);
721 this->columns.push_back(std::make_unique<VectorType>(column));
722
723 // now a Gram-Schmidt part: orthonormalize the new column
724 // against everything we have so far:
725 auto &last_col = *this->columns.back();
726 for (unsigned int i = 0; i < this->current_size; ++i)
727 {
728 const auto &i_col = *this->columns[i];
729 this->R(i, this->current_size) = i_col * last_col;
730 last_col.add(-this->R(i, this->current_size), i_col);
731 }
732
733 this->R(this->current_size, this->current_size) = last_col.l2_norm();
734
735 Assert(this->R(this->current_size, this->current_size) > 0.,
737 last_col *= 1. / this->R(this->current_size, this->current_size);
738
739 ++this->current_size;
740 return true;
741}
742
743
744
745template <typename VectorType>
746void
747QR<VectorType>::apply_givens_rotation(const unsigned int i,
748 const unsigned int k)
749{
750 AssertIndexRange(i, k);
751 AssertIndexRange(k, this->current_size);
752 const std::array<Number, 3> csr =
753 ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
754 this->R(k, k));
755
756 // first, set k'th column:
757 this->R(i, k) = csr[2];
758 this->R(k, k) = 0.;
759 // now do the rest:
760 for (unsigned int j = 0; j < this->R.n(); ++j)
761 if (j != k)
762 {
763 const Number t = this->R(i, j);
764 this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
765 this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
766 }
767
768 // now adjust i,k columns due to multiplication with the
769 // transpose Givens matrix from right:
770 auto &col_i = *this->columns[i];
771 auto &col_k = *this->columns[k];
772 // save column i:
773 tmp = col_i;
774 col_i.sadd(csr[0], csr[1], col_k);
775 col_k.sadd(csr[0], -csr[1], tmp);
776
777 if (!this->givens_signal.empty())
778 this->givens_signal(i, k, csr);
779}
780
781
782
783template <typename VectorType>
784void
785QR<VectorType>::remove_column(const unsigned int k)
786{
787 AssertIndexRange(k, this->current_size);
788 Assert(this->current_size > 0,
789 ExcMessage("Can not remove a column if QR is empty"));
790 // apply a sequence of Givens rotations
791 // see section 6.5 "Updating matrix factorizations" in Golub 2013, Matrix
792 // computations
793
794 // So we want to have QR for \tilde A \in R^{m*(n-1)}
795 // if we remove the column k, we end up with upper Hessenberg matrix
796 // x x x x x
797 // x x x x
798 // H = x x x
799 // x x x
800 // x x
801 // x
802 // where k = 2 (3rd column), m = 7, n = 6
803 //
804 // before actually removing a column from Q and resizing R,
805 // apply givens rotations to bring H into upper triangular form:
806 for (unsigned int j = k + 1; j < this->R.n(); ++j)
807 {
808 const unsigned int i = j - 1;
809 apply_givens_rotation(i, j);
810 }
811
812 // now we can throw away the column from Q and adjust R
813 // since we do thin-QR, after Givens rotations we need to throw
814 // away the last column:
815 const unsigned int size_minus_1 = this->columns.size() - 1;
816 this->columns.erase(this->columns.begin() + size_minus_1);
817
818 // remove last row and k-th column
819 --this->current_size;
820 this->R.remove_row_and_column(this->current_size, k);
821}
822
823
824
825template <typename VectorType>
826void
827QR<VectorType>::multiply_with_Q(VectorType &y, const Vector<Number> &x) const
828{
830}
831
832
833
834template <typename VectorType>
835void
836QR<VectorType>::multiply_with_QT(Vector<Number> &y, const VectorType &x) const
837{
839}
840
841
842
843template <typename VectorType>
844void
845QR<VectorType>::multiply_with_A(VectorType &y, const Vector<Number> &x) const
846{
847 Vector<Number> x1 = x;
848 const int N = this->current_size;
849 const int lda = N;
850 const int incx = 1;
852 'U', 'N', 'N', N, &this->R(0, 0), lda, &x1[0], incx);
853
854 multiply_with_Q(y, x1);
855}
856
857
858
859template <typename VectorType>
860void
861QR<VectorType>::multiply_with_AT(Vector<Number> &y, const VectorType &x) const
862{
863 multiply_with_QT(y, x);
864
865 const int N = this->current_size;
866 const int lda = N;
867 const int incx = 1;
869 'U', 'T', 'N', N, &this->R(0, 0), lda, &y[0], incx);
870}
871
872#endif // no DOXYGEN
873
875
876#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:316
boost::signals2::signal< bool(const Vector< Number > &u, const Number &rho, const Number &col_norm_sqr)> column_signal
Definition qr.h:391
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:202
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const
typename VectorType::value_type Number
Definition qr.h:207
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:287
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
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 > &)