00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__eigen_h
00013 #define __deal2__eigen_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/lac/shifted_matrix.h>
00018 #include <deal.II/lac/solver.h>
00019 #include <deal.II/lac/solver_control.h>
00020 #include <deal.II/lac/solver_cg.h>
00021 #include <deal.II/lac/solver_gmres.h>
00022 #include <deal.II/lac/solver_minres.h>
00023 #include <deal.II/lac/vector_memory.h>
00024 #include <deal.II/lac/precondition.h>
00025
00026 #include <cmath>
00027
00028 DEAL_II_NAMESPACE_OPEN
00029
00030
00033
00049 template <class VECTOR = Vector<double> >
00050 class EigenPower : private Solver<VECTOR>
00051 {
00052 public:
00058 struct AdditionalData
00059 {
00066 double shift;
00070 AdditionalData (const double shift = 0.)
00071 :
00072 shift(shift)
00073 {}
00074
00075 };
00076
00080 EigenPower (SolverControl &cn,
00081 VectorMemory<VECTOR> &mem,
00082 const AdditionalData &data=AdditionalData());
00083
00087 virtual ~EigenPower ();
00088
00100 template <class MATRIX>
00101 void
00102 solve (double &value,
00103 const MATRIX &A,
00104 VECTOR &x);
00105
00106 protected:
00110 AdditionalData additional_data;
00111 };
00112
00138 template <class VECTOR = Vector<double> >
00139 class EigenInverse : private Solver<VECTOR>
00140 {
00141 public:
00147 struct AdditionalData
00148 {
00152 double relaxation;
00153
00158 unsigned int start_adaption;
00162 bool use_residual;
00166 AdditionalData (double relaxation = 1.,
00167 unsigned int start_adaption = 6,
00168 bool use_residual = true)
00169 :
00170 relaxation(relaxation),
00171 start_adaption(start_adaption),
00172 use_residual(use_residual)
00173 {}
00174
00175 };
00176
00180 EigenInverse (SolverControl &cn,
00181 VectorMemory<VECTOR> &mem,
00182 const AdditionalData &data=AdditionalData());
00183
00184
00188 virtual ~EigenInverse ();
00189
00203 template <class MATRIX>
00204 void
00205 solve (double &value,
00206 const MATRIX &A,
00207 VECTOR &x);
00208
00209 protected:
00213 AdditionalData additional_data;
00214 };
00215
00217
00218
00219
00220 template <class VECTOR>
00221 EigenPower<VECTOR>::EigenPower (SolverControl &cn,
00222 VectorMemory<VECTOR> &mem,
00223 const AdditionalData &data)
00224 :
00225 Solver<VECTOR>(cn, mem),
00226 additional_data(data)
00227 {}
00228
00229
00230
00231 template <class VECTOR>
00232 EigenPower<VECTOR>::~EigenPower ()
00233 {}
00234
00235
00236
00237 template <class VECTOR>
00238 template <class MATRIX>
00239 void
00240 EigenPower<VECTOR>::solve (double &value,
00241 const MATRIX &A,
00242 VECTOR &x)
00243 {
00244 SolverControl::State conv=SolverControl::iterate;
00245
00246 deallog.push("Power method");
00247
00248 VECTOR* Vy = this->memory.alloc (); VECTOR& y = *Vy; y.reinit (x);
00249 VECTOR* Vr = this->memory.alloc (); VECTOR& r = *Vr; r.reinit (x);
00250
00251 double length = x.l2_norm ();
00252 double old_length = 0.;
00253 x.scale(1./length);
00254
00255 A.vmult (y,x);
00256
00257
00258 for(int iter=0; conv==SolverControl::iterate; iter++)
00259 {
00260 y.add(additional_data.shift, x);
00261
00262
00263 old_length = length;
00264 length = y.l2_norm ();
00265
00266
00267
00268 double entry = 0.;
00269 unsigned int i = 0;
00270 double thresh = length/x.size();
00271 do
00272 {
00273 Assert (i<x.size(), ExcInternalError());
00274 entry = y (i++);
00275 }
00276 while (std::fabs(entry) < thresh);
00277
00278 --i;
00279
00280
00281 value = (entry * x (i) < 0.) ? -length : length;
00282 value -= additional_data.shift;
00283
00284
00285 x.equ (1/length, y);
00286
00287
00288 A.vmult (y,x);
00289
00290
00291
00292 conv = this->control().check (iter, std::fabs(1./length-1./old_length));
00293 }
00294
00295 this->memory.free(Vy);
00296 this->memory.free(Vr);
00297
00298 deallog.pop();
00299
00300
00301
00302 if (this->control().last_check() != SolverControl::success)
00303 throw SolverControl::NoConvergence (this->control().last_step(),
00304 this->control().last_value());
00305
00306 }
00307
00308
00309
00310 template <class VECTOR>
00311 EigenInverse<VECTOR>::EigenInverse (SolverControl &cn,
00312 VectorMemory<VECTOR> &mem,
00313 const AdditionalData &data)
00314 :
00315 Solver<VECTOR>(cn, mem),
00316 additional_data(data)
00317 {}
00318
00319
00320
00321 template <class VECTOR>
00322 EigenInverse<VECTOR>::~EigenInverse ()
00323 {}
00324
00325
00326
00327 template <class VECTOR>
00328 template <class MATRIX>
00329 void
00330 EigenInverse<VECTOR>::solve (double &value,
00331 const MATRIX &A,
00332 VECTOR &x)
00333 {
00334 deallog.push("Wielandt");
00335
00336 SolverControl::State conv=SolverControl::iterate;
00337
00338
00339 ShiftedMatrix <MATRIX> A_s(A, -value);
00340
00341
00342 ReductionControl inner_control (5000, 1.e-16, 1.e-5, false, false);
00343 PreconditionIdentity prec;
00344 SolverGMRES<VECTOR>
00345 solver(inner_control, this->memory);
00346
00347
00348 unsigned int goal = additional_data.start_adaption;
00349
00350
00351 VECTOR* Vy = this->memory.alloc (); VECTOR& y = *Vy; y.reinit (x);
00352 VECTOR* Vr = this->memory.alloc (); VECTOR& r = *Vr; r.reinit (x);
00353
00354 double length = x.l2_norm ();
00355 double old_value = value;
00356
00357 x.scale(1./length);
00358
00359
00360 for (unsigned int iter=0; conv==SolverControl::iterate; iter++)
00361 {
00362 solver.solve (A_s, y, x, prec);
00363
00364
00365 length = y.l2_norm ();
00366
00367
00368
00369 double entry = 0.;
00370 unsigned int i = 0;
00371 double thresh = length/x.size();
00372 do
00373 {
00374 Assert (i<x.size(), ExcInternalError());
00375 entry = y (i++);
00376 }
00377 while (std::fabs(entry) < thresh);
00378
00379 --i;
00380
00381
00382 value = (entry * x (i) < 0.) ? -length : length;
00383 value = 1./value;
00384 value -= A_s.shift ();
00385
00386 if (iter==goal)
00387 {
00388 const double new_shift = - additional_data.relaxation * value
00389 + (1.-additional_data.relaxation) * A_s.shift();
00390 A_s.shift(new_shift);
00391 ++goal;
00392 }
00393
00394
00395 x.equ (1./length, y);
00396
00397 if (additional_data.use_residual)
00398 {
00399 y.equ (value, x);
00400 A.vmult(r,x);
00401 r.sadd(-1., value, x);
00402 double res = r.l2_norm();
00403
00404 conv = this->control().check (iter, res);
00405 } else {
00406 conv = this->control().check (iter, std::fabs(1./value-1./old_value));
00407 }
00408 old_value = value;
00409 }
00410
00411 this->memory.free(Vy);
00412 this->memory.free(Vr);
00413
00414 deallog.pop();
00415
00416
00417
00418 if (this->control().last_check() != SolverControl::success)
00419 throw SolverControl::NoConvergence (this->control().last_step(),
00420 this->control().last_value());
00421
00422 }
00423
00424 DEAL_II_NAMESPACE_CLOSE
00425
00426 #endif
00427