]> https://gitweb.dealii.org/ - dealii.git/commitdiff
works!
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 6 Apr 2005 18:01:07 +0000 (18:01 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 6 Apr 2005 18:01:07 +0000 (18:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@10392 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/householder.h

index 9669d0a217960dea48651532964452308ad41d9c..361a1c12e29ebc1b14bc00a8a169404cd78371f5 100644 (file)
@@ -90,26 +90,39 @@ Householder<number>::Householder(const FullMatrix<number2>& M)
   {
     number2 sigma = 0;
     unsigned int i;
-    for (i=j ; i<m() ; ++i) sigma += this->el(i,j)*this->el(i,j);
+                                    // sigma = ||v||^2
+    for (i=j ; i<m() ; ++i)
+      sigma += this->el(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 ; k<n() ; ++k)
-    {
-      number2 sum = 0.;
-      for (i=j ; i<m() ; ++i) sum += this->el(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 ; i<m() ; ++i)
+      this->el(i,j) *= beta;
 
-      for (i=j ; i<m() ; ++i) this->el(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 ; k<n() ; ++k)
+    {
+      number2 sum = diagonal[j]*this->el(j,k);
+      for (i=j+1 ; i<m() ; ++i)
+       sum += this->el(i,j)*this->el(i,k);
+
+      this->el(j,k) -= sum*this->diagonal[j];
+      for (i=j+1 ; i<m() ; ++i)
+       this->el(i,k) -= sum*this->el(i,j);
+    }    
   }
 }
 
@@ -123,18 +136,19 @@ Householder<number>::least_squares (Vector<number2>& 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;j<n();++j)
     {
+                                      // sum = v_i^T src
       number2 sum = diagonal[j]*src(j);
       for (unsigned int i=j+1 ; i<m() ; ++i)
        sum += this->el(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 ; i<m() ; ++i)
-       src(i) += sum*this->el(i,j);
+       src(i) -= sum*this->el(i,j);
     }
   
   backward(dst, src);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.