00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__solver_minres_h
00013 #define __deal2__solver_minres_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/lac/solver.h>
00018 #include <deal.II/lac/solver_control.h>
00019 #include <deal.II/base/logstream.h>
00020 #include <cmath>
00021 #include <deal.II/base/subscriptor.h>
00022
00023 DEAL_II_NAMESPACE_OPEN
00024
00027
00057 template <class VECTOR = Vector<double> >
00058 class SolverMinRes : public Solver<VECTOR>
00059 {
00060 public:
00067 struct AdditionalData
00068 {
00069 };
00070
00074 SolverMinRes (SolverControl &cn,
00075 VectorMemory<VECTOR> &mem,
00076 const AdditionalData &data=AdditionalData());
00077
00083 SolverMinRes (SolverControl &cn,
00084 const AdditionalData &data=AdditionalData());
00085
00089 virtual ~SolverMinRes ();
00090
00095 template<class MATRIX, class PRECONDITIONER>
00096 void
00097 solve (const MATRIX &A,
00098 VECTOR &x,
00099 const VECTOR &b,
00100 const PRECONDITIONER &precondition);
00101
00108 DeclException0 (ExcPreconditionerNotDefinite);
00110
00111 protected:
00116 virtual double criterion();
00126 virtual void print_vectors(const unsigned int step,
00127 const VECTOR& x,
00128 const VECTOR& r,
00129 const VECTOR& d) const;
00130
00137 VECTOR *Vu0, *Vu1, *Vu2;
00138 VECTOR *Vm0, *Vm1, *Vm2;
00139 VECTOR *Vv;
00140
00151 double res2;
00152 };
00153
00155
00156
00157 #ifndef DOXYGEN
00158
00159 template<class VECTOR>
00160 SolverMinRes<VECTOR>::SolverMinRes (SolverControl &cn,
00161 VectorMemory<VECTOR> &mem,
00162 const AdditionalData &)
00163 :
00164 Solver<VECTOR>(cn,mem)
00165 {}
00166
00167
00168
00169 template<class VECTOR>
00170 SolverMinRes<VECTOR>::SolverMinRes (SolverControl &cn,
00171 const AdditionalData &)
00172 :
00173 Solver<VECTOR>(cn)
00174 {}
00175
00176
00177 template<class VECTOR>
00178 SolverMinRes<VECTOR>::~SolverMinRes ()
00179 {}
00180
00181
00182
00183 template<class VECTOR>
00184 double
00185 SolverMinRes<VECTOR>::criterion()
00186 {
00187 return res2;
00188 }
00189
00190
00191 template<class VECTOR>
00192 void
00193 SolverMinRes<VECTOR>::print_vectors(const unsigned int,
00194 const VECTOR&,
00195 const VECTOR&,
00196 const VECTOR&) const
00197 {}
00198
00199
00200
00201 template<class VECTOR>
00202 template<class MATRIX, class PRECONDITIONER>
00203 void
00204 SolverMinRes<VECTOR>::solve (const MATRIX &A,
00205 VECTOR &x,
00206 const VECTOR &b,
00207 const PRECONDITIONER &precondition)
00208 {
00209 SolverControl::State conv=SolverControl::iterate;
00210
00211 deallog.push("minres");
00212
00213
00214 Vu0 = this->memory.alloc();
00215 Vu1 = this->memory.alloc();
00216 Vu2 = this->memory.alloc();
00217 Vv = this->memory.alloc();
00218 Vm0 = this->memory.alloc();
00219 Vm1 = this->memory.alloc();
00220 Vm2 = this->memory.alloc();
00221
00222 typedef VECTOR *vecptr;
00223 vecptr u[3] = {Vu0, Vu1, Vu2};
00224 vecptr m[3] = {Vm0, Vm1, Vm2};
00225 VECTOR &v = *Vv;
00226
00227
00228
00229 u[0]->reinit(b,true);
00230 u[1]->reinit(b,true);
00231 u[2]->reinit(b,true);
00232 m[0]->reinit(b,true);
00233 m[1]->reinit(b,true);
00234 m[2]->reinit(b,true);
00235 v.reinit(b,true);
00236
00237
00238 double delta[3] = { 0, 0, 0 };
00239 double f[2] = { 0, 0 };
00240 double e[2] = { 0, 0 };
00241
00242 double r_l2 = 0;
00243 double r0 = 0;
00244 double tau = 0;
00245 double c = 0;
00246 double gamma = 0;
00247 double s = 0;
00248 double d_ = 0;
00249 double d = 0;
00250
00251
00252 unsigned int j = 1;
00253
00254
00255
00256 A.vmult(*m[0],x);
00257 *u[1] = b;
00258 *u[1] -= *m[0];
00259
00260
00261
00262
00263
00264 precondition.vmult (v,*u[1]);
00265
00266 delta[1] = v * (*u[1]);
00267
00268 Assert (delta[1]>=0, ExcPreconditionerNotDefinite());
00269
00270 r0 = std::sqrt(delta[1]);
00271 r_l2 = r0;
00272
00273
00274 u[0]->reinit(b);
00275 delta[0] = 1.;
00276 m[0]->reinit(b);
00277 m[1]->reinit(b);
00278 m[2]->reinit(b);
00279
00280 conv = this->control().check(0,r_l2);
00281
00282 while (conv==SolverControl::iterate)
00283 {
00284 if (delta[1]!=0)
00285 v *= 1./std::sqrt(delta[1]);
00286 else
00287 v.reinit(b);
00288
00289 A.vmult(*u[2],v);
00290 u[2]->add (-std::sqrt(delta[1]/delta[0]), *u[0]);
00291
00292 gamma = *u[2] * v;
00293 u[2]->add (-gamma / std::sqrt(delta[1]), *u[1]);
00294 *m[0] = v;
00295
00296
00297
00298
00299 precondition.vmult(v,*u[2]);
00300
00301 delta[2] = v * (*u[2]);
00302
00303 Assert (delta[2]>=0, ExcPreconditionerNotDefinite());
00304
00305 if (j==1)
00306 {
00307 d_ = gamma;
00308 e[1] = std::sqrt(delta[2]);
00309 }
00310 if (j>1)
00311 {
00312 d_ = s * e[0] - c * gamma;
00313 e[0] = c * e[0] + s * gamma;
00314 f[1] = s * std::sqrt(delta[2]);
00315 e[1] = -c * std::sqrt(delta[2]);
00316 }
00317
00318 d = std::sqrt (d_*d_ + delta[2]);
00319
00320 if (j>1) tau *= s / c;
00321 c = d_ / d;
00322 tau *= c;
00323
00324 s = std::sqrt(delta[2]) / d;
00325
00326 if (j==1)
00327 tau = r0 * c;
00328
00329 m[0]->add (-e[0], *m[1]);
00330 if (j>1)
00331 m[0]->add (-f[0], *m[2]);
00332 *m[0] *= 1./d;
00333 x.add (tau, *m[0]);
00334 r_l2 *= std::fabs(s);
00335
00336 conv = this->control().check(j,r_l2);
00337
00338
00339 ++j;
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350 swap (*m[2], *m[1]);
00351 swap (*m[1], *m[0]);
00352
00353
00354
00355
00356 swap (*u[0], *u[1]);
00357 swap (*u[1], *u[2]);
00358
00359
00360
00361 f[0] = f[1];
00362 e[0] = e[1];
00363 delta[0] = delta[1];
00364 delta[1] = delta[2];
00365 }
00366
00367
00368 this->memory.free(Vu0);
00369 this->memory.free(Vu1);
00370 this->memory.free(Vu2);
00371 this->memory.free(Vv);
00372 this->memory.free(Vm0);
00373 this->memory.free(Vm1);
00374 this->memory.free(Vm2);
00375
00376 deallog.pop ();
00377
00378
00379
00380 if (this->control().last_check() != SolverControl::success)
00381 throw SolverControl::NoConvergence (this->control().last_step(),
00382 this->control().last_value());
00383
00384 }
00385
00386 #endif // DOXYGEN
00387
00388 DEAL_II_NAMESPACE_CLOSE
00389
00390 #endif
00391