00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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 ,
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