NGSolve  4.9
linalg/cg.hpp
00001 #ifndef FILE_CG_NEW
00002 #define FILE_CG_NEW
00003 
00004 
00005 /**************************************************************************/
00006 /* File:   cg.hh                                                          */
00007 /* Author: Joachim Schoeberl                                              */
00008 /* Date:   5. Jul. 96                                                     */
00009 /**************************************************************************/
00010 
00011 namespace ngla
00012 {
00013 
00014 
00018   class NGS_DLL_HEADER KrylovSpaceSolver : public BaseMatrix
00019   {
00020   protected:
00022     const BaseMatrix *a, *c;
00024     double prec;
00026     int maxsteps;
00028     int steps;
00030     int initialize;
00032     bool stop_absolute;
00034     int printrates;
00036     int absoluteRes;
00038     bool useseed;
00039 
00041     const BaseStatusHandler * sh;
00042 
00043   public:
00045     KrylovSpaceSolver ();
00047     KrylovSpaceSolver (const BaseMatrix & aa);
00049     KrylovSpaceSolver (const BaseMatrix & aa, const BaseMatrix & ac);
00051     void SetMatrix (const BaseMatrix & aa)
00052     { a = &aa; }
00054     void SetPrecond (const BaseMatrix & ac)
00055     { c = ∾ }
00057     void SetMaxSteps (int amaxsteps)
00058     { maxsteps = amaxsteps; }
00059 
00061     void SetPrecision (double aprec)
00062     { prec = aprec; stop_absolute = 0; }
00064     void SetAbsolutePrecision (double aprec)
00065     {  prec = aprec; stop_absolute = 1; }
00067     void SetRelativePrecision (double aprec)
00068     {  prec = aprec; stop_absolute = 0; }
00069 
00070 
00071     void SetPrintRates (int pr = 1)
00072     { printrates = pr; }
00074     void SetInitialize (int ai)
00075     { initialize = ai; }
00077     void SetStatusHandler (const BaseStatusHandler & stha)
00078     { sh = &stha; }
00079     
00080 
00081     void UseSeed(const bool useit = true)
00082     { useseed = useit; }
00083 
00085     int GetSteps () const
00086     { return steps; }
00088     virtual void Mult (const BaseVector & v, BaseVector & prod) const = 0;
00090     virtual BaseVector * CreateVector () const;
00091 
00092 
00093     virtual int VHeight() const {return a->VWidth();}
00094     virtual int VWidth() const {return a->VHeight();}
00095   };
00096   
00097 
00099   template <class IPTYPE>
00100   class NGS_DLL_HEADER CGSolver : public KrylovSpaceSolver
00101   {
00102   protected:
00104     void MultiMult (const BaseVector & f, BaseVector & u, const int dim) const;
00106     void MultiMultSeed (const BaseVector & f, BaseVector & u, const int dim) const;
00107   public:
00108     typedef typename SCAL_TRAIT<IPTYPE>::SCAL SCAL;
00110     CGSolver () 
00111       : KrylovSpaceSolver () { ; }
00113     CGSolver (const BaseMatrix & aa)
00114       : KrylovSpaceSolver (aa) { ; }
00115 
00117     CGSolver (const BaseMatrix & aa, const BaseMatrix & ac)
00118       : KrylovSpaceSolver (aa, ac) { ; }
00119 
00121     virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00122   };
00123 
00124 
00126   template <class IPTYPE>
00127   class BiCGStabSolver : public KrylovSpaceSolver
00128   {
00129   public:
00130     typedef typename SCAL_TRAIT<IPTYPE>::SCAL SCAL;
00132     BiCGStabSolver () 
00133       : KrylovSpaceSolver () { ; }
00135     BiCGStabSolver (const BaseMatrix & aa)
00136       : KrylovSpaceSolver (aa) { ; }
00137 
00139     BiCGStabSolver (const BaseMatrix & aa, const BaseMatrix & ac)
00140       : KrylovSpaceSolver (aa, ac) { ; }
00141 
00143     virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00144   };
00145   
00146 
00147 
00148 
00149 
00150   //   Simple iteration solver
00151   template <class IPTYPE>
00152   class SimpleIterationSolver : public KrylovSpaceSolver
00153   {
00154   public:
00155     typedef typename SCAL_TRAIT<IPTYPE>::SCAL SCAL;
00156   private:
00157     SCAL tau;
00158   public:
00160     SimpleIterationSolver ()
00161       : KrylovSpaceSolver() { tau = 1; }
00163     SimpleIterationSolver (const BaseMatrix & aa)
00164       : KrylovSpaceSolver (aa) { tau = 1; }
00166     SimpleIterationSolver (const BaseMatrix & aa, const BaseMatrix & ac)
00167       : KrylovSpaceSolver (aa, ac) { tau = 1; }
00169     virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00170 
00171     void SetTau (SCAL atau) { tau = atau; }
00172   };
00173 
00174 
00175 
00177   template <class IPTYPE>
00178   class NGS_DLL_HEADER GMRESSolver : public KrylovSpaceSolver
00179   {
00180   public:
00181     typedef typename SCAL_TRAIT<IPTYPE>::SCAL SCAL;
00183     GMRESSolver () 
00184       : KrylovSpaceSolver () { ; }
00186     GMRESSolver (const BaseMatrix & aa)
00187       : KrylovSpaceSolver (aa) { ; }
00188 
00190     GMRESSolver (const BaseMatrix & aa, const BaseMatrix & ac)
00191       : KrylovSpaceSolver (aa, ac) { ; }
00192 
00194     virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00195   };
00196   
00197 
00198 
00199 
00200 
00202   template <class IPTYPE>
00203   class QMRSolver : public KrylovSpaceSolver
00204   {
00205     int status;
00206     const BaseMatrix * c2;
00207   public:
00208     typedef typename SCAL_TRAIT<IPTYPE>::SCAL SCAL;
00210     QMRSolver () 
00211       : KrylovSpaceSolver (), c2(0) { ; }
00213     QMRSolver (const BaseMatrix & aa)
00214       : KrylovSpaceSolver (aa), c2(0) { ; }
00215 
00217     QMRSolver (const BaseMatrix & aa, const BaseMatrix & ac)
00218       : KrylovSpaceSolver (aa, ac), c2(0) { ; }
00219 
00221     virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00222   };
00223   
00224 
00225 
00226 
00227 
00228   /*
00229 
00230   // Conjugate residual solver for symmetric, indefinite matrices
00231 
00232 
00233   class CRSolver : public CGSolver
00234   {
00235   public:
00237   CRSolver ();
00239   CRSolver (const BaseMatrix & aa);
00241   CRSolver (const BaseMatrix & aa, const BaseMatrix & ac);
00243   ~CRSolver ();
00245   virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00246   };
00247 
00248 
00249 
00250   //  CG Solver for normal problem 
00251 
00252 
00253   class NormalCGSolver : public CGSolver
00254   {
00255   public:
00257   NormalCGSolver ();
00259   NormalCGSolver (const BaseMatrix & aa);
00261   NormalCGSolver (const BaseMatrix & aa, const BaseMatrix & ac);
00263   ~NormalCGSolver ();
00265   virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00266   };
00267 
00268 
00269 
00271   class MinResGradientSolver : public CGSolver
00272   {
00273   public:
00275   MinResGradientSolver ();
00277   MinResGradientSolver (const BaseMatrix & aa);
00279   MinResGradientSolver (const BaseMatrix & aa, const BaseMatrix & ac);
00281   ~MinResGradientSolver ();
00283   virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00284   };
00285 
00286 
00287 
00288   //  Bi-CG Stab solver for non-symmetric matrices
00289 
00290   class BiCGStabSolver : public CGSolver
00291   {
00292   public:
00294   BiCGStabSolver ();
00296   BiCGStabSolver (const BaseMatrix & aa);
00298   BiCGStabSolver (const BaseMatrix & aa, const BaseMatrix & ac);
00300   ~BiCGStabSolver ();
00302   virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00303   };
00304 
00305 
00306 
00307 
00308 
00309 
00310 
00311 
00312 
00313   //  QMR solver for non-symmetric matrices
00314 
00315   class QMRSolver : public CGSolver
00316   {
00317   int status;
00319   const BaseMatrix * c2;
00320   public:
00322   QMRSolver ();
00324   QMRSolver (const BaseMatrix & aa);
00326   QMRSolver (const BaseMatrix & aa, const BaseMatrix & ac);
00328   ~QMRSolver ();
00330   virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00331   };
00332 
00333 
00334 
00335 
00336 
00337 
00338   //  complex QMR solver for non-symmetric matrices
00339 
00340   class ComplexQMRSolver : public CGSolver
00341   {
00342   public:
00344   int status;
00346   const BaseMatrix * c2;
00347   public:
00349   ComplexQMRSolver ();
00351   ComplexQMRSolver (const BaseMatrix & aa);
00353   ComplexQMRSolver (const BaseMatrix & aa, const BaseMatrix & ac);
00355   ~ComplexQMRSolver ();
00357   virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00358   };
00359 
00360   */
00361 
00362 
00363 }
00364 
00365 
00366 #endif