25#include <boost/signals2/signal.hpp>
42template <
typename VectorType>
48 using Number =
typename VectorType::value_type;
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);
159 std::vector<std::unique_ptr<VectorType>>
columns;
175 boost::signals2::signal<void(
const unsigned int i,
176 const unsigned int j,
177 const std::array<Number, 3> &)>
238template <
typename VectorType>
245 using Number =
typename VectorType::value_type;
347template <
typename VectorType>
354 using Number =
typename VectorType::value_type;
405 boost::signals2::connection
409 const Number &col_norm_sqr)> &slot);
428 const Number &col_norm_sqr)>
437 namespace QRImplementation
444 template <
typename Number>
455 template <
typename Number>
472template <
typename VectorType>
481template <
typename VectorType>
490template <
typename VectorType>
499template <
typename VectorType>
514 const int lda = this->current_size;
515 const int ldb = this->current_size;
516 const int N = this->current_size;
533template <
typename VectorType>
542 for (
unsigned int j = 0; j < this->current_size; ++j)
543 y.add(x[j], *this->columns[j]);
548template <
typename VectorType>
551 const VectorType &x)
const
556 for (
unsigned int j = 0; j < this->current_size; ++j)
557 y[j] = (*this->columns[j]) * x;
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)
569 return givens_signal.connect(slot);
574template <
typename VectorType>
575boost::signals2::connection
579 const Number &col_norm_sqr)> &slot)
581 return column_signal.connect(slot);
586template <
typename VectorType>
593template <
typename VectorType>
597 if (this->current_size == 0)
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;
608 this->multiply_with_AT(u, column);
611 const int lda = this->current_size;
612 const int ldb = this->current_size;
613 const int N = this->current_size;
617 'U',
'T',
'N', N, n_rhs, &this->R(0, 0), lda, &u(0), ldb, &info);
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();
628 if (!linearly_independent)
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);
639 this->current_size++;
647template <
typename VectorType>
650 const unsigned int k)
654 const std::array<Number, 3> csr =
655 ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
659 this->R(i, k) = csr[2];
662 for (
unsigned int j = 0; j < this->R.n(); ++j)
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);
670 if (!this->givens_signal.empty())
671 this->givens_signal(i, k, csr);
676template <
typename VectorType>
682 for (
unsigned int j = k + 1; j < this->R.n(); ++j)
684 const unsigned int i = j - 1;
685 apply_givens_rotation(i, j);
689 --this->current_size;
690 this->R.remove_row_and_column(this->current_size, k);
693 this->columns.erase(this->columns.begin() + k);
698template <
typename VectorType>
707 multiply_with_A(y, x1);
712template <
typename VectorType>
715 const VectorType &x)
const
720 multiply_with_AT(y, x);
726template <
typename VectorType>
736template <
typename VectorType>
739 const VectorType &x)
const
746template <
typename VectorType>
753template <
typename VectorType>
758 this->R.grow_or_shrink(this->current_size + 1);
759 this->columns.push_back(std::make_unique<VectorType>(column));
763 auto &last_col = *this->columns.back();
764 for (
unsigned int i = 0; i < this->current_size; ++i)
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);
771 this->R(this->current_size, this->current_size) = last_col.l2_norm();
773 Assert(this->R(this->current_size, this->current_size) > 0.,
775 last_col *= 1. / this->R(this->current_size, this->current_size);
777 ++this->current_size;
783template <
typename VectorType>
786 const unsigned int k)
790 const std::array<Number, 3> csr =
791 ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
795 this->R(i, k) = csr[2];
798 for (
unsigned int j = 0; j < this->R.n(); ++j)
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);
808 auto &col_i = *this->columns[i];
809 auto &col_k = *this->columns[k];
812 col_i.sadd(csr[0], csr[1], col_k);
813 col_k.sadd(csr[0], -csr[1], tmp);
815 if (!this->givens_signal.empty())
816 this->givens_signal(i, k, csr);
821template <
typename VectorType>
826 Assert(this->current_size > 0,
827 ExcMessage(
"Can not remove a column if QR is empty"));
844 for (
unsigned int j = k + 1; j < this->R.n(); ++j)
846 const unsigned int i = j - 1;
847 apply_givens_rotation(i, j);
853 const unsigned int size_minus_1 = this->columns.size() - 1;
854 this->columns.erase(this->columns.begin() + size_minus_1);
857 --this->current_size;
858 this->R.remove_row_and_column(this->current_size, k);
863template <
typename VectorType>
872template <
typename VectorType>
881template <
typename VectorType>
886 const int N = this->current_size;
890 'U',
'N',
'N', N, &this->R(0, 0), lda, &x1[0], incx);
892 multiply_with_Q(y, x1);
897template <
typename VectorType>
901 multiply_with_QT(y, x);
903 const int N = this->current_size;
907 'U',
'T',
'N', N, &this->R(0, 0), lda, &y[0], incx);
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
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
typename VectorType::value_type Number
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
void multiply_with_cols(VectorType &y, const Vector< Number > &x) const
unsigned int current_size
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
boost::signals2::signal< bool(const Vector< Number > &u, const Number &rho, const Number &col_norm_sqr)> column_signal
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)
virtual void multiply_with_QT(Vector< Number > &y, const VectorType &x) const
typename VectorType::value_type Number
virtual bool append_column(const VectorType &column)
void apply_givens_rotation(const unsigned int i, const unsigned int k)
virtual void multiply_with_A(VectorType &y, const Vector< Number > &x) const
virtual void multiply_with_AT(Vector< Number > &y, const VectorType &x) const
virtual void remove_column(const unsigned int k=0)
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
#define DEAL_II_NAMESPACE_CLOSE
#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)
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)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)