From fb71b5ef058de5e6ae4153d9371d349349328e39 Mon Sep 17 00:00:00 2001
From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Wed, 6 Apr 2005 16:04:32 +0000
Subject: [PATCH] File compiles, but results are nonsense. What the heck did I
 do there?

git-svn-id: https://svn.dealii.org/trunk@10389 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/lac/include/lac/householder.h | 152 ++++++++++++++++++++++++++
 1 file changed, 152 insertions(+)
 create mode 100644 deal.II/lac/include/lac/householder.h

diff --git a/deal.II/lac/include/lac/householder.h b/deal.II/lac/include/lac/householder.h
new file mode 100644
index 0000000000..9669d0a217
--- /dev/null
+++ b/deal.II/lac/include/lac/householder.h
@@ -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
+
-- 
2.39.5