00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__solver_qmrs_h
00013 #define __deal2__solver_qmrs_h
00014
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/lac/solver.h>
00017 #include <deal.II/lac/solver_control.h>
00018 #include <deal.II/base/logstream.h>
00019 #include <cmath>
00020 #include <deal.II/base/subscriptor.h>
00021
00022 #include <cmath>
00023
00024 DEAL_II_NAMESPACE_OPEN
00025
00028
00061 template <class VECTOR = Vector<double> >
00062 class SolverQMRS : public Solver<VECTOR>
00063 {
00064 public:
00085 struct AdditionalData
00086 {
00094 AdditionalData(bool exact_residual = false,
00095 double breakdown=1.e-16) :
00096 exact_residual(exact_residual),
00097 breakdown(breakdown)
00098 {}
00099
00103 bool exact_residual;
00104
00108 double breakdown;
00109 };
00110
00114 SolverQMRS (SolverControl &cn,
00115 VectorMemory<VECTOR> &mem,
00116 const AdditionalData &data=AdditionalData());
00117
00123 SolverQMRS (SolverControl &cn,
00124 const AdditionalData &data=AdditionalData());
00125
00130 template<class MATRIX, class PRECONDITIONER>
00131 void
00132 solve (const MATRIX &A,
00133 VECTOR &x,
00134 const VECTOR &b,
00135 const PRECONDITIONER &precondition);
00136
00146 virtual void print_vectors(const unsigned int step,
00147 const VECTOR& x,
00148 const VECTOR& r,
00149 const VECTOR& d) const;
00150 protected:
00155 virtual double criterion();
00156
00163 VECTOR *Vv;
00164 VECTOR *Vp;
00165 VECTOR *Vq;
00166 VECTOR *Vt;
00167 VECTOR *Vd;
00171 VECTOR *Vx;
00175 const VECTOR *Vb;
00176
00187 double res2;
00191 AdditionalData additional_data;
00192 private:
00196 template<class MATRIX, class PRECONDITIONER>
00197 bool
00198 iterate(const MATRIX& A, const PRECONDITIONER& precondition);
00202 unsigned int step;
00203 };
00204
00206
00207
00208 #ifndef DOXYGEN
00209
00210 template<class VECTOR>
00211 SolverQMRS<VECTOR>::SolverQMRS(SolverControl &cn,
00212 VectorMemory<VECTOR> &mem,
00213 const AdditionalData &data)
00214 :
00215 Solver<VECTOR>(cn,mem),
00216 additional_data(data)
00217 {}
00218
00219
00220
00221 template<class VECTOR>
00222 SolverQMRS<VECTOR>::SolverQMRS(SolverControl &cn,
00223 const AdditionalData &data)
00224 :
00225 Solver<VECTOR>(cn),
00226 additional_data(data)
00227 {}
00228
00229
00230
00231 template<class VECTOR>
00232 double
00233 SolverQMRS<VECTOR>::criterion()
00234 {
00235 return std::sqrt(res2);
00236 }
00237
00238
00239
00240 template<class VECTOR>
00241 void
00242 SolverQMRS<VECTOR>::print_vectors(const unsigned int,
00243 const VECTOR&,
00244 const VECTOR&,
00245 const VECTOR&) const
00246 {}
00247
00248
00249
00250 template<class VECTOR>
00251 template<class MATRIX, class PRECONDITIONER>
00252 void
00253 SolverQMRS<VECTOR>::solve (const MATRIX &A,
00254 VECTOR &x,
00255 const VECTOR &b,
00256 const PRECONDITIONER &precondition)
00257 {
00258 deallog.push("QMRS");
00259
00260
00261 Vv = this->memory.alloc();
00262 Vp = this->memory.alloc();
00263 Vq = this->memory.alloc();
00264 Vt = this->memory.alloc();
00265 Vd = this->memory.alloc();
00266
00267 Vx = &x;
00268 Vb = &b;
00269
00270
00271
00272 Vv->reinit(x, true);
00273 Vp->reinit(x, true);
00274 Vq->reinit(x, true);
00275 Vt->reinit(x, true);
00276
00277 step = 0;
00278
00279 bool state;
00280
00281 do
00282 {
00283 if (step)
00284 deallog << "Restart step " << step << std::endl;
00285 state = iterate(A, precondition);
00286 }
00287 while (state);
00288
00289
00290 this->memory.free(Vv);
00291 this->memory.free(Vp);
00292 this->memory.free(Vq);
00293 this->memory.free(Vt);
00294 this->memory.free(Vd);
00295
00296
00297 deallog.pop();
00298
00299
00300
00301 if (this->control().last_check() != SolverControl::success)
00302 throw SolverControl::NoConvergence (this->control().last_step(),
00303 this->control().last_value());
00304
00305 }
00306
00307
00308
00309 template<class VECTOR>
00310 template<class MATRIX, class PRECONDITIONER>
00311 bool
00312 SolverQMRS<VECTOR>::iterate(const MATRIX &A,
00313 const PRECONDITIONER &precondition)
00314 {
00315
00316
00317
00318
00319
00320
00321 SolverControl::State state = SolverControl::iterate;
00322
00323
00324 VECTOR& v = *Vv;
00325 VECTOR& p = *Vp;
00326 VECTOR& q = *Vq;
00327 VECTOR& t = *Vt;
00328 VECTOR& d = *Vd;
00329 VECTOR& x = *Vx;
00330 const VECTOR& b = *Vb;
00331
00332 int it=0;
00333
00334 double tau, rho, theta=0, sigma, alpha, psi, theta_old, rho_old, beta;
00335 double res;
00336
00337 d.reinit(x);
00338
00339
00340 precondition.vmult(q,x);
00341
00342 A.vmult(v,q);
00343 v.sadd(-1.,1.,b);
00344 res = v.l2_norm();
00345
00346 if (this->control().check(step, res) == SolverControl::success)
00347 return false;
00348
00349 p = v;
00350
00351 precondition.vmult(q,p);
00352
00353 tau = v.norm_sqr();
00354 rho = q*v;
00355
00356 while (state == SolverControl::iterate)
00357 {
00358 step++; it++;
00359
00360 A.vmult(t,q);
00361
00362 sigma = q*t;
00363
00364
00365 if (std::fabs(sigma/rho) < additional_data.breakdown)
00366 return true;
00367
00368 alpha = rho/sigma;
00369
00370 v.add(-alpha,t);
00371
00372 theta_old = theta;
00373 theta = v*v/tau;
00374 psi = 1./(1.+theta);
00375 tau *= theta*psi;
00376
00377 d.sadd(psi*theta_old, psi*alpha, p);
00378 x.add(d);
00379
00380 print_vectors(step,x,v,d);
00381
00382 if (additional_data.exact_residual)
00383 {
00384 A.vmult(q,x);
00385 q.sadd(-1.,1.,b);
00386 res = q.l2_norm();
00387 }
00388 else
00389 res = std::sqrt((it+1)*tau);
00390 state = this->control().check(step,res);
00391 if ((state == SolverControl::success)
00392 || (state == SolverControl::failure))
00393 return false;
00394
00395 if (std::fabs(rho) < additional_data.breakdown)
00396 return true;
00397
00398 rho_old = rho;
00399 precondition.vmult(q,v);
00400 rho = q*v;
00401
00402 beta = rho/rho_old;
00403 p.sadd(beta,v);
00404 precondition.vmult(q,p);
00405 }
00406 return false;
00407 }
00408
00409 #endif // DOXYGEN
00410
00411 DEAL_II_NAMESPACE_CLOSE
00412
00413 #endif
00414