include/deal.II/lac/solver_minres.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: solver_minres.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 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_minres_h
00013 #define __deal2__solver_minres_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/lac/solver.h>
00018 #include <deal.II/lac/solver_control.h>
00019 #include <deal.II/base/logstream.h>
00020 #include <cmath>
00021 #include <deal.II/base/subscriptor.h>
00022 
00023 DEAL_II_NAMESPACE_OPEN
00024 
00027 
00057 template <class VECTOR = Vector<double> >
00058 class SolverMinRes : public Solver<VECTOR>
00059 {
00060   public:
00067     struct AdditionalData
00068     {
00069     };
00070 
00074     SolverMinRes (SolverControl &cn,
00075                   VectorMemory<VECTOR> &mem,
00076                   const AdditionalData &data=AdditionalData());
00077 
00083     SolverMinRes (SolverControl        &cn,
00084                   const AdditionalData &data=AdditionalData());
00085 
00089     virtual ~SolverMinRes ();
00090 
00095     template<class MATRIX, class PRECONDITIONER>
00096     void
00097     solve (const MATRIX         &A,
00098            VECTOR               &x,
00099            const VECTOR         &b,
00100            const PRECONDITIONER &precondition);
00101 
00108     DeclException0 (ExcPreconditionerNotDefinite);
00110 
00111   protected:
00116     virtual double criterion();
00126     virtual void print_vectors(const unsigned int step,
00127                                const VECTOR& x,
00128                                const VECTOR& r,
00129                                const VECTOR& d) const;
00130 
00137     VECTOR *Vu0, *Vu1, *Vu2;
00138     VECTOR *Vm0, *Vm1, *Vm2;
00139     VECTOR *Vv;
00140 
00151     double res2;
00152 };
00153 
00155 /*------------------------- Implementation ----------------------------*/
00156 
00157 #ifndef DOXYGEN
00158 
00159 template<class VECTOR>
00160 SolverMinRes<VECTOR>::SolverMinRes (SolverControl &cn,
00161                                     VectorMemory<VECTOR> &mem,
00162                                     const AdditionalData &)
00163                 :
00164                 Solver<VECTOR>(cn,mem)
00165 {}
00166 
00167 
00168 
00169 template<class VECTOR>
00170 SolverMinRes<VECTOR>::SolverMinRes (SolverControl &cn,
00171                                     const AdditionalData &)
00172                 :
00173                 Solver<VECTOR>(cn)
00174 {}
00175 
00176 
00177 template<class VECTOR>
00178 SolverMinRes<VECTOR>::~SolverMinRes ()
00179 {}
00180 
00181 
00182 
00183 template<class VECTOR>
00184 double
00185 SolverMinRes<VECTOR>::criterion()
00186 {
00187   return res2;
00188 }
00189 
00190 
00191 template<class VECTOR>
00192 void
00193 SolverMinRes<VECTOR>::print_vectors(const unsigned int,
00194                                     const VECTOR&,
00195                                     const VECTOR&,
00196                                     const VECTOR&) const
00197 {}
00198 
00199 
00200 
00201 template<class VECTOR>
00202 template<class MATRIX, class PRECONDITIONER>
00203 void
00204 SolverMinRes<VECTOR>::solve (const MATRIX         &A,
00205                              VECTOR               &x,
00206                              const VECTOR         &b,
00207                              const PRECONDITIONER &precondition)
00208 {
00209   SolverControl::State conv=SolverControl::iterate;
00210 
00211   deallog.push("minres");
00212 
00213                                    // Memory allocation
00214   Vu0  = this->memory.alloc();
00215   Vu1  = this->memory.alloc();
00216   Vu2  = this->memory.alloc();
00217   Vv   = this->memory.alloc();
00218   Vm0  = this->memory.alloc();
00219   Vm1  = this->memory.alloc();
00220   Vm2  = this->memory.alloc();
00221                                    // define some aliases for simpler access
00222   typedef VECTOR *vecptr;
00223   vecptr u[3] = {Vu0, Vu1, Vu2};
00224   vecptr m[3] = {Vm0, Vm1, Vm2};
00225   VECTOR &v   = *Vv;
00226                                    // resize the vectors, but do not set
00227                                    // the values since they'd be overwritten
00228                                    // soon anyway.
00229   u[0]->reinit(b,true);
00230   u[1]->reinit(b,true);
00231   u[2]->reinit(b,true);
00232   m[0]->reinit(b,true);
00233   m[1]->reinit(b,true);
00234   m[2]->reinit(b,true);
00235   v.reinit(b,true);
00236 
00237                                    // some values needed
00238   double delta[3] = { 0, 0, 0 };
00239   double f[2] = { 0, 0 };
00240   double e[2] = { 0, 0 };
00241 
00242   double r_l2 = 0;
00243   double r0   = 0;
00244   double tau = 0;
00245   double c    = 0;
00246   double gamma = 0;
00247   double s = 0;
00248   double d_ = 0;
00249   double d = 0;
00250 
00251                                    // The iteration step.
00252   unsigned int j = 1;
00253 
00254 
00255                                    // Start of the solution process
00256   A.vmult(*m[0],x);
00257   *u[1] = b;
00258   *u[1] -= *m[0];
00259                                    // Precondition is applied.
00260                                    // The preconditioner has to be
00261                                    // positiv definite and symmetric
00262 
00263                                    // M v = u[1]
00264   precondition.vmult (v,*u[1]);
00265 
00266   delta[1] = v * (*u[1]);
00267                                    // Preconditioner positive
00268   Assert (delta[1]>=0, ExcPreconditionerNotDefinite());
00269 
00270   r0 = std::sqrt(delta[1]);
00271   r_l2 = r0;
00272 
00273 
00274   u[0]->reinit(b);
00275   delta[0] = 1.;
00276   m[0]->reinit(b);
00277   m[1]->reinit(b);
00278   m[2]->reinit(b);
00279 
00280   conv = this->control().check(0,r_l2);
00281 
00282   while (conv==SolverControl::iterate)
00283     {
00284       if (delta[1]!=0)
00285         v *= 1./std::sqrt(delta[1]);
00286       else
00287         v.reinit(b);
00288 
00289       A.vmult(*u[2],v);
00290       u[2]->add (-std::sqrt(delta[1]/delta[0]), *u[0]);
00291 
00292       gamma = *u[2] * v;
00293       u[2]->add (-gamma / std::sqrt(delta[1]), *u[1]);
00294       *m[0] = v;
00295 
00296                                        // precondition: solve M v = u[2]
00297                                        // Preconditioner has to be positiv
00298                                        // definite and symmetric.
00299       precondition.vmult(v,*u[2]);
00300 
00301       delta[2] = v * (*u[2]);
00302 
00303       Assert (delta[2]>=0, ExcPreconditionerNotDefinite());
00304 
00305       if (j==1)
00306         {
00307           d_ = gamma;
00308           e[1] = std::sqrt(delta[2]);
00309         }
00310       if (j>1)
00311         {
00312           d_ = s * e[0] - c * gamma;
00313           e[0] = c * e[0] + s * gamma;
00314           f[1] = s * std::sqrt(delta[2]);
00315           e[1] = -c * std::sqrt(delta[2]);
00316         }
00317 
00318       d = std::sqrt (d_*d_ + delta[2]);
00319 
00320       if (j>1) tau *= s / c;
00321       c = d_ / d;
00322       tau *= c;
00323 
00324       s = std::sqrt(delta[2]) / d;
00325 
00326       if (j==1)
00327         tau = r0 * c;
00328 
00329       m[0]->add (-e[0], *m[1]);
00330       if (j>1)
00331         m[0]->add (-f[0], *m[2]);
00332       *m[0] *= 1./d;
00333       x.add (tau, *m[0]);
00334       r_l2 *= std::fabs(s);
00335 
00336       conv = this->control().check(j,r_l2);
00337 
00338                                        // next iteration step
00339       ++j;
00340                                        // All vectors have to be shifted
00341                                        // one iteration step.
00342                                        // This should be changed one time.
00343                                        //
00344                                        // the previous code was like this:
00345                                        //   m[2] = m[1];
00346                                        //   m[1] = m[0];
00347                                        // but it can be made more efficient,
00348                                        // since the value of m[0] is no more
00349                                        // needed in the next iteration
00350       swap (*m[2], *m[1]);
00351       swap (*m[1], *m[0]);
00352 
00353                                        // likewise, but reverse direction:
00354                                        //   u[0] = u[1];
00355                                        //   u[1] = u[2];
00356       swap (*u[0], *u[1]);
00357       swap (*u[1], *u[2]);
00358 
00359                                        // these are scalars, so need
00360                                        // to bother
00361       f[0] = f[1];
00362       e[0] = e[1];
00363       delta[0] = delta[1];
00364       delta[1] = delta[2];
00365     }
00366 
00367                                    // Deallocation of Memory
00368   this->memory.free(Vu0);
00369   this->memory.free(Vu1);
00370   this->memory.free(Vu2);
00371   this->memory.free(Vv);
00372   this->memory.free(Vm0);
00373   this->memory.free(Vm1);
00374   this->memory.free(Vm2);
00375                                    // Output
00376   deallog.pop ();
00377 
00378                                    // in case of failure: throw
00379                                    // exception
00380   if (this->control().last_check() != SolverControl::success)
00381     throw SolverControl::NoConvergence (this->control().last_step(),
00382                                         this->control().last_value());
00383                                    // otherwise exit as normal
00384 }
00385 
00386 #endif // DOXYGEN
00387 
00388 DEAL_II_NAMESPACE_CLOSE
00389 
00390 #endif
00391 
 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