]> https://gitweb.dealii.org/ - dealii.git/commitdiff
File compiles, but results are nonsense. What the heck did I do there?
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 6 Apr 2005 16:04:32 +0000 (16:04 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 6 Apr 2005 16:04:32 +0000 (16:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@10389 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/householder.h [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/householder.h b/deal.II/lac/include/lac/householder.h
new file mode 100644 (file)
index 0000000..9669d0a
--- /dev/null
@@ -0,0 +1,152 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__householder_h
+#define __deal2__householder_h
+
+
+#include <base/config.h>
+#include <lac/full_matrix.h>
+
+#include <vector>
+
+// forward declarations
+template<typename number> class Vector;
+
+
+/*! @addtogroup Matrix2
+ *@{
+ */
+
+
+/**
+ * QR-decomposition of a full matrix.
+ *
+ * @ref Instantiations: some (<tt>@<float@> @<double@></tt>)
+ *
+ * @author Guido Kanschat, 2005
+ */
+template<typename number>
+class Householder : private FullMatrix<number>
+{
+  public:
+                                    /**
+                                     * Create an object holding the
+                                     * QR-decomposition of a matrix.
+                                     */
+    template<typename number2>
+    Householder (const FullMatrix<number2>&);
+
+                                    /**
+                                     * Solve the least-squares
+                                     * problem for the right hand
+                                     * side <tt>src</tt>.  The return
+                                     * value is the Euclidean norm of
+                                     * the approximation error.
+                                     */
+    template<typename number2>
+    double least_squares (Vector<number2> &dst,
+                         Vector<number2> &src);
+
+  private:
+                                    /**
+                                     * Storage for the diagonal
+                                     * elements of the orthogonal
+                                     * transformation.
+                                     */
+    std::vector<number> diagonal;
+};
+
+/*@}*/
+
+/// @if NoDoc
+/*-------------------------Inline functions -------------------------------*/
+
+// QR-transformation cf. Stoer 1 4.8.2 (p. 191)
+
+template <typename number>
+template <typename number2>
+Householder<number>::Householder(const FullMatrix<number2>& M)
+               :
+               FullMatrix<number>(M),
+               diagonal(M.n_rows())
+{
+//  Assert (!this->empty(), ExcEmptyMatrix());
+  
+                                  // m > n, src.n() = m
+  Assert (this->n_cols() <= this->n_rows(),
+         ExcDimensionMismatch(this->n_cols(), this->n_rows()));
+
+  for (unsigned int j=0 ; j<n() ; ++j)
+  {
+    number2 sigma = 0;
+    unsigned int i;
+    for (i=j ; i<m() ; ++i) sigma += this->el(i,j)*this->el(i,j);
+    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;
+
+      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;    
+  }
+}
+
+
+template <typename number>
+template <typename number2>
+double
+Householder<number>::least_squares (Vector<number2>& dst,
+                                   Vector<number2>& src)
+{
+//  Assert (!this->empty(), ExcEmptyMatrix());
+  
+                                  // m > n, m = src.n, n = dst.n
+  
+  for (unsigned int j=0;j<n();++j)
+    {
+      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];
+      for (unsigned int i=j+1 ; i<m() ; ++i)
+       src(i) += sum*this->el(i,j);
+    }
+  
+  backward(dst, src);
+  
+  number2 sum = 0.;
+  for (unsigned int i=n() ; i<m() ; ++i) sum += src(i) * src(i);
+  return std::sqrt(sum);
+}
+
+
+
+///@endif
+
+#endif
+

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.