include/deal.II/lac/solver_gmres.h

00001 //----------------------------  solver_gmres.h  ---------------------------
00002 //    @f$Id: solver_gmres.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 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 //----------------------------  solver_gmres.h  ---------------------------
00012 #ifndef __deal2__solver_gmres_h
00013 #define __deal2__solver_gmres_h
00014 
00015 
00016 
00017 #include <deal.II/base/config.h>
00018 #include <deal.II/base/subscriptor.h>
00019 #include <deal.II/base/logstream.h>
00020 #include <deal.II/lac/householder.h>
00021 #include <deal.II/lac/solver.h>
00022 #include <deal.II/lac/solver_control.h>
00023 #include <deal.II/lac/full_matrix.h>
00024 #include <deal.II/lac/vector.h>
00025 
00026 #include <vector>
00027 #include <cmath>
00028 
00029 DEAL_II_NAMESPACE_OPEN
00030 
00033 
00034 namespace internal
00035 {
00040   namespace SolverGMRES
00041   {
00054     template <class VECTOR>
00055     class TmpVectors
00056     {
00057       public:
00063         TmpVectors(const unsigned int    max_size,
00064                    VectorMemory<VECTOR> &vmem);
00065 
00069         ~TmpVectors();
00070 
00077         VECTOR& operator[] (const unsigned int i) const;
00078 
00089         VECTOR& operator() (const unsigned int i,
00090                             const VECTOR      &temp);
00091 
00092       private:
00097         VectorMemory<VECTOR> &mem;
00098 
00103         std::vector<VECTOR*> data;
00104 
00111         unsigned int offset;
00112     };
00113   }
00114 }
00115 
00167 template <class VECTOR = Vector<double> >
00168 class SolverGMRES : public Solver<VECTOR>
00169 {
00170   public:
00176     struct AdditionalData
00177     {
00187         AdditionalData (const unsigned int max_n_tmp_vectors = 30,
00188                         const bool right_preconditioning = false,
00189                         const bool use_default_residual = true);
00190 
00200         unsigned int    max_n_tmp_vectors;
00201 
00214         bool right_preconditioning;
00215 
00221         bool use_default_residual;
00222     };
00223 
00227     SolverGMRES (SolverControl        &cn,
00228                  VectorMemory<VECTOR> &mem,
00229                  const AdditionalData &data=AdditionalData());
00230 
00236     SolverGMRES (SolverControl        &cn,
00237                  const AdditionalData &data=AdditionalData());
00238 
00243     template<class MATRIX, class PRECONDITIONER>
00244     void
00245     solve (const MATRIX         &A,
00246            VECTOR               &x,
00247            const VECTOR         &b,
00248            const PRECONDITIONER &precondition);
00249 
00250     DeclException1 (ExcTooFewTmpVectors,
00251                     int,
00252                     << "The number of temporary vectors you gave ("
00253                     << arg1 << ") is too small. It should be at least 10 for "
00254                     << "any results, and much more for reasonable ones.");
00255 
00256   protected:
00261     AdditionalData additional_data;
00262 
00267     virtual double criterion();
00268 
00275     void givens_rotation (Vector<double> &h,  Vector<double> &b,
00276                           Vector<double> &ci, Vector<double> &si,
00277                           int col) const;
00281     FullMatrix<double> H;
00285     FullMatrix<double> H1;
00286 
00287   private:
00291     SolverGMRES (const SolverGMRES<VECTOR>&);
00292 };
00293 
00313 template <class VECTOR = Vector<double> >
00314 class SolverFGMRES : public Solver<VECTOR>
00315 {
00316   public:
00322     struct AdditionalData
00323     {
00334         AdditionalData(const unsigned int max_basis_size = 30,
00335                        const bool /*use_default_residual*/ = true)
00336                         :
00337                         max_basis_size(max_basis_size)
00338           {}
00339 
00344         unsigned int    max_basis_size;
00345     };
00346 
00350     SolverFGMRES (SolverControl        &cn,
00351                   VectorMemory<VECTOR> &mem,
00352                   const AdditionalData &data=AdditionalData());
00353 
00359     SolverFGMRES (SolverControl        &cn,
00360                   const AdditionalData &data=AdditionalData());
00361 
00366     template<class MATRIX, class PRECONDITIONER>
00367     void
00368     solve (const MATRIX         &A,
00369            VECTOR               &x,
00370            const VECTOR         &b,
00371            const PRECONDITIONER &precondition);
00372 
00373   private:
00377     AdditionalData additional_data;
00381     FullMatrix<double> H;
00385     FullMatrix<double> H1;
00386 };
00387 
00389 /* --------------------- Inline and template functions ------------------- */
00390 
00391 
00392 #ifndef DOXYGEN
00393 namespace internal
00394 {
00395   namespace SolverGMRES
00396   {
00397     template <class VECTOR>
00398     inline
00399     TmpVectors<VECTOR>::
00400     TmpVectors (const unsigned int    max_size,
00401                 VectorMemory<VECTOR> &vmem)
00402                     :
00403                     mem(vmem),
00404                     data (max_size, 0),
00405                     offset(0)
00406     {}
00407 
00408 
00409     template <class VECTOR>
00410     inline
00411     TmpVectors<VECTOR>::~TmpVectors ()
00412     {
00413       for (typename std::vector<VECTOR*>::iterator v = data.begin();
00414            v != data.end(); ++v)
00415         if (*v != 0)
00416           mem.free(*v);
00417     }
00418 
00419 
00420     template <class VECTOR>
00421     inline VECTOR&
00422     TmpVectors<VECTOR>::operator[] (const unsigned int i) const
00423     {
00424       Assert (i+offset<data.size(),
00425               ExcIndexRange(i, -offset, data.size()-offset));
00426 
00427       Assert (data[i-offset] != 0, ExcNotInitialized());
00428       return *data[i-offset];
00429     }
00430 
00431 
00432     template <class VECTOR>
00433     inline VECTOR&
00434     TmpVectors<VECTOR>::operator() (const unsigned int i,
00435                                     const VECTOR      &temp)
00436     {
00437       Assert (i+offset<data.size(),
00438               ExcIndexRange(i,-offset, data.size()-offset));
00439       if (data[i-offset] == 0)
00440         {
00441           data[i-offset] = mem.alloc();
00442           data[i-offset]->reinit(temp);
00443         }
00444       return *data[i-offset];
00445     }
00446   }
00447 }
00448 
00449 
00450 
00451 template <class VECTOR>
00452 inline
00453 SolverGMRES<VECTOR>::AdditionalData::
00454 AdditionalData (const unsigned int max_n_tmp_vectors,
00455                 const bool         right_preconditioning,
00456                 const bool         use_default_residual)
00457                 :
00458                 max_n_tmp_vectors(max_n_tmp_vectors),
00459                 right_preconditioning(right_preconditioning),
00460                 use_default_residual(use_default_residual)
00461 {}
00462 
00463 
00464 template <class VECTOR>
00465 SolverGMRES<VECTOR>::SolverGMRES (SolverControl        &cn,
00466                                   VectorMemory<VECTOR> &mem,
00467                                   const AdditionalData &data)
00468                 :
00469                 Solver<VECTOR> (cn,mem),
00470                 additional_data(data)
00471 {}
00472 
00473 
00474 
00475 template <class VECTOR>
00476 SolverGMRES<VECTOR>::SolverGMRES (SolverControl        &cn,
00477                                   const AdditionalData &data) :
00478                 Solver<VECTOR> (cn),
00479                 additional_data(data)
00480 {}
00481 
00482 
00483 
00484 template <class VECTOR>
00485 inline
00486 void
00487 SolverGMRES<VECTOR>::givens_rotation (Vector<double> &h,
00488                                       Vector<double> &b,
00489                                       Vector<double> &ci,
00490                                       Vector<double> &si,
00491                                       int     col) const
00492 {
00493   for (int i=0 ; i<col ; i++)
00494     {
00495       const double s = si(i);
00496       const double c = ci(i);
00497       const double dummy = h(i);
00498       h(i)   =  c*dummy + s*h(i+1);
00499       h(i+1) = -s*dummy + c*h(i+1);
00500     };
00501 
00502   const double r = 1./std::sqrt(h(col)*h(col) + h(col+1)*h(col+1));
00503   si(col) = h(col+1) *r;
00504   ci(col) = h(col)   *r;
00505   h(col)  =  ci(col)*h(col) + si(col)*h(col+1);
00506   b(col+1)= -si(col)*b(col);
00507   b(col) *=  ci(col);
00508 }
00509 
00510 
00511 
00512 template<class VECTOR>
00513 template<class MATRIX, class PRECONDITIONER>
00514 void
00515 SolverGMRES<VECTOR>::solve (const MATRIX         &A,
00516                             VECTOR               &x,
00517                             const VECTOR         &b,
00518                             const PRECONDITIONER &precondition)
00519 {
00520                                    // this code was written a very
00521                                    // long time ago by people not
00522                                    // associated with deal.II. we
00523                                    // don't make any guarantees to its
00524                                    // optimality or that it even works
00525                                    // as expected...
00526 
00527 //TODO:[?] Check, why there are two different start residuals.
00528 //TODO:[GK] Make sure the parameter in the constructor means maximum basis size
00529 
00530   deallog.push("GMRES");
00531   const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors;
00532 
00533                                    // Generate an object where basis
00534                                    // vectors are stored.
00535   internal::SolverGMRES::TmpVectors<VECTOR> tmp_vectors (n_tmp_vectors, this->memory);
00536 
00537                                    // number of the present iteration; this
00538                                    // number is not reset to zero upon a
00539                                    // restart
00540   unsigned int accumulated_iterations = 0;
00541 
00542                                    // matrix used for the orthogonalization
00543                                    // process later
00544   H.reinit(n_tmp_vectors, n_tmp_vectors-1);
00545 
00546                                    // some additional vectors, also used
00547                                    // in the orthogonalization
00548   ::Vector<double>
00549     gamma(n_tmp_vectors),
00550     ci   (n_tmp_vectors-1),
00551     si   (n_tmp_vectors-1),
00552     h    (n_tmp_vectors-1);
00553 
00554 
00555   unsigned int dim = 0;
00556 
00557   SolverControl::State iteration_state = SolverControl::iterate;
00558 
00559                                    // switch to determine whether we want a
00560                                    // left or a right preconditioner. at
00561                                    // present, left is default, but both
00562                                    // ways are implemented
00563   const bool left_precondition = !additional_data.right_preconditioning;
00564                                    // Per default the left
00565                                    // preconditioned GMRes uses the
00566                                    // preconditioned residual and the
00567                                    // right preconditioned GMRes uses
00568                                    // the unpreconditioned residual as
00569                                    // stopping criterion.
00570   const bool use_default_residual = additional_data.use_default_residual;
00571 
00572                                    // define two aliases
00573   VECTOR &v = tmp_vectors(0, x);
00574   VECTOR &p = tmp_vectors(n_tmp_vectors-1, x);
00575 
00576                                    // Following vectors are needed
00577                                    // when not the default residuals
00578                                    // are used as stopping criterion
00579   VECTOR *r=0;
00580   VECTOR *x_=0;
00581   ::Vector<double> *gamma_=0;
00582   if (!use_default_residual)
00583     {
00584       r=this->memory.alloc();
00585       x_=this->memory.alloc();
00586       r->reinit(x);
00587       x_->reinit(x);
00588 
00589       gamma_ = new ::Vector<double> (gamma.size());
00590     }
00591 
00593                                    // outer iteration: loop until we
00594                                    // either reach convergence or the
00595                                    // maximum number of iterations is
00596                                    // exceeded. each cycle of this
00597                                    // loop amounts to one restart
00598   do
00599     {
00600                                        // reset this vector to the
00601                                        // right size
00602       h.reinit (n_tmp_vectors-1);
00603 
00604       if (left_precondition)
00605         {
00606           A.vmult(p,x);
00607           p.sadd(-1.,1.,b);
00608           precondition.vmult(v,p);
00609         }
00610       else
00611         {
00612           A.vmult(v,x);
00613           v.sadd(-1.,1.,b);
00614         };
00615 
00616       double rho = v.l2_norm();
00617 
00618                                        // check the residual here as
00619                                        // well since it may be that we
00620                                        // got the exact (or an almost
00621                                        // exact) solution vector at
00622                                        // the outset. if we wouldn't
00623                                        // check here, the next scaling
00624                                        // operation would produce
00625                                        // garbage
00626       if (use_default_residual)
00627         {
00628           iteration_state = this->control().check (
00629             accumulated_iterations, rho);
00630 
00631           if (iteration_state != SolverControl::iterate)
00632             break;
00633         }
00634       else
00635         {
00636           deallog << "default_res=" << rho << std::endl;
00637 
00638           if (left_precondition)
00639             {
00640               A.vmult(*r,x);
00641               r->sadd(-1.,1.,b);
00642             }
00643           else
00644             precondition.vmult(*r,v);
00645 
00646           double res = r->l2_norm();
00647           iteration_state = this->control().check (
00648             accumulated_iterations, res);
00649 
00650           if (iteration_state != SolverControl::iterate)
00651             {
00652               this->memory.free(r);
00653               this->memory.free(x_);
00654 
00655               delete gamma_;
00656               break;
00657             }
00658         }
00659 
00660       gamma(0) = rho;
00661 
00662       v *= 1./rho;
00663 
00664                                        // inner iteration doing at
00665                                        // most as many steps as there
00666                                        // are temporary vectors. the
00667                                        // number of steps actually
00668                                        // been done is propagated
00669                                        // outside through the @p dim
00670                                        // variable
00671       for (unsigned int inner_iteration=0;
00672            ((inner_iteration < n_tmp_vectors-2)
00673             &&
00674             (iteration_state==SolverControl::iterate));
00675            ++inner_iteration)
00676         {
00677           ++accumulated_iterations;
00678           // yet another alias
00679           VECTOR& vv = tmp_vectors(inner_iteration+1, x);
00680 
00681           if (left_precondition)
00682             {
00683               A.vmult(p, tmp_vectors[inner_iteration]);
00684               precondition.vmult(vv,p);
00685             } else {
00686               precondition.vmult(p, tmp_vectors[inner_iteration]);
00687               A.vmult(vv,p);
00688             };
00689 
00690           dim = inner_iteration+1;
00691 
00692                                            /* Orthogonalization */
00693           for (unsigned int i=0 ; i<dim ; ++i)
00694             {
00695               h(i) = vv * tmp_vectors[i];
00696               vv.add(-h(i), tmp_vectors[i]);
00697             };
00698 
00699                                            /* Re-orthogonalization */
00700           for (unsigned int i=0 ; i<dim ; ++i)
00701             {
00702               double htmp = vv * tmp_vectors[i];
00703               h(i) += htmp;
00704               vv.add(-htmp, tmp_vectors[i]);
00705             }
00706 
00707           const double s = vv.l2_norm();
00708           h(inner_iteration+1) = s;
00709 //TODO: s=0 is a lucky breakdown. Handle this somehow decently
00710 
00711           vv *= 1./s;
00712 
00713                                            /*  Transformation into
00714                                                triagonal structure  */
00715           givens_rotation(h,gamma,ci,si,inner_iteration);
00716 
00717                                            /*  append vector on matrix  */
00718           for (unsigned int i=0; i<dim; ++i)
00719             H(i,inner_iteration) = h(i);
00720 
00721                                            /*  default residual  */
00722           rho = std::fabs(gamma(dim));
00723 
00724           if (use_default_residual)
00725             iteration_state = this->control().check (
00726               accumulated_iterations, rho);
00727           else
00728             {
00729               deallog << "default_res=" << rho << std::endl;
00730 
00731               ::Vector<double> h_(dim);
00732               *x_=x;
00733               *gamma_=gamma;
00734               H1.reinit(dim+1,dim);
00735 
00736               for (unsigned int i=0; i<dim+1; ++i)
00737                 for (unsigned int j=0; j<dim; ++j)
00738                   H1(i,j) = H(i,j);
00739 
00740               H1.backward(h_,*gamma_);
00741 
00742               if (left_precondition)
00743                 for (unsigned int i=0 ; i<dim; ++i)
00744                   x_->add(h_(i), tmp_vectors[i]);
00745               else
00746                 {
00747                   p = 0.;
00748                   for (unsigned int i=0; i<dim; ++i)
00749                     p.add(h_(i), tmp_vectors[i]);
00750                   precondition.vmult(*r,p);
00751                   x_->add(1.,*r);
00752                 };
00753               A.vmult(*r,*x_);
00754               r->sadd(-1.,1.,b);
00755                                                // Now *r contains the
00756                                                // unpreconditioned
00757                                                // residual!!
00758               if (left_precondition)
00759                 {
00760                   const double res=r->l2_norm();
00761 
00762                   iteration_state = this->control().check (
00763                     accumulated_iterations, res);
00764                 }
00765               else
00766                 {
00767                   precondition.vmult(*x_, *r);
00768                   const double preconditioned_res=x_->l2_norm();
00769 
00770                   iteration_state = this->control().check (
00771                     accumulated_iterations, preconditioned_res);
00772                 }
00773             }
00774         };
00775                                        // end of inner iteration. now
00776                                        // calculate the solution from
00777                                        // the temporary vectors
00778       h.reinit(dim);
00779       H1.reinit(dim+1,dim);
00780 
00781       for (unsigned int i=0; i<dim+1; ++i)
00782         for (unsigned int j=0; j<dim; ++j)
00783           H1(i,j) = H(i,j);
00784 
00785       H1.backward(h,gamma);
00786 
00787       if (left_precondition)
00788         for (unsigned int i=0 ; i<dim; ++i)
00789           x.add(h(i), tmp_vectors[i]);
00790       else
00791         {
00792           p = 0.;
00793           for (unsigned int i=0; i<dim; ++i)
00794             p.add(h(i), tmp_vectors[i]);
00795           precondition.vmult(v,p);
00796           x.add(1.,v);
00797         };
00798                                        // end of outer iteration. restart if
00799                                        // no convergence and the number of
00800                                        // iterations is not exceeded
00801     }
00802   while (iteration_state == SolverControl::iterate);
00803 
00804   if (!use_default_residual)
00805     {
00806       this->memory.free(r);
00807       this->memory.free(x_);
00808 
00809       delete gamma_;
00810     }
00811 
00812   deallog.pop();
00813                                    // in case of failure: throw
00814                                    // exception
00815   if (this->control().last_check() != SolverControl::success)
00816     throw SolverControl::NoConvergence (this->control().last_step(),
00817                                         this->control().last_value());
00818                                    // otherwise exit as normal
00819 }
00820 
00821 
00822 
00823 template<class VECTOR>
00824 double
00825 SolverGMRES<VECTOR>::criterion ()
00826 {
00827                                    // dummy implementation. this function is
00828                                    // not needed for the present implementation
00829                                    // of gmres
00830   Assert (false, ExcInternalError());
00831   return 0;
00832 }
00833 
00834 
00835 //----------------------------------------------------------------------//
00836 
00837 template <class VECTOR>
00838 SolverFGMRES<VECTOR>::SolverFGMRES (SolverControl        &cn,
00839                                     VectorMemory<VECTOR> &mem,
00840                                     const AdditionalData &data)
00841                 :
00842                 Solver<VECTOR> (cn, mem),
00843                 additional_data(data)
00844 {}
00845 
00846 
00847 
00848 template <class VECTOR>
00849 SolverFGMRES<VECTOR>::SolverFGMRES (SolverControl        &cn,
00850                                     const AdditionalData &data)
00851                 :
00852                 Solver<VECTOR> (cn),
00853                 additional_data(data)
00854 {}
00855 
00856 
00857 
00858 template<class VECTOR>
00859 template<class MATRIX, class PRECONDITIONER>
00860 void
00861 SolverFGMRES<VECTOR>::solve (
00862   const MATRIX& A,
00863   VECTOR& x,
00864   const VECTOR& b,
00865   const PRECONDITIONER& precondition)
00866 {
00867   deallog.push("FGMRES");
00868 
00869   SolverControl::State iteration_state = SolverControl::iterate;
00870 
00871   const unsigned int basis_size = additional_data.max_basis_size;
00872 
00873                                    // Generate an object where basis
00874                                    // vectors are stored.
00875   typename internal::SolverGMRES::TmpVectors<VECTOR> v (basis_size, this->memory);
00876   typename internal::SolverGMRES::TmpVectors<VECTOR> z (basis_size, this->memory);
00877 
00878                                    // number of the present iteration; this
00879                                    // number is not reset to zero upon a
00880                                    // restart
00881   unsigned int accumulated_iterations = 0;
00882 
00883                                    // matrix used for the orthogonalization
00884                                    // process later
00885   H.reinit(basis_size+1, basis_size);
00886 
00887                                    // Vectors for projected system
00888   Vector<double> projected_rhs;
00889   Vector<double> y;
00890 
00891   // Iteration starts here
00892 
00893   VECTOR* aux = this->memory.alloc();
00894   aux->reinit(x);
00895   do
00896     {
00897       A.vmult(*aux, x);
00898       aux->sadd(-1., 1., b);
00899 
00900       double beta = aux->l2_norm();
00901       if (this->control().check(accumulated_iterations,beta)
00902           == SolverControl::success)
00903         break;
00904 
00905       H.reinit(basis_size+1, basis_size);
00906       double a = beta;
00907 
00908       for (unsigned int j=0;j<basis_size;++j)
00909         {
00910           v(j,x).equ(1./a, *aux);
00911 
00912           precondition.vmult(z(j,x), v[j]);
00913           A.vmult(*aux, z[j]);
00914 
00915                                            // Gram-Schmidt
00916           for (unsigned int i=0;i<=j;++i)
00917             {
00918               H(i,j) = *aux * v[i];
00919               aux->add(-H(i,j), v[i]);
00920             }
00921           H(j+1,j) = a = aux->l2_norm();
00922 
00923                                            // Compute projected solution
00924 
00925           if (j>0)
00926             {
00927               H1.reinit(j+1,j);
00928               projected_rhs.reinit(j+1);
00929               y.reinit(j);
00930               projected_rhs(0) = beta;
00931               H1.fill(H);
00932               Householder<double> house(H1);
00933               double res = house.least_squares(y, projected_rhs);
00934               iteration_state = this->control().check(++accumulated_iterations, res);
00935               if (iteration_state != SolverControl::iterate)
00936                 break;
00937             }
00938         }
00939                                        // Update solution vector
00940       for (unsigned int j=0;j<y.size();++j)
00941         x.add(y(j), z[j]);
00942 
00943     } while (iteration_state == SolverControl::iterate);
00944 
00945   this->memory.free(aux);
00946 
00947   deallog.pop();
00948                                    // in case of failure: throw
00949                                    // exception
00950   if (this->control().last_check() != SolverControl::success)
00951     throw SolverControl::NoConvergence (this->control().last_step(),
00952                                         this->control().last_value());
00953 }
00954 
00955 #endif // DOXYGEN
00956 
00957 DEAL_II_NAMESPACE_CLOSE
00958 
00959 #endif
00960 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:39 by doxygen 1.7.3