include/deal.II/multigrid/mg_coarse.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: mg_coarse.h 25460 2012-04-27 09:19:06Z bangerth @f$
00003 //
00004 //    Copyright (C) 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 //---------------------------------------------------------------------------
00012 #ifndef __deal2__mg_coarse_h
00013 #define __deal2__mg_coarse_h
00014 
00015 
00016 #include <deal.II/lac/full_matrix.h>
00017 #include <deal.II/lac/matrix_lib.h>
00018 #include <deal.II/lac/householder.h>
00019 #include <deal.II/multigrid/mg_base.h>
00020 
00021 DEAL_II_NAMESPACE_OPEN
00022 
00025 
00037 template<class SOLVER, class VECTOR = Vector<double> >
00038 class MGCoarseGridLACIteration :  public MGCoarseGridBase<VECTOR>
00039 {
00040   public:
00044     MGCoarseGridLACIteration ();
00045 
00052     template<class MATRIX, class PRECOND>
00053     MGCoarseGridLACIteration (SOLVER        &,
00054                               const MATRIX  &,
00055                               const PRECOND &);
00056 
00060     ~MGCoarseGridLACIteration ();
00061 
00065     template<class MATRIX, class PRECOND>
00066     void initialize (SOLVER        &,
00067                      const MATRIX  &,
00068                      const PRECOND &);
00069 
00073     void clear ();
00074 
00082     void operator() (const unsigned int   level,
00083                      VECTOR       &dst,
00084                      const VECTOR &src) const;
00085 
00092     template <class MATRIX>
00093     void set_matrix (const MATRIX &);
00094 
00095   private:
00099     SmartPointer<SOLVER,MGCoarseGridLACIteration<SOLVER,VECTOR> > solver;
00100 
00104     PointerMatrixBase<VECTOR>* matrix;
00105 
00109     PointerMatrixBase<VECTOR>* precondition;
00110 };
00111 
00112 
00113 
00123 template<typename number = double, class VECTOR = Vector<number> >
00124 class MGCoarseGridHouseholder : public MGCoarseGridBase<VECTOR>
00125 {
00126   public:
00131     MGCoarseGridHouseholder (const FullMatrix<number>* A = 0);
00132 
00136     void initialize (const FullMatrix<number>& A);
00137 
00138     void operator() (const unsigned int   level,
00139                      VECTOR       &dst,
00140                      const VECTOR &src) const;
00141 
00142   private:
00146     Householder<number> householder;
00147 };
00148 
00157 template<typename number = double, class VECTOR = Vector<number> >
00158 class MGCoarseGridSVD : public MGCoarseGridBase<VECTOR>
00159 {
00160   public:
00165     MGCoarseGridSVD ();
00166 
00172     void initialize (const FullMatrix<number>& A, const double threshold = 0);
00173 
00174     void operator() (const unsigned int   level,
00175                      VECTOR       &dst,
00176                      const VECTOR &src) const;
00177 
00181     void log () const;
00182 
00183   private:
00184 
00188     LAPACKFullMatrix<number> matrix;
00189 };
00190 
00193 #ifndef DOXYGEN
00194 /* ------------------ Functions for MGCoarseGridLACIteration ------------ */
00195 
00196 
00197 template<class SOLVER, class VECTOR>
00198 MGCoarseGridLACIteration<SOLVER, VECTOR>
00199 ::MGCoarseGridLACIteration()
00200                 :
00201                 solver(0, typeid(*this).name()),
00202                 matrix(0),
00203                 precondition(0)
00204 {}
00205 
00206 
00207 template<class SOLVER, class VECTOR>
00208 template<class MATRIX, class PRECOND>
00209 MGCoarseGridLACIteration<SOLVER, VECTOR>
00210 ::MGCoarseGridLACIteration(SOLVER& s,
00211                            const MATRIX  &m,
00212                            const PRECOND &p)
00213                 :
00214                 solver(&s, typeid(*this).name())
00215 {
00216   matrix = new PointerMatrix<MATRIX, VECTOR>(&m);
00217   precondition = new PointerMatrix<PRECOND, VECTOR>(&p);
00218 }
00219 
00220 
00221 template<class SOLVER, class VECTOR>
00222 MGCoarseGridLACIteration<SOLVER, VECTOR>
00223 ::~MGCoarseGridLACIteration()
00224 {
00225   clear();
00226 }
00227 
00228 
00229 template<class SOLVER, class VECTOR>
00230 template<class MATRIX, class PRECOND>
00231 void
00232 MGCoarseGridLACIteration<SOLVER, VECTOR>
00233 ::initialize(SOLVER& s,
00234              const MATRIX  &m,
00235              const PRECOND &p)
00236 {
00237   solver = &s;
00238   if (matrix)
00239     delete matrix;
00240   matrix = new PointerMatrix<MATRIX, VECTOR>(&m);
00241   if (precondition)
00242     delete precondition;
00243   precondition = new PointerMatrix<PRECOND, VECTOR>(&p);
00244 }
00245 
00246 
00247 template<class SOLVER, class VECTOR>
00248 void
00249 MGCoarseGridLACIteration<SOLVER, VECTOR>
00250 ::clear()
00251 {
00252   solver = 0;
00253   if (matrix)
00254     delete matrix;
00255   matrix = 0;
00256   if (precondition)
00257     delete precondition;
00258   precondition = 0;
00259 }
00260 
00261 
00262 template<class SOLVER, class VECTOR>
00263 void
00264 MGCoarseGridLACIteration<SOLVER, VECTOR>
00265 ::operator() (const unsigned int    /* level */,
00266               VECTOR       &dst,
00267               const VECTOR &src) const
00268 {
00269   Assert(solver!=0, ExcNotInitialized());
00270   Assert(matrix!=0, ExcNotInitialized());
00271   Assert(precondition!=0, ExcNotInitialized());
00272   solver->solve(*matrix, dst, src, *precondition);
00273 }
00274 
00275 
00276 template<class SOLVER, class VECTOR>
00277 template<class MATRIX>
00278 void
00279 MGCoarseGridLACIteration<SOLVER, VECTOR>
00280 ::set_matrix(const MATRIX &m)
00281 {
00282   if (matrix)
00283     delete matrix;
00284   matrix = new PointerMatrix<MATRIX, VECTOR>(&m);
00285 }
00286 
00287 //---------------------------------------------------------------------------
00288 
00289 template<typename number, class VECTOR>
00290 MGCoarseGridHouseholder<number, VECTOR>::MGCoarseGridHouseholder(
00291   const FullMatrix<number>* A)
00292 {
00293   if (A != 0) householder.initialize(*A);
00294 }
00295 
00296 
00297 
00298 template<typename number, class VECTOR>
00299 void
00300 MGCoarseGridHouseholder<number, VECTOR>::initialize(
00301   const FullMatrix<number>& A)
00302 {
00303   householder.initialize(A);
00304 }
00305 
00306 
00307 
00308 template<typename number, class VECTOR>
00309 void
00310 MGCoarseGridHouseholder<number, VECTOR>::operator() (
00311   const unsigned int,
00312   VECTOR       &dst,
00313   const VECTOR &src) const
00314 {
00315   householder.least_squares(dst, src);
00316 }
00317 
00318 //---------------------------------------------------------------------------
00319 
00320 template<typename number, class VECTOR>
00321 inline
00322 MGCoarseGridSVD<number, VECTOR>::MGCoarseGridSVD()
00323 {}
00324 
00325 
00326 
00327 template<typename number, class VECTOR>
00328 void
00329 MGCoarseGridSVD<number, VECTOR>::initialize(
00330   const FullMatrix<number>& A,
00331   double threshold)
00332 {
00333   matrix.reinit(A.n_rows(), A.n_cols());
00334   matrix = A;
00335   matrix.compute_inverse_svd(threshold);
00336 }
00337 
00338 
00339 template<typename number, class VECTOR>
00340 void
00341 MGCoarseGridSVD<number, VECTOR>::operator() (
00342   const unsigned int,
00343   VECTOR       &dst,
00344   const VECTOR &src) const
00345 {
00346   matrix.vmult(dst, src);
00347 }
00348 
00349 
00350 template<typename number, class VECTOR>
00351 void
00352 MGCoarseGridSVD<number, VECTOR>::log() const
00353 {
00354   const unsigned int n = std::min(matrix.n_rows(), matrix.n_cols());
00355 
00356   for (unsigned int i=0;i<n;++i)
00357     deallog << ' ' << matrix.singular_value(i);
00358   deallog << std::endl;
00359 }
00360 
00361 
00362 #endif // DOXYGEN
00363 
00364 DEAL_II_NAMESPACE_CLOSE
00365 
00366 #endif
00367 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

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