From: Guido Kanschat Date: Wed, 6 Apr 2005 18:01:07 +0000 (+0000) Subject: works! X-Git-Tag: v8.0.0~14155 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=73f4b6a93bd5a7026e196ee4a26438ee71a58a5d;p=dealii.git works! git-svn-id: https://svn.dealii.org/trunk@10392 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/householder.h b/deal.II/lac/include/lac/householder.h index 9669d0a217..361a1c12e2 100644 --- a/deal.II/lac/include/lac/householder.h +++ b/deal.II/lac/include/lac/householder.h @@ -90,26 +90,39 @@ Householder::Householder(const FullMatrix& M) { number2 sigma = 0; unsigned int i; - for (i=j ; iel(i,j)*this->el(i,j); + // sigma = ||v||^2 + for (i=j ; iel(i,j)*this->el(i,j); + // We are ready if the column is + // empty. Are we? if (std::fabs(sigma) < 1.e-15) return; - number2 s = this->el(j,j); - s = (s<0) ? std::sqrt(sigma) : -std::sqrt(sigma); - number2 dj = s; - - number2 beta = 1./(s*this->el(j,j)-sigma); - this->el(j,j) -= s; - for (unsigned int k=j+1 ; kel(i,j)*this->el(i,k); - sum *= beta; + number2 s = (this->el(j,j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma); + // + number2 beta = std::sqrt(1./(sigma-s*this->el(j,j))); + + // Make column j the Householder + // vector, store first entry in + // diagonal + diagonal[j] = beta*(this->el(j,j) - s); + this->el(j,j) = s; + + for (i=j+1 ; iel(i,j) *= beta; - for (i=j ; iel(i,k) += sum*this->el(i,j); - } - diagonal[j] = this->el(j,j); - this->el(j,j) = dj; + // For all subsequent columns do + // the Householder reflexion + for (unsigned int k=j+1 ; kel(j,k); + for (i=j+1 ; iel(i,j)*this->el(i,k); + + this->el(j,k) -= sum*this->diagonal[j]; + for (i=j+1 ; iel(i,k) -= sum*this->el(i,j); + } } } @@ -123,18 +136,19 @@ Householder::least_squares (Vector& dst, // Assert (!this->empty(), ExcEmptyMatrix()); // 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)*src(i); -// F*** what is beta??? -// sum *= beta; - - src(j) += sum*diagonal[j]; + // src -= v * sum + src(j) -= sum*diagonal[j]; for (unsigned int i=j+1 ; iel(i,j); + src(i) -= sum*this->el(i,j); } backward(dst, src);