]> https://gitweb.dealii.org/ - dealii.git/commitdiff
copy the functions least_squares for BlockVectors
authorBaerbel Jannsen <baerbel.janssen@gmail.com>
Tue, 5 Jan 2010 10:02:05 +0000 (10:02 +0000)
committerBaerbel Jannsen <baerbel.janssen@gmail.com>
Tue, 5 Jan 2010 10:02:05 +0000 (10:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@20285 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 28833bf57966c39c2a750e62e3cfd8f82ec8c355..76ee424b66e5d9b7ffe84c1691b84d59b6e9a562 100644 (file)
@@ -93,6 +93,10 @@ class Householder : private FullMatrix<number>
     double least_squares (Vector<number2> &dst,
                          const Vector<number2> &src) const;
 
+    template<typename number2>
+    double least_squares (BlockVector<number2> &dst,
+                         const BlockVector<number2> &src) const;
+
   private:
                                     /**
                                      * Storage for the diagonal
@@ -222,6 +226,54 @@ Householder<number>::least_squares (Vector<number2>& dst,
   return std::sqrt(sum);
 }
 
+template <typename number>
+template <typename number2>
+double
+Householder<number>::least_squares (BlockVector<number2>& dst,
+                                   const BlockVector<number2>& src) const
+{
+  Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
+  AssertDimension(dst.size(), this->n());
+  AssertDimension(src.size(), this->m());
+
+  const unsigned int m = this->m(), n = this->n();
+  
+  GrowingVectorMemory<BlockVector<number2> > mem;
+  BlockVector<number2>* 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;j<n;++j)
+    {
+                                      // sum = v_i^T dst
+      number2 sum = diagonal[j]* (*aux)(j);
+      for (unsigned int i=j+1 ; i<m ; ++i)
+       sum += this->el(i,j)*(*aux)(i);
+                                      // dst -= v * sum
+      (*aux)(j) -= sum*diagonal[j];
+      for (unsigned int i=j+1 ; i<m ; ++i)
+       (*aux)(i) -= sum*this->el(i,j);
+    }
+                                  // Compute norm of residual
+  number2 sum = 0.;
+  for (unsigned int i=n ; i<m ; ++i)
+    sum += (*aux)(i) * (*aux)(i);
+                                  // Compute solution
+                                   //backward works for 
+                                   //Vectors only, so copy 
+                                   //them before
+  Vector<number2> 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

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.