00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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
00304 Vr = this->memory.alloc();
00305 Vp = this->memory.alloc();
00306 Vz = this->memory.alloc();
00307 VAp = this->memory.alloc();
00308
00309
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
00316
00317 std::vector<double> diagonal;
00318 std::vector<double> offdiagonal;
00319
00320 try {
00321
00322 VECTOR& g = *Vr;
00323 VECTOR& h = *Vp;
00324 VECTOR& d = *Vz;
00325 VECTOR& Ad = *VAp;
00326
00327
00328
00329 g.reinit(x, true);
00330 h.reinit(x, true);
00331 d.reinit(x, true);
00332 Ad.reinit(x, true);
00333
00334
00335 int it=0;
00336 double res,gh,alpha,beta;
00337
00338
00339
00340
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
00392
00393
00394
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
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
00450 cleanup();
00451
00452
00453 if (this->control().last_check() != SolverControl::success)
00454 throw SolverControl::NoConvergence (this->control().last_step(),
00455 this->control().last_value());
00456
00457 }
00458
00459 #endif // DOXYGEN
00460
00461 DEAL_II_NAMESPACE_CLOSE
00462
00463 #endif
00464