17 #ifndef __deal2__solver_gmres_h 18 #define __deal2__solver_gmres_h 22 #include <deal.II/base/config.h> 23 #include <deal.II/base/subscriptor.h> 24 #include <deal.II/base/logstream.h> 25 #include <deal.II/lac/householder.h> 26 #include <deal.II/lac/solver.h> 27 #include <deal.II/lac/solver_control.h> 28 #include <deal.II/lac/full_matrix.h> 29 #include <deal.II/lac/vector.h> 54 template <
class VECTOR>
73 VECTOR &
operator[] (
const unsigned int i)
const;
93 std::vector<VECTOR *>
data;
155 template <
class VECTOR = Vector<
double> >
171 const bool right_preconditioning =
false,
172 const bool use_default_residual =
true,
173 const bool force_re_orthogonalization =
false);
222 template<
class MATRIX,
class PRECONDITIONER>
227 const PRECONDITIONER &precondition);
231 <<
"The number of temporary vectors you gave (" 232 << arg1 <<
") is too small. It should be at least 10 for " 233 <<
"any results, and much more for reasonable ones.");
244 virtual double criterion();
265 const unsigned int dim,
266 const unsigned int accumulated_iterations,
269 bool &re_orthogonalize);
308 template <
class VECTOR = Vector<
double> >
325 max_basis_size(max_basis_size)
351 template<
class MATRIX,
class PRECONDITIONER>
356 const PRECONDITIONER &precondition);
382 template <
class VECTOR>
394 template <
class VECTOR>
398 for (
typename std::vector<VECTOR *>::iterator v =
data.begin();
399 v !=
data.end(); ++v)
405 template <
class VECTOR>
417 template <
class VECTOR>
436 template <
class VECTOR>
440 const bool right_preconditioning,
441 const bool use_default_residual,
442 const bool force_re_orthogonalization)
444 max_n_tmp_vectors(max_n_tmp_vectors),
445 right_preconditioning(right_preconditioning),
446 use_default_residual(use_default_residual),
447 force_re_orthogonalization(force_re_orthogonalization)
451 template <
class VECTOR>
457 additional_data(data)
462 template <
class VECTOR>
466 additional_data(data)
471 template <
class VECTOR>
480 for (
int i=0 ; i<col ; i++)
482 const double s = si(i);
483 const double c = ci(i);
484 const double dummy = h(i);
485 h(i) = c*dummy + s*h(i+1);
486 h(i+1) = -s*dummy + c*h(i+1);
489 const double r = 1./std::sqrt(h(col)*h(col) + h(col+1)*h(col+1));
490 si(col) = h(col+1) *r;
492 h(col) = ci(col)*h(col) + si(col)*h(col+1);
493 b(col+1)= -si(col)*b(col);
499 template <
class VECTOR>
503 const unsigned int dim,
504 const unsigned int accumulated_iterations,
507 bool &re_orthogonalize)
509 const unsigned int inner_iteration = dim - 1;
512 double norm_vv_start = 0;
513 if (re_orthogonalize ==
false && inner_iteration % 5 == 4)
514 norm_vv_start = vv.l2_norm();
517 for (
unsigned int i=0 ; i<dim ; ++i)
519 h(i) = vv * orthogonal_vectors[i];
520 vv.add(-h(i), orthogonal_vectors[i]);
530 if (re_orthogonalize ==
false && inner_iteration % 5 == 4)
532 const double norm_vv = vv.l2_norm();
533 if (norm_vv > 10. * norm_vv_start *
534 std::sqrt(std::numeric_limits<typename VECTOR::value_type>::epsilon()))
539 re_orthogonalize =
true;
540 deallog <<
"Re-orthogonalization enabled at step " 541 << accumulated_iterations << std::endl;
545 if (re_orthogonalize ==
true)
546 for (
unsigned int i=0 ; i<dim ; ++i)
548 double htmp = vv * orthogonal_vectors[i];
550 vv.add(-htmp, orthogonal_vectors[i]);
558 template<
class VECTOR>
559 template<
class MATRIX,
class PRECONDITIONER>
564 const PRECONDITIONER &precondition)
573 deallog.
push(
"GMRES");
574 const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors;
582 unsigned int accumulated_iterations = 0;
585 H.reinit(n_tmp_vectors, n_tmp_vectors-1);
589 gamma(n_tmp_vectors),
590 ci (n_tmp_vectors-1),
591 si (n_tmp_vectors-1),
595 unsigned int dim = 0;
601 const bool left_precondition = !additional_data.right_preconditioning;
606 const bool use_default_residual = additional_data.use_default_residual;
609 VECTOR &v = tmp_vectors(0, x);
610 VECTOR &p = tmp_vectors(n_tmp_vectors-1, x);
618 if (!use_default_residual)
620 r=this->memory.alloc();
621 x_=this->memory.alloc();
625 gamma_ = new ::Vector<double> (gamma.size());
628 bool re_orthogonalize = additional_data.force_re_orthogonalization;
637 h.
reinit (n_tmp_vectors-1);
639 if (left_precondition)
643 precondition.vmult(v,p);
651 double rho = v.l2_norm();
656 if (use_default_residual)
658 iteration_state = this->control().check (
659 accumulated_iterations, rho);
666 deallog <<
"default_res=" << rho << std::endl;
668 if (left_precondition)
674 precondition.vmult(*r,v);
676 double res = r->l2_norm();
677 iteration_state = this->control().check (
678 accumulated_iterations, res);
682 this->memory.free(r);
683 this->memory.free(x_);
697 for (
unsigned int inner_iteration=0;
698 ((inner_iteration < n_tmp_vectors-2)
703 ++accumulated_iterations;
705 VECTOR &vv = tmp_vectors(inner_iteration+1, x);
707 if (left_precondition)
709 A.
vmult(p, tmp_vectors[inner_iteration]);
710 precondition.vmult(vv,p);
714 precondition.vmult(p, tmp_vectors[inner_iteration]);
718 dim = inner_iteration+1;
720 const double s = modified_gram_schmidt(tmp_vectors, dim,
721 accumulated_iterations,
722 vv, h, re_orthogonalize);
723 h(inner_iteration+1) = s;
731 givens_rotation(h,gamma,ci,si,inner_iteration);
734 for (
unsigned int i=0; i<dim; ++i)
735 H(i,inner_iteration) = h(i);
738 rho = std::fabs(gamma(dim));
740 if (use_default_residual)
741 iteration_state = this->control().check (
742 accumulated_iterations, rho);
745 deallog <<
"default_res=" << rho << std::endl;
752 for (
unsigned int i=0; i<dim+1; ++i)
753 for (
unsigned int j=0; j<dim; ++j)
756 H1.backward(h_,*gamma_);
758 if (left_precondition)
759 for (
unsigned int i=0 ; i<dim; ++i)
760 x_->add(h_(i), tmp_vectors[i]);
764 for (
unsigned int i=0; i<dim; ++i)
765 p.add(h_(i), tmp_vectors[i]);
766 precondition.vmult(*r,p);
772 if (left_precondition)
774 const double res=r->l2_norm();
776 iteration_state = this->control().check (
777 accumulated_iterations, res);
781 precondition.vmult(*x_, *r);
782 const double preconditioned_res=x_->l2_norm();
784 iteration_state = this->control().check (
785 accumulated_iterations, preconditioned_res);
792 H1.reinit(dim+1,dim);
794 for (
unsigned int i=0; i<dim+1; ++i)
795 for (
unsigned int j=0; j<dim; ++j)
798 H1.backward(h,gamma);
800 if (left_precondition)
801 for (
unsigned int i=0 ; i<dim; ++i)
802 x.add(h(i), tmp_vectors[i]);
806 for (
unsigned int i=0; i<dim; ++i)
807 p.add(h(i), tmp_vectors[i]);
808 precondition.vmult(v,p);
816 if (!use_default_residual)
818 this->memory.free(r);
819 this->memory.free(x_);
828 this->control().last_value()));
834 template<
class VECTOR>
847 template <
class VECTOR>
853 additional_data(data)
858 template <
class VECTOR>
863 additional_data(data)
868 template<
class VECTOR>
869 template<
class MATRIX,
class PRECONDITIONER>
875 const PRECONDITIONER &precondition)
877 deallog.
push(
"FGMRES");
881 const unsigned int basis_size = additional_data.max_basis_size;
889 unsigned int accumulated_iterations = 0;
892 H.reinit(basis_size+1, basis_size);
900 VECTOR *aux = this->memory.alloc();
905 aux->sadd(-1., 1., b);
907 double beta = aux->l2_norm();
908 if (this->control().check(accumulated_iterations,beta)
912 H.reinit(basis_size+1, basis_size);
915 for (
unsigned int j=0; j<basis_size; ++j)
918 v(j,x).equ(1./a, *aux);
923 precondition.vmult(z(j,x), v[j]);
927 for (
unsigned int i=0; i<=j; ++i)
929 H(i,j) = *aux * v[i];
930 aux->add(-H(i,j), v[i]);
932 H(j+1,j) = a = aux->l2_norm();
939 projected_rhs.
reinit(j+1);
941 projected_rhs(0) = beta;
945 iteration_state = this->control().check(++accumulated_iterations, res);
951 for (
unsigned int j=0; j<y.
size(); ++j)
957 this->memory.free(aux);
963 this->control().last_value()));
968 DEAL_II_NAMESPACE_CLOSE
AdditionalData additional_data
bool right_preconditioning
void vmult(VECTOR &u, const VECTOR &v) const
virtual double criterion()
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
virtual VECTOR * alloc()=0
#define AssertThrow(cond, exc)
SolverFGMRES(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
void givens_rotation(Vector< double > &h, Vector< double > &b, Vector< double > &ci, Vector< double > &si, int col) const
bool is_finite(const double x)
bool use_default_residual
#define DeclException1(Exception1, type1, outsequence)
AdditionalData(const unsigned int max_basis_size=30, const bool=true)
virtual void reinit(const size_type N, const bool fast=false)
Stop iteration, goal reached.
#define Assert(cond, exc)
AdditionalData additional_data
SolverGMRES(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
bool force_re_orthogonalization
VECTOR & operator()(const unsigned int i, const VECTOR &temp)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
std::vector< VECTOR * > data
TmpVectors(const unsigned int max_size, VectorMemory< VECTOR > &vmem)
VECTOR & operator[](const unsigned int i) const
virtual void free(const VECTOR *const)=0
static double modified_gram_schmidt(const internal::SolverGMRES::TmpVectors< VECTOR > &orthogonal_vectors, const unsigned int dim, const unsigned int accumulated_iterations, VECTOR &vv, Vector< double > &h, bool &re_orthogonalize)
unsigned int max_basis_size
unsigned int max_n_tmp_vectors
void push(const std::string &text)
::ExceptionBase & ExcNotInitialized()
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
::ExceptionBase & ExcInternalError()
AdditionalData(const unsigned int max_n_tmp_vectors=30, const bool right_preconditioning=false, const bool use_default_residual=true, const bool force_re_orthogonalization=false)
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
VectorMemory< VECTOR > & mem