From: kanschat Date: Mon, 5 May 2008 22:36:28 +0000 (+0000) Subject: use an auxiliary vector in Householder, such that rectangular systems can be treated X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=eb99c094146b848b441fce007f088cba7eb73b04;p=dealii-svn.git use an auxiliary vector in Householder, such that rectangular systems can be treated git-svn-id: https://svn.dealii.org/trunk@16031 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/householder.h b/deal.II/lac/include/lac/householder.h index 943548816d..3162e71f39 100644 --- a/deal.II/lac/include/lac/householder.h +++ b/deal.II/lac/include/lac/householder.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 2005, 2006, 2007, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -17,6 +17,7 @@ #include #include #include +#include #include @@ -182,7 +183,13 @@ Householder::least_squares (Vector& dst, const Vector& src) const { Assert (!this->empty(), typename FullMatrix::ExcEmptyMatrix()); - dst = src; + AssertDimension(dst.size(), this->n()); + AssertDimension(src.size(), this->m()); + + GrowingVectorMemory > mem; + Vector* 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 @@ -190,20 +197,22 @@ Householder::least_squares (Vector& dst, for (unsigned int j=0;jn();++j) { // sum = v_i^T dst - number2 sum = diagonal[j]*dst(j); + number2 sum = diagonal[j]* (*aux)(j); for (unsigned int i=j+1 ; im() ; ++i) - sum += this->el(i,j)*dst(i); + sum += this->el(i,j)*(*aux)(i); // dst -= v * sum - dst(j) -= sum*diagonal[j]; + (*aux)(j) -= sum*diagonal[j]; for (unsigned int i=j+1 ; im() ; ++i) - dst(i) -= sum*this->el(i,j); + (*aux)(i) -= sum*this->el(i,j); } // Compute norm of residual number2 sum = 0.; for (unsigned int i=this->n() ; im() ; ++i) - sum += dst(i) * dst(i); + sum += (*aux)(i) * (*aux)(i); // Compute solution - this->backward(dst, dst); + this->backward(dst, *aux); + + mem.free(aux); return std::sqrt(sum); }