From: Guido Kanschat Date: Wed, 6 Apr 2005 16:04:32 +0000 (+0000) Subject: File compiles, but results are nonsense. What the heck did I do there? X-Git-Tag: v8.0.0~14158 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fb71b5ef058de5e6ae4153d9371d349349328e39;p=dealii.git 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 --- 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 +#include + +#include + +// forward declarations +template class Vector; + + +/*! @addtogroup Matrix2 + *@{ + */ + + +/** + * QR-decomposition of a full matrix. + * + * @ref Instantiations: some (@ @) + * + * @author Guido Kanschat, 2005 + */ +template +class Householder : private FullMatrix +{ + public: + /** + * Create an object holding the + * QR-decomposition of a matrix. + */ + template + Householder (const FullMatrix&); + + /** + * Solve the least-squares + * problem for the right hand + * side src. The return + * value is the Euclidean norm of + * the approximation error. + */ + template + double least_squares (Vector &dst, + Vector &src); + + private: + /** + * Storage for the diagonal + * elements of the orthogonal + * transformation. + */ + std::vector diagonal; +}; + +/*@}*/ + +/// @if NoDoc +/*-------------------------Inline functions -------------------------------*/ + +// QR-transformation cf. Stoer 1 4.8.2 (p. 191) + +template +template +Householder::Householder(const FullMatrix& M) + : + FullMatrix(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 ; jel(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 ; kel(i,j)*this->el(i,k); + sum *= beta; + + for (i=j ; iel(i,k) += sum*this->el(i,j); + } + + diagonal[j] = this->el(j,j); + this->el(j,j) = dj; + } +} + + +template +template +double +Householder::least_squares (Vector& dst, + Vector& src) +{ +// Assert (!this->empty(), ExcEmptyMatrix()); + + // m > n, m = src.n, n = dst.n + + for (unsigned int j=0;jel(i,j)*src(i); +// F*** what is beta??? +// sum *= beta; + + src(j) += sum*diagonal[j]; + for (unsigned int i=j+1 ; iel(i,j); + } + + backward(dst, src); + + number2 sum = 0.; + for (unsigned int i=n() ; i