From 0171282f3b87b2c86b1822d693a71297755fad5c Mon Sep 17 00:00:00 2001 From: Baerbel Jannsen Date: Tue, 5 Jan 2010 10:02:05 +0000 Subject: [PATCH] copy the functions least_squares for BlockVectors git-svn-id: https://svn.dealii.org/trunk@20285 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/householder.h | 52 +++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/deal.II/lac/include/lac/householder.h b/deal.II/lac/include/lac/householder.h index 28833bf579..76ee424b66 100644 --- a/deal.II/lac/include/lac/householder.h +++ b/deal.II/lac/include/lac/householder.h @@ -93,6 +93,10 @@ class Householder : private FullMatrix double least_squares (Vector &dst, const Vector &src) const; + template + double least_squares (BlockVector &dst, + const BlockVector &src) const; + private: /** * Storage for the diagonal @@ -222,6 +226,54 @@ Householder::least_squares (Vector& dst, return std::sqrt(sum); } +template +template +double +Householder::least_squares (BlockVector& dst, + const BlockVector& src) const +{ + Assert (!this->empty(), typename FullMatrix::ExcEmptyMatrix()); + AssertDimension(dst.size(), this->n()); + AssertDimension(src.size(), this->m()); + + const unsigned int m = this->m(), n = this->n(); + + GrowingVectorMemory > mem; + BlockVector* aux = mem.alloc(); + aux->reinit(src, true); + *aux = src; + // m > n, m = src.n, n = dst.n + + // Multiply Q_n ... Q_2 Q_1 src + // Where Q_i = I-v_i v_i^T + for (unsigned int j=0;jel(i,j)*(*aux)(i); + // dst -= v * sum + (*aux)(j) -= sum*diagonal[j]; + for (unsigned int i=j+1 ; iel(i,j); + } + // Compute norm of residual + number2 sum = 0.; + for (unsigned int i=n ; i v_dst, v_aux; + v_dst = dst; + v_aux = *aux; + this->backward(v_dst, v_aux); + + mem.free(aux); + + return std::sqrt(sum); +} #endif // DOXYGEN -- 2.39.5