include/deal.II/lac/solver_cg.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: solver_cg.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 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__solver_cg_h
00013 #define __deal2__solver_cg_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/lac/tridiagonal_matrix.h>
00018 #include <deal.II/lac/solver.h>
00019 #include <deal.II/lac/solver_control.h>
00020 #include <deal.II/base/exceptions.h>
00021 #include <deal.II/base/logstream.h>
00022 #include <deal.II/base/subscriptor.h>
00023 #include <cmath>
00024 
00025 DEAL_II_NAMESPACE_OPEN
00026 
00029 
00070 template <class VECTOR = Vector<double> >
00071 class SolverCG : public Solver<VECTOR>
00072 {
00073   public:
00078     struct AdditionalData
00079     {
00085         bool log_coefficients;
00086 
00094         bool compute_condition_number;
00095 
00103         bool compute_all_condition_numbers;
00104 
00111         bool compute_eigenvalues;
00112 
00118         AdditionalData (const bool log_coefficients = false,
00119                         const bool compute_condition_number = false,
00120                         const bool compute_all_condition_numbers = false,
00121                         const bool compute_eigenvalues = false);
00122     };
00123 
00127     SolverCG (SolverControl        &cn,
00128               VectorMemory<VECTOR> &mem,
00129               const AdditionalData &data = AdditionalData());
00130 
00136     SolverCG (SolverControl        &cn,
00137               const AdditionalData &data=AdditionalData());
00138 
00142     virtual ~SolverCG ();
00143 
00148     template <class MATRIX, class PRECONDITIONER>
00149     void
00150     solve (const MATRIX         &A,
00151            VECTOR               &x,
00152            const VECTOR         &b,
00153            const PRECONDITIONER &precondition);
00154 
00155   protected:
00162     virtual double criterion();
00163 
00173     virtual void print_vectors(const unsigned int step,
00174                                const VECTOR& x,
00175                                const VECTOR& r,
00176                                const VECTOR& d) const;
00177 
00184     VECTOR *Vr;
00185     VECTOR *Vp;
00186     VECTOR *Vz;
00187     VECTOR *VAp;
00188 
00199     double res2;
00200 
00204     AdditionalData additional_data;
00205 
00206   private:
00207     void cleanup();
00208 };
00209 
00212 /*------------------------- Implementation ----------------------------*/
00213 
00214 #ifndef DOXYGEN
00215 
00216 template <class VECTOR>
00217 inline
00218 SolverCG<VECTOR>::AdditionalData::
00219 AdditionalData (const bool log_coefficients,
00220                 const bool compute_condition_number,
00221                 const bool compute_all_condition_numbers,
00222                 const bool compute_eigenvalues)
00223                 :
00224                 log_coefficients (log_coefficients),
00225                 compute_condition_number(compute_condition_number),
00226                 compute_all_condition_numbers(compute_all_condition_numbers),
00227                 compute_eigenvalues(compute_eigenvalues)
00228 {}
00229 
00230 
00231 
00232 template <class VECTOR>
00233 SolverCG<VECTOR>::SolverCG (SolverControl        &cn,
00234                             VectorMemory<VECTOR> &mem,
00235                             const AdditionalData &data)
00236                 :
00237                 Solver<VECTOR>(cn,mem),
00238                 additional_data(data)
00239 {}
00240 
00241 
00242 
00243 template <class VECTOR>
00244 SolverCG<VECTOR>::SolverCG (SolverControl        &cn,
00245                             const AdditionalData &data)
00246                 :
00247                 Solver<VECTOR>(cn),
00248                 additional_data(data)
00249 {}
00250 
00251 
00252 
00253 template <class VECTOR>
00254 SolverCG<VECTOR>::~SolverCG ()
00255 {}
00256 
00257 
00258 
00259 template <class VECTOR>
00260 double
00261 SolverCG<VECTOR>::criterion()
00262 {
00263   return std::sqrt(res2);
00264 }
00265 
00266 
00267 
00268 template <class VECTOR>
00269 void
00270 SolverCG<VECTOR>::cleanup()
00271 {
00272   this->memory.free(Vr);
00273   this->memory.free(Vp);
00274   this->memory.free(Vz);
00275   this->memory.free(VAp);
00276   deallog.pop();
00277 }
00278 
00279 
00280 
00281 template <class VECTOR>
00282 void
00283 SolverCG<VECTOR>::print_vectors(const unsigned int,
00284                                 const VECTOR&,
00285                                 const VECTOR&,
00286                                 const VECTOR&) const
00287 {}
00288 
00289 
00290 
00291 template <class VECTOR>
00292 template <class MATRIX, class PRECONDITIONER>
00293 void
00294 SolverCG<VECTOR>::solve (const MATRIX         &A,
00295                          VECTOR               &x,
00296                          const VECTOR         &b,
00297                          const PRECONDITIONER &precondition)
00298 {
00299   SolverControl::State conv=SolverControl::iterate;
00300 
00301   deallog.push("cg");
00302 
00303                                    // Memory allocation
00304   Vr  = this->memory.alloc();
00305   Vp  = this->memory.alloc();
00306   Vz  = this->memory.alloc();
00307   VAp = this->memory.alloc();
00308                                    // Should we build the matrix for
00309                                    // eigenvalue computations?
00310   bool do_eigenvalues = additional_data.compute_condition_number
00311                         | additional_data.compute_all_condition_numbers
00312                         | additional_data.compute_eigenvalues;
00313   double eigen_beta_alpha = 0;
00314 
00315                                    // vectors used for eigenvalue
00316                                    // computations
00317   std::vector<double> diagonal;
00318   std::vector<double> offdiagonal;
00319 
00320   try {
00321                                      // define some aliases for simpler access
00322     VECTOR& g  = *Vr;
00323     VECTOR& h  = *Vp;
00324     VECTOR& d  = *Vz;
00325     VECTOR& Ad = *VAp;
00326                                      // resize the vectors, but do not set
00327                                      // the values since they'd be overwritten
00328                                      // soon anyway.
00329     g.reinit(x, true);
00330     h.reinit(x, true);
00331     d.reinit(x, true);
00332     Ad.reinit(x, true);
00333                                      // Implementation taken from the DEAL
00334                                      // library
00335     int  it=0;
00336     double res,gh,alpha,beta;
00337 
00338                                      // compute residual. if vector is
00339                                      // zero, then short-circuit the
00340                                      // full computation
00341     if (!x.all_zero())
00342       {
00343         A.vmult(g,x);
00344         g.add(-1.,b);
00345       }
00346     else
00347       g.equ(-1.,b);
00348     res = g.l2_norm();
00349 
00350     conv = this->control().check(0,res);
00351     if (conv)
00352       {
00353         cleanup();
00354         return;
00355       }
00356 
00357     precondition.vmult(h,g);
00358 
00359     d.equ(-1.,h);
00360 
00361     gh = g*h;
00362 
00363     while (conv == SolverControl::iterate)
00364       {
00365         it++;
00366         A.vmult(Ad,d);
00367 
00368         alpha = d*Ad;
00369         Assert(alpha != 0., ExcDivideByZero());
00370         alpha = gh/alpha;
00371 
00372         g.add(alpha,Ad);
00373         x.add(alpha,d );
00374         res = g.l2_norm();
00375 
00376         print_vectors(it, x, g, d);
00377 
00378         conv = this->control().check(it,res);
00379         if (conv != SolverControl::iterate)
00380           break;
00381 
00382         precondition.vmult(h,g);
00383 
00384         beta = gh;
00385         Assert(beta != 0., ExcDivideByZero());
00386         gh   = g*h;
00387         beta = gh/beta;
00388 
00389         if (additional_data.log_coefficients)
00390           deallog << "alpha-beta:" << alpha << '\t' << beta << std::endl;
00391                                          // set up the vectors
00392                                          // containing the diagonal
00393                                          // and the off diagonal of
00394                                          // the projected matrix.
00395         if (do_eigenvalues)
00396           {
00397             diagonal.push_back(1./alpha + eigen_beta_alpha);
00398             eigen_beta_alpha = beta/alpha;
00399             offdiagonal.push_back(std::sqrt(beta)/alpha);
00400           }
00401 
00402         if (additional_data.compute_all_condition_numbers && (diagonal.size()>1))
00403           {
00404             TridiagonalMatrix<double> T(diagonal.size(), true);
00405             for (unsigned int i=0;i<diagonal.size();++i)
00406               {
00407                 T(i,i) = diagonal[i];
00408                 if (i< diagonal.size()-1)
00409                   T(i,i+1) = offdiagonal[i];
00410               }
00411             T.compute_eigenvalues();
00412             deallog << "Condition number estimate: " <<
00413               T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
00414           }
00415 
00416         d.sadd(beta,-1.,h);
00417       }
00418   }
00419   catch (...)
00420     {
00421       cleanup();
00422       throw;
00423     }
00424 
00425                                    // Write eigenvalues or condition number
00426   if (do_eigenvalues)
00427     {
00428       TridiagonalMatrix<double> T(diagonal.size(), true);
00429       for (unsigned int i=0;i<diagonal.size();++i)
00430         {
00431           T(i,i) = diagonal[i];
00432           if (i< diagonal.size()-1)
00433             T(i,i+1) = offdiagonal[i];
00434         }
00435       T.compute_eigenvalues();
00436       if (additional_data.compute_condition_number
00437           && ! additional_data.compute_all_condition_numbers
00438           && (diagonal.size() > 1))
00439         deallog << "Condition number estimate: " <<
00440           T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
00441       if (additional_data.compute_eigenvalues)
00442         {
00443           for (unsigned int i=0;i<T.n();++i)
00444             deallog << ' ' << T.eigenvalue(i);
00445           deallog << std::endl;
00446         }
00447     }
00448 
00449                                    // Deallocate Memory
00450   cleanup();
00451                                    // in case of failure: throw
00452                                    // exception
00453   if (this->control().last_check() != SolverControl::success)
00454     throw SolverControl::NoConvergence (this->control().last_step(),
00455                                         this->control().last_value());
00456                                    // otherwise exit as normal
00457 }
00458 
00459 #endif // DOXYGEN
00460 
00461 DEAL_II_NAMESPACE_CLOSE
00462 
00463 #endif
00464 
 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