include/deal.II/lac/householder.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: householder.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2005, 2006, 2007, 2008, 2010, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__householder_h
00013 #define __deal2__householder_h
00014 
00015 
00016 #include <cmath>
00017 #include <deal.II/base/config.h>
00018 #include <deal.II/lac/full_matrix.h>
00019 #include <deal.II/lac/vector_memory.h>
00020 
00021 #include <vector>
00022 
00023 DEAL_II_NAMESPACE_OPEN
00024 
00025 
00026 // forward declarations
00027 template<typename number> class Vector;
00028 
00029 
00050 template<typename number>
00051 class Householder : private FullMatrix<number>
00052 {
00053   public:
00057     Householder ();
00058 
00063     template<typename number2>
00064     Householder (const FullMatrix<number2>&);
00065 
00070     template<typename number2>
00071     void
00072     initialize (const FullMatrix<number2>&);
00073 
00091     template<typename number2>
00092     double least_squares (Vector<number2> &dst,
00093                           const Vector<number2> &src) const;
00094 
00099     template<typename number2>
00100     double least_squares (BlockVector<number2> &dst,
00101                           const BlockVector<number2> &src) const;
00102 
00108     template<class VECTOR>
00109     void vmult (VECTOR &dst, const VECTOR &src) const;
00110 
00111     template<class VECTOR>
00112     void Tvmult (VECTOR &dst, const VECTOR &src) const;
00113 
00114 
00115   private:
00121     std::vector<number> diagonal;
00122 };
00123 
00126 #ifndef DOXYGEN
00127 /*-------------------------Inline functions -------------------------------*/
00128 
00129 // QR-transformation cf. Stoer 1 4.8.2 (p. 191)
00130 
00131 template <typename number>
00132 Householder<number>::Householder()
00133 {}
00134 
00135 
00136 
00137 template <typename number>
00138 template <typename number2>
00139 void
00140 Householder<number>::initialize(const FullMatrix<number2>& M)
00141 {
00142   const unsigned int m = M.n_rows(), n = M.n_cols();
00143   this->reinit(m, n);
00144   this->fill(M);
00145   Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
00146   diagonal.resize(m);
00147 
00148                                    // m > n, src.n() = m
00149   Assert (this->n_cols() <= this->n_rows(),
00150           ExcDimensionMismatch(this->n_cols(), this->n_rows()));
00151 
00152   for (unsigned int j=0 ; j<n ; ++j)
00153   {
00154     number2 sigma = 0;
00155     unsigned int i;
00156                                      // sigma = ||v||^2
00157     for (i=j ; i<m ; ++i)
00158       sigma += this->el(i,j)*this->el(i,j);
00159                                      // We are ready if the column is
00160                                      // empty. Are we?
00161     if (std::fabs(sigma) < 1.e-15) return;
00162 
00163     number2 s = (this->el(j,j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
00164                                      //
00165     number2 beta = std::sqrt(1./(sigma-s*this->el(j,j)));
00166 
00167                                      // Make column j the Householder
00168                                      // vector, store first entry in
00169                                      // diagonal
00170     diagonal[j] = beta*(this->el(j,j) - s);
00171     this->el(j,j) = s;
00172 
00173     for (i=j+1 ; i<m ; ++i)
00174       this->el(i,j) *= beta;
00175 
00176 
00177                                      // For all subsequent columns do
00178                                      // the Householder reflexion
00179     for (unsigned int k=j+1 ; k<n ; ++k)
00180     {
00181       number2 sum = diagonal[j]*this->el(j,k);
00182       for (i=j+1 ; i<m ; ++i)
00183         sum += this->el(i,j)*this->el(i,k);
00184 
00185       this->el(j,k) -= sum*this->diagonal[j];
00186       for (i=j+1 ; i<m ; ++i)
00187         this->el(i,k) -= sum*this->el(i,j);
00188     }
00189   }
00190 }
00191 
00192 
00193 template <typename number>
00194 template <typename number2>
00195 Householder<number>::Householder(const FullMatrix<number2>& M)
00196 {
00197   initialize(M);
00198 }
00199 
00200 
00201 template <typename number>
00202 template <typename number2>
00203 double
00204 Householder<number>::least_squares (Vector<number2>& dst,
00205                                     const Vector<number2>& src) const
00206 {
00207   Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
00208   AssertDimension(dst.size(), this->n());
00209   AssertDimension(src.size(), this->m());
00210 
00211   const unsigned int m = this->m(), n = this->n();
00212 
00213   GrowingVectorMemory<Vector<number2> > mem;
00214   Vector<number2>* aux = mem.alloc();
00215   aux->reinit(src, true);
00216   *aux = src;
00217                                    // m > n, m = src.n, n = dst.n
00218 
00219                                    // Multiply Q_n ... Q_2 Q_1 src
00220                                    // Where Q_i = I-v_i v_i^T
00221   for (unsigned int j=0;j<n;++j)
00222     {
00223                                        // sum = v_i^T dst
00224       number2 sum = diagonal[j]* (*aux)(j);
00225       for (unsigned int i=j+1 ; i<m ; ++i)
00226         sum += this->el(i,j)*(*aux)(i);
00227                                        // dst -= v * sum
00228       (*aux)(j) -= sum*diagonal[j];
00229       for (unsigned int i=j+1 ; i<m ; ++i)
00230         (*aux)(i) -= sum*this->el(i,j);
00231     }
00232                                    // Compute norm of residual
00233   number2 sum = 0.;
00234   for (unsigned int i=n ; i<m ; ++i)
00235     sum += (*aux)(i) * (*aux)(i);
00236                                    // Compute solution
00237   this->backward(dst, *aux);
00238 
00239   mem.free(aux);
00240 
00241   return std::sqrt(sum);
00242 }
00243 
00244 template <typename number>
00245 template <typename number2>
00246 double
00247 Householder<number>::least_squares (BlockVector<number2>& dst,
00248                                     const BlockVector<number2>& src) const
00249 {
00250   Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
00251   AssertDimension(dst.size(), this->n());
00252   AssertDimension(src.size(), this->m());
00253 
00254   const unsigned int m = this->m(), n = this->n();
00255 
00256   GrowingVectorMemory<BlockVector<number2> > mem;
00257   BlockVector<number2>* aux = mem.alloc();
00258   aux->reinit(src, true);
00259   *aux = src;
00260                                    // m > n, m = src.n, n = dst.n
00261 
00262                                    // Multiply Q_n ... Q_2 Q_1 src
00263                                    // Where Q_i = I-v_i v_i^T
00264   for (unsigned int j=0;j<n;++j)
00265     {
00266                                        // sum = v_i^T dst
00267       number2 sum = diagonal[j]* (*aux)(j);
00268       for (unsigned int i=j+1 ; i<m ; ++i)
00269         sum += this->el(i,j)*(*aux)(i);
00270                                        // dst -= v * sum
00271       (*aux)(j) -= sum*diagonal[j];
00272       for (unsigned int i=j+1 ; i<m ; ++i)
00273         (*aux)(i) -= sum*this->el(i,j);
00274     }
00275                                    // Compute norm of residual
00276   number2 sum = 0.;
00277   for (unsigned int i=n ; i<m ; ++i)
00278     sum += (*aux)(i) * (*aux)(i);
00279                                    //backward works for
00280                                    //Vectors only, so copy
00281                                    //them before
00282   Vector<number2> v_dst, v_aux;
00283   v_dst = dst;
00284   v_aux = *aux;
00285                                    // Compute solution
00286   this->backward(v_dst, v_aux);
00287                                    //copy the result back
00288                                    //to the BlockVector
00289   dst = v_dst;
00290 
00291   mem.free(aux);
00292 
00293   return std::sqrt(sum);
00294 }
00295 
00296 
00297 template <typename number>
00298 template <class VECTOR>
00299 void
00300 Householder<number>::vmult (VECTOR &dst, const VECTOR &src) const
00301 {
00302   least_squares (dst, src);
00303 }
00304 
00305 
00306 template <typename number>
00307 template <class VECTOR>
00308 void
00309 Householder<number>::Tvmult (VECTOR &, const VECTOR &) const
00310 {
00311   Assert(false, ExcNotImplemented());
00312 }
00313 
00314 
00315 
00316 #endif // DOXYGEN
00317 
00318 DEAL_II_NAMESPACE_CLOSE
00319 
00320 #endif
00321 
00322 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Tue May 22 2012 12:06:11 by doxygen 1.7.3