00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__solver_gmres_h
00013 #define __deal2__solver_gmres_h
00014
00015
00016
00017 #include <deal.II/base/config.h>
00018 #include <deal.II/base/subscriptor.h>
00019 #include <deal.II/base/logstream.h>
00020 #include <deal.II/lac/householder.h>
00021 #include <deal.II/lac/solver.h>
00022 #include <deal.II/lac/solver_control.h>
00023 #include <deal.II/lac/full_matrix.h>
00024 #include <deal.II/lac/vector.h>
00025
00026 #include <vector>
00027 #include <cmath>
00028
00029 DEAL_II_NAMESPACE_OPEN
00030
00033
00034 namespace internal
00035 {
00040 namespace SolverGMRES
00041 {
00054 template <class VECTOR>
00055 class TmpVectors
00056 {
00057 public:
00063 TmpVectors(const unsigned int max_size,
00064 VectorMemory<VECTOR> &vmem);
00065
00069 ~TmpVectors();
00070
00077 VECTOR& operator[] (const unsigned int i) const;
00078
00089 VECTOR& operator() (const unsigned int i,
00090 const VECTOR &temp);
00091
00092 private:
00097 VectorMemory<VECTOR> &mem;
00098
00103 std::vector<VECTOR*> data;
00104
00111 unsigned int offset;
00112 };
00113 }
00114 }
00115
00167 template <class VECTOR = Vector<double> >
00168 class SolverGMRES : public Solver<VECTOR>
00169 {
00170 public:
00176 struct AdditionalData
00177 {
00187 AdditionalData (const unsigned int max_n_tmp_vectors = 30,
00188 const bool right_preconditioning = false,
00189 const bool use_default_residual = true);
00190
00200 unsigned int max_n_tmp_vectors;
00201
00214 bool right_preconditioning;
00215
00221 bool use_default_residual;
00222 };
00223
00227 SolverGMRES (SolverControl &cn,
00228 VectorMemory<VECTOR> &mem,
00229 const AdditionalData &data=AdditionalData());
00230
00236 SolverGMRES (SolverControl &cn,
00237 const AdditionalData &data=AdditionalData());
00238
00243 template<class MATRIX, class PRECONDITIONER>
00244 void
00245 solve (const MATRIX &A,
00246 VECTOR &x,
00247 const VECTOR &b,
00248 const PRECONDITIONER &precondition);
00249
00250 DeclException1 (ExcTooFewTmpVectors,
00251 int,
00252 << "The number of temporary vectors you gave ("
00253 << arg1 << ") is too small. It should be at least 10 for "
00254 << "any results, and much more for reasonable ones.");
00255
00256 protected:
00261 AdditionalData additional_data;
00262
00267 virtual double criterion();
00268
00275 void givens_rotation (Vector<double> &h, Vector<double> &b,
00276 Vector<double> &ci, Vector<double> &si,
00277 int col) const;
00281 FullMatrix<double> H;
00285 FullMatrix<double> H1;
00286
00287 private:
00291 SolverGMRES (const SolverGMRES<VECTOR>&);
00292 };
00293
00313 template <class VECTOR = Vector<double> >
00314 class SolverFGMRES : public Solver<VECTOR>
00315 {
00316 public:
00322 struct AdditionalData
00323 {
00334 AdditionalData(const unsigned int max_basis_size = 30,
00335 const bool = true)
00336 :
00337 max_basis_size(max_basis_size)
00338 {}
00339
00344 unsigned int max_basis_size;
00345 };
00346
00350 SolverFGMRES (SolverControl &cn,
00351 VectorMemory<VECTOR> &mem,
00352 const AdditionalData &data=AdditionalData());
00353
00359 SolverFGMRES (SolverControl &cn,
00360 const AdditionalData &data=AdditionalData());
00361
00366 template<class MATRIX, class PRECONDITIONER>
00367 void
00368 solve (const MATRIX &A,
00369 VECTOR &x,
00370 const VECTOR &b,
00371 const PRECONDITIONER &precondition);
00372
00373 private:
00377 AdditionalData additional_data;
00381 FullMatrix<double> H;
00385 FullMatrix<double> H1;
00386 };
00387
00389
00390
00391
00392 #ifndef DOXYGEN
00393 namespace internal
00394 {
00395 namespace SolverGMRES
00396 {
00397 template <class VECTOR>
00398 inline
00399 TmpVectors<VECTOR>::
00400 TmpVectors (const unsigned int max_size,
00401 VectorMemory<VECTOR> &vmem)
00402 :
00403 mem(vmem),
00404 data (max_size, 0),
00405 offset(0)
00406 {}
00407
00408
00409 template <class VECTOR>
00410 inline
00411 TmpVectors<VECTOR>::~TmpVectors ()
00412 {
00413 for (typename std::vector<VECTOR*>::iterator v = data.begin();
00414 v != data.end(); ++v)
00415 if (*v != 0)
00416 mem.free(*v);
00417 }
00418
00419
00420 template <class VECTOR>
00421 inline VECTOR&
00422 TmpVectors<VECTOR>::operator[] (const unsigned int i) const
00423 {
00424 Assert (i+offset<data.size(),
00425 ExcIndexRange(i, -offset, data.size()-offset));
00426
00427 Assert (data[i-offset] != 0, ExcNotInitialized());
00428 return *data[i-offset];
00429 }
00430
00431
00432 template <class VECTOR>
00433 inline VECTOR&
00434 TmpVectors<VECTOR>::operator() (const unsigned int i,
00435 const VECTOR &temp)
00436 {
00437 Assert (i+offset<data.size(),
00438 ExcIndexRange(i,-offset, data.size()-offset));
00439 if (data[i-offset] == 0)
00440 {
00441 data[i-offset] = mem.alloc();
00442 data[i-offset]->reinit(temp);
00443 }
00444 return *data[i-offset];
00445 }
00446 }
00447 }
00448
00449
00450
00451 template <class VECTOR>
00452 inline
00453 SolverGMRES<VECTOR>::AdditionalData::
00454 AdditionalData (const unsigned int max_n_tmp_vectors,
00455 const bool right_preconditioning,
00456 const bool use_default_residual)
00457 :
00458 max_n_tmp_vectors(max_n_tmp_vectors),
00459 right_preconditioning(right_preconditioning),
00460 use_default_residual(use_default_residual)
00461 {}
00462
00463
00464 template <class VECTOR>
00465 SolverGMRES<VECTOR>::SolverGMRES (SolverControl &cn,
00466 VectorMemory<VECTOR> &mem,
00467 const AdditionalData &data)
00468 :
00469 Solver<VECTOR> (cn,mem),
00470 additional_data(data)
00471 {}
00472
00473
00474
00475 template <class VECTOR>
00476 SolverGMRES<VECTOR>::SolverGMRES (SolverControl &cn,
00477 const AdditionalData &data) :
00478 Solver<VECTOR> (cn),
00479 additional_data(data)
00480 {}
00481
00482
00483
00484 template <class VECTOR>
00485 inline
00486 void
00487 SolverGMRES<VECTOR>::givens_rotation (Vector<double> &h,
00488 Vector<double> &b,
00489 Vector<double> &ci,
00490 Vector<double> &si,
00491 int col) const
00492 {
00493 for (int i=0 ; i<col ; i++)
00494 {
00495 const double s = si(i);
00496 const double c = ci(i);
00497 const double dummy = h(i);
00498 h(i) = c*dummy + s*h(i+1);
00499 h(i+1) = -s*dummy + c*h(i+1);
00500 };
00501
00502 const double r = 1./std::sqrt(h(col)*h(col) + h(col+1)*h(col+1));
00503 si(col) = h(col+1) *r;
00504 ci(col) = h(col) *r;
00505 h(col) = ci(col)*h(col) + si(col)*h(col+1);
00506 b(col+1)= -si(col)*b(col);
00507 b(col) *= ci(col);
00508 }
00509
00510
00511
00512 template<class VECTOR>
00513 template<class MATRIX, class PRECONDITIONER>
00514 void
00515 SolverGMRES<VECTOR>::solve (const MATRIX &A,
00516 VECTOR &x,
00517 const VECTOR &b,
00518 const PRECONDITIONER &precondition)
00519 {
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530 deallog.push("GMRES");
00531 const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors;
00532
00533
00534
00535 internal::SolverGMRES::TmpVectors<VECTOR> tmp_vectors (n_tmp_vectors, this->memory);
00536
00537
00538
00539
00540 unsigned int accumulated_iterations = 0;
00541
00542
00543
00544 H.reinit(n_tmp_vectors, n_tmp_vectors-1);
00545
00546
00547
00548 ::Vector<double>
00549 gamma(n_tmp_vectors),
00550 ci (n_tmp_vectors-1),
00551 si (n_tmp_vectors-1),
00552 h (n_tmp_vectors-1);
00553
00554
00555 unsigned int dim = 0;
00556
00557 SolverControl::State iteration_state = SolverControl::iterate;
00558
00559
00560
00561
00562
00563 const bool left_precondition = !additional_data.right_preconditioning;
00564
00565
00566
00567
00568
00569
00570 const bool use_default_residual = additional_data.use_default_residual;
00571
00572
00573 VECTOR &v = tmp_vectors(0, x);
00574 VECTOR &p = tmp_vectors(n_tmp_vectors-1, x);
00575
00576
00577
00578
00579 VECTOR *r=0;
00580 VECTOR *x_=0;
00581 ::Vector<double> *gamma_=0;
00582 if (!use_default_residual)
00583 {
00584 r=this->memory.alloc();
00585 x_=this->memory.alloc();
00586 r->reinit(x);
00587 x_->reinit(x);
00588
00589 gamma_ = new ::Vector<double> (gamma.size());
00590 }
00591
00593
00594
00595
00596
00597
00598 do
00599 {
00600
00601
00602 h.reinit (n_tmp_vectors-1);
00603
00604 if (left_precondition)
00605 {
00606 A.vmult(p,x);
00607 p.sadd(-1.,1.,b);
00608 precondition.vmult(v,p);
00609 }
00610 else
00611 {
00612 A.vmult(v,x);
00613 v.sadd(-1.,1.,b);
00614 };
00615
00616 double rho = v.l2_norm();
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626 if (use_default_residual)
00627 {
00628 iteration_state = this->control().check (
00629 accumulated_iterations, rho);
00630
00631 if (iteration_state != SolverControl::iterate)
00632 break;
00633 }
00634 else
00635 {
00636 deallog << "default_res=" << rho << std::endl;
00637
00638 if (left_precondition)
00639 {
00640 A.vmult(*r,x);
00641 r->sadd(-1.,1.,b);
00642 }
00643 else
00644 precondition.vmult(*r,v);
00645
00646 double res = r->l2_norm();
00647 iteration_state = this->control().check (
00648 accumulated_iterations, res);
00649
00650 if (iteration_state != SolverControl::iterate)
00651 {
00652 this->memory.free(r);
00653 this->memory.free(x_);
00654
00655 delete gamma_;
00656 break;
00657 }
00658 }
00659
00660 gamma(0) = rho;
00661
00662 v *= 1./rho;
00663
00664
00665
00666
00667
00668
00669
00670
00671 for (unsigned int inner_iteration=0;
00672 ((inner_iteration < n_tmp_vectors-2)
00673 &&
00674 (iteration_state==SolverControl::iterate));
00675 ++inner_iteration)
00676 {
00677 ++accumulated_iterations;
00678
00679 VECTOR& vv = tmp_vectors(inner_iteration+1, x);
00680
00681 if (left_precondition)
00682 {
00683 A.vmult(p, tmp_vectors[inner_iteration]);
00684 precondition.vmult(vv,p);
00685 } else {
00686 precondition.vmult(p, tmp_vectors[inner_iteration]);
00687 A.vmult(vv,p);
00688 };
00689
00690 dim = inner_iteration+1;
00691
00692
00693 for (unsigned int i=0 ; i<dim ; ++i)
00694 {
00695 h(i) = vv * tmp_vectors[i];
00696 vv.add(-h(i), tmp_vectors[i]);
00697 };
00698
00699
00700 for (unsigned int i=0 ; i<dim ; ++i)
00701 {
00702 double htmp = vv * tmp_vectors[i];
00703 h(i) += htmp;
00704 vv.add(-htmp, tmp_vectors[i]);
00705 }
00706
00707 const double s = vv.l2_norm();
00708 h(inner_iteration+1) = s;
00709
00710
00711 vv *= 1./s;
00712
00713
00714
00715 givens_rotation(h,gamma,ci,si,inner_iteration);
00716
00717
00718 for (unsigned int i=0; i<dim; ++i)
00719 H(i,inner_iteration) = h(i);
00720
00721
00722 rho = std::fabs(gamma(dim));
00723
00724 if (use_default_residual)
00725 iteration_state = this->control().check (
00726 accumulated_iterations, rho);
00727 else
00728 {
00729 deallog << "default_res=" << rho << std::endl;
00730
00731 ::Vector<double> h_(dim);
00732 *x_=x;
00733 *gamma_=gamma;
00734 H1.reinit(dim+1,dim);
00735
00736 for (unsigned int i=0; i<dim+1; ++i)
00737 for (unsigned int j=0; j<dim; ++j)
00738 H1(i,j) = H(i,j);
00739
00740 H1.backward(h_,*gamma_);
00741
00742 if (left_precondition)
00743 for (unsigned int i=0 ; i<dim; ++i)
00744 x_->add(h_(i), tmp_vectors[i]);
00745 else
00746 {
00747 p = 0.;
00748 for (unsigned int i=0; i<dim; ++i)
00749 p.add(h_(i), tmp_vectors[i]);
00750 precondition.vmult(*r,p);
00751 x_->add(1.,*r);
00752 };
00753 A.vmult(*r,*x_);
00754 r->sadd(-1.,1.,b);
00755
00756
00757
00758 if (left_precondition)
00759 {
00760 const double res=r->l2_norm();
00761
00762 iteration_state = this->control().check (
00763 accumulated_iterations, res);
00764 }
00765 else
00766 {
00767 precondition.vmult(*x_, *r);
00768 const double preconditioned_res=x_->l2_norm();
00769
00770 iteration_state = this->control().check (
00771 accumulated_iterations, preconditioned_res);
00772 }
00773 }
00774 };
00775
00776
00777
00778 h.reinit(dim);
00779 H1.reinit(dim+1,dim);
00780
00781 for (unsigned int i=0; i<dim+1; ++i)
00782 for (unsigned int j=0; j<dim; ++j)
00783 H1(i,j) = H(i,j);
00784
00785 H1.backward(h,gamma);
00786
00787 if (left_precondition)
00788 for (unsigned int i=0 ; i<dim; ++i)
00789 x.add(h(i), tmp_vectors[i]);
00790 else
00791 {
00792 p = 0.;
00793 for (unsigned int i=0; i<dim; ++i)
00794 p.add(h(i), tmp_vectors[i]);
00795 precondition.vmult(v,p);
00796 x.add(1.,v);
00797 };
00798
00799
00800
00801 }
00802 while (iteration_state == SolverControl::iterate);
00803
00804 if (!use_default_residual)
00805 {
00806 this->memory.free(r);
00807 this->memory.free(x_);
00808
00809 delete gamma_;
00810 }
00811
00812 deallog.pop();
00813
00814
00815 if (this->control().last_check() != SolverControl::success)
00816 throw SolverControl::NoConvergence (this->control().last_step(),
00817 this->control().last_value());
00818
00819 }
00820
00821
00822
00823 template<class VECTOR>
00824 double
00825 SolverGMRES<VECTOR>::criterion ()
00826 {
00827
00828
00829
00830 Assert (false, ExcInternalError());
00831 return 0;
00832 }
00833
00834
00835
00836
00837 template <class VECTOR>
00838 SolverFGMRES<VECTOR>::SolverFGMRES (SolverControl &cn,
00839 VectorMemory<VECTOR> &mem,
00840 const AdditionalData &data)
00841 :
00842 Solver<VECTOR> (cn, mem),
00843 additional_data(data)
00844 {}
00845
00846
00847
00848 template <class VECTOR>
00849 SolverFGMRES<VECTOR>::SolverFGMRES (SolverControl &cn,
00850 const AdditionalData &data)
00851 :
00852 Solver<VECTOR> (cn),
00853 additional_data(data)
00854 {}
00855
00856
00857
00858 template<class VECTOR>
00859 template<class MATRIX, class PRECONDITIONER>
00860 void
00861 SolverFGMRES<VECTOR>::solve (
00862 const MATRIX& A,
00863 VECTOR& x,
00864 const VECTOR& b,
00865 const PRECONDITIONER& precondition)
00866 {
00867 deallog.push("FGMRES");
00868
00869 SolverControl::State iteration_state = SolverControl::iterate;
00870
00871 const unsigned int basis_size = additional_data.max_basis_size;
00872
00873
00874
00875 typename internal::SolverGMRES::TmpVectors<VECTOR> v (basis_size, this->memory);
00876 typename internal::SolverGMRES::TmpVectors<VECTOR> z (basis_size, this->memory);
00877
00878
00879
00880
00881 unsigned int accumulated_iterations = 0;
00882
00883
00884
00885 H.reinit(basis_size+1, basis_size);
00886
00887
00888 Vector<double> projected_rhs;
00889 Vector<double> y;
00890
00891
00892
00893 VECTOR* aux = this->memory.alloc();
00894 aux->reinit(x);
00895 do
00896 {
00897 A.vmult(*aux, x);
00898 aux->sadd(-1., 1., b);
00899
00900 double beta = aux->l2_norm();
00901 if (this->control().check(accumulated_iterations,beta)
00902 == SolverControl::success)
00903 break;
00904
00905 H.reinit(basis_size+1, basis_size);
00906 double a = beta;
00907
00908 for (unsigned int j=0;j<basis_size;++j)
00909 {
00910 v(j,x).equ(1./a, *aux);
00911
00912 precondition.vmult(z(j,x), v[j]);
00913 A.vmult(*aux, z[j]);
00914
00915
00916 for (unsigned int i=0;i<=j;++i)
00917 {
00918 H(i,j) = *aux * v[i];
00919 aux->add(-H(i,j), v[i]);
00920 }
00921 H(j+1,j) = a = aux->l2_norm();
00922
00923
00924
00925 if (j>0)
00926 {
00927 H1.reinit(j+1,j);
00928 projected_rhs.reinit(j+1);
00929 y.reinit(j);
00930 projected_rhs(0) = beta;
00931 H1.fill(H);
00932 Householder<double> house(H1);
00933 double res = house.least_squares(y, projected_rhs);
00934 iteration_state = this->control().check(++accumulated_iterations, res);
00935 if (iteration_state != SolverControl::iterate)
00936 break;
00937 }
00938 }
00939
00940 for (unsigned int j=0;j<y.size();++j)
00941 x.add(y(j), z[j]);
00942
00943 } while (iteration_state == SolverControl::iterate);
00944
00945 this->memory.free(aux);
00946
00947 deallog.pop();
00948
00949
00950 if (this->control().last_check() != SolverControl::success)
00951 throw SolverControl::NoConvergence (this->control().last_step(),
00952 this->control().last_value());
00953 }
00954
00955 #endif // DOXYGEN
00956
00957 DEAL_II_NAMESPACE_CLOSE
00958
00959 #endif
00960