28 #include <boost/signals2/signal.hpp>
43 template <
typename VectorType>
49 using Number =
typename VectorType::value_type;
96 solve(Vector<Number> & x,
97 const Vector<Number> &y,
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);
160 std::vector<std::unique_ptr<VectorType>>
columns;
176 boost::signals2::signal<void(
const unsigned int i,
177 const unsigned int j,
178 const std::array<Number, 3> &)>
239 template <
typename VectorType>
246 using Number =
typename VectorType::value_type;
256 virtual ~QR() =
default;
348 template <
typename VectorType>
355 using Number =
typename VectorType::value_type;
406 boost::signals2::connection
408 const std::function<
bool(
const Vector<Number> &u,
410 const Number & col_norm_sqr)> &slot);
427 boost::signals2::signal<
bool(
const Vector<Number> &u,
429 const Number & col_norm_sqr)>
438 template <
typename VectorType>
447 template <
typename VectorType>
456 template <
typename VectorType>
465 template <
typename VectorType>
468 const Vector<Number> &y,
471 Assert(x.size() == this->current_size,
473 Assert(y.size() == this->current_size,
480 const int lda = this->current_size;
481 const int ldb = this->current_size;
482 const int N = this->current_size;
499 template <
typename VectorType>
502 const Vector<Number> &x)
const
504 Assert(x.size() == this->current_size,
508 for (
unsigned int j = 0; j < this->current_size; ++j)
509 y.add(x[j], *this->columns[j]);
514 template <
typename VectorType>
519 Assert(y.size() == this->current_size,
522 for (
unsigned int j = 0; j < this->current_size; ++j)
523 y[j] = (*this->columns[j]) * x;
528 template <
class VectorType>
529 boost::signals2::connection
531 const std::function<
void(
const unsigned int i,
532 const unsigned int j,
533 const std::array<Number, 3> &)> &slot)
535 return givens_signal.connect(slot);
540 template <
class VectorType>
541 boost::signals2::connection
543 const std::function<
bool(
const Vector<Number> &u,
545 const Number & col_norm_sqr)> &slot)
547 return column_signal.connect(slot);
552 template <
typename VectorType>
559 template <
typename VectorType>
563 if (this->current_size == 0)
565 this->R.grow_or_shrink(this->current_size + 1);
566 this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
567 this->R(0, 0) = column.l2_norm();
568 ++this->current_size;
573 Vector<Number> u(this->current_size);
574 this->multiply_with_AT(u, column);
577 const int lda = this->current_size;
578 const int ldb = this->current_size;
579 const int N = this->current_size;
583 "U",
"T",
"N", &
N, &n_rhs, &this->R(0, 0), &lda, &u(0), &ldb, &info);
587 const Number column_norm_sqr = column.norm_sqr();
588 const Number rho2 = column_norm_sqr - u.norm_sqr();
589 const bool linearly_independent =
590 column_signal.empty() ? rho2 > 0 :
591 column_signal(u, rho2, column_norm_sqr).get();
594 if (!linearly_independent)
599 this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
600 this->R.grow_or_shrink(this->current_size + 1);
601 this->R(this->current_size, this->current_size) =
std::sqrt(rho2);
602 for (
unsigned int i = 0; i < this->current_size; ++i)
603 this->R(i, this->current_size) = u(i);
605 this->current_size++;
613 template <
typename VectorType>
616 const unsigned int k)
620 const std::array<Number, 3> csr =
621 ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
625 this->R(i, k) = csr[2];
628 for (
unsigned int j = 0; j < this->R.n(); ++j)
631 const Number t = this->R(i, j);
632 this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
633 this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
636 if (!this->givens_signal.empty())
637 this->givens_signal(i, k, csr);
642 template <
typename VectorType>
648 for (
unsigned int j = k + 1; j < this->R.n(); ++j)
650 const unsigned int i = j - 1;
651 apply_givens_rotation(i, j);
655 --this->current_size;
656 this->R.remove_row_and_column(this->current_size, k);
659 this->columns.erase(this->columns.begin() + k);
664 template <
typename VectorType>
667 const Vector<Number> &x)
const
671 Vector<Number> x1 = x;
673 multiply_with_A(y, x1);
678 template <
typename VectorType>
686 multiply_with_AT(y, x);
692 template <
typename VectorType>
695 const Vector<Number> &x)
const
702 template <
typename VectorType>
712 template <
typename VectorType>
719 template <
typename VectorType>
724 this->R.grow_or_shrink(this->current_size + 1);
725 this->columns.push_back(std_cxx14::make_unique<VectorType>(column));
729 auto &last_col = *this->columns.back();
730 for (
unsigned int i = 0; i < this->current_size; ++i)
732 const auto &i_col = *this->columns[i];
733 this->R(i, this->current_size) = i_col * last_col;
734 last_col.add(-this->R(i, this->current_size), i_col);
737 this->R(this->current_size, this->current_size) = last_col.l2_norm();
739 Assert(this->R(this->current_size, this->current_size) > 0.,
741 last_col *= 1. / this->R(this->current_size, this->current_size);
743 ++this->current_size;
749 template <
typename VectorType>
752 const unsigned int k)
756 const std::array<Number, 3> csr =
757 ::Utilities::LinearAlgebra::givens_rotation<Number>(this->R(i, k),
761 this->R(i, k) = csr[2];
764 for (
unsigned int j = 0; j < this->R.n(); ++j)
767 const Number t = this->R(i, j);
768 this->R(i, j) = csr[0] * this->R(i, j) + csr[1] * this->R(k, j);
769 this->R(k, j) = -csr[1] * t + csr[0] * this->R(k, j);
774 auto &col_i = *this->columns[i];
775 auto &col_k = *this->columns[k];
778 col_i.sadd(csr[0], csr[1], col_k);
779 col_k.sadd(csr[0], -csr[1], tmp);
781 if (!this->givens_signal.empty())
782 this->givens_signal(i, k, csr);
787 template <
typename VectorType>
792 Assert(this->current_size > 0,
793 ExcMessage(
"Can not remove a column if QR is empty"));
810 for (
unsigned int j = k + 1; j < this->R.n(); ++j)
812 const unsigned int i = j - 1;
813 apply_givens_rotation(i, j);
819 const unsigned int size_minus_1 = this->columns.size() - 1;
820 this->columns.erase(this->columns.begin() + size_minus_1);
823 --this->current_size;
824 this->R.remove_row_and_column(this->current_size, k);
829 template <
typename VectorType>
838 template <
typename VectorType>
847 template <
typename VectorType>
851 Vector<Number> x1 = x;
852 const int N = this->current_size;
855 trmv(
"U",
"N",
"N", &
N, &this->R(0, 0), &lda, &x1[0], &incx);
857 multiply_with_Q(y, x1);
862 template <
typename VectorType>
866 multiply_with_QT(y, x);
868 const int N = this->current_size;
871 trmv(
"U",
"T",
"N", &
N, &this->R(0, 0), &lda, &y[0], &incx);