From: tcclevenger Date: Mon, 30 Sep 2019 07:08:26 +0000 (-0400) Subject: add idr(s) solver X-Git-Tag: v9.2.0-rc1~812^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=09590c27f49ea549db2edca9d153deaddba3ade3;p=dealii.git add idr(s) solver --- diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h index f6e1ad966e..f707151be3 100644 --- a/include/deal.II/lac/solver_gmres.h +++ b/include/deal.II/lac/solver_gmres.h @@ -95,7 +95,7 @@ namespace internal private: /** - * Pool were vectors are obtained from. + * Pool where vectors are obtained from. */ VectorMemory &mem; diff --git a/include/deal.II/lac/solver_idr.h b/include/deal.II/lac/solver_idr.h new file mode 100644 index 0000000000..617a69392d --- /dev/null +++ b/include/deal.II/lac/solver_idr.h @@ -0,0 +1,478 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_solver_idr_h +#define dealii_solver_idr_h + + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +/*!@addtogroup Solvers */ +/*@{*/ + +namespace internal +{ + /** + * A namespace for a helper class to the IDR(s) solver. + */ + namespace SolverIDRImplementation + { + /** + * Class to hold temporary vectors whose size depends on + * the solver parameter s. + */ + template + class TmpVectors + { + public: + /** + * Constructor. Prepares an array of @p VectorType of length @p s_param. + */ + TmpVectors(const unsigned int s_param, VectorMemory &vmem); + + /** + * Destructor. Delete all allocated vectors. + */ + ~TmpVectors() = default; + + /** + * Get vector number @p i. If this vector was unused before, an error + * occurs. + */ + VectorType &operator[](const unsigned int i) const; + + /** + * Get vector number @p i. Allocate it if necessary. + * + * If a vector must be allocated, @p temp is used to reinit it to the + * proper dimensions. + */ + VectorType & + operator()(const unsigned int i, const VectorType &temp); + + private: + /** + * Pool where vectors are obtained from. + */ + VectorMemory &mem; + + /** + * Field for storing the vectors. + */ + std::vector::Pointer> data; + }; + } // namespace SolverIDRImplementation +} // namespace internal + +/** + * This class implements the IDR(s) method used for solving nonsymmetric, + * indefinite linear systems, developed in + * IDR(s): A Family of Simple and Fast Algorithms for Solving Large + * Nonsymmetric Systems of Linear Equations by Martin B. van Gijzen and Peter + * Sonneveld . The implementation here is the preconditioned version from + * Algorithm 913: An Elegant IDR(s) Variant that Efficiently Exploits + * Biorthogonality Properties + * by Martin B. van Gijzen and Peter Sonneveld. The local structure + * @p AdditionalData takes the value for the parameter s which can be any + * integer greater than or equal to 1. For s=1, this method has + * similar convergence to BiCGStab. + * + * @note Each iteration of IDR(s) requires s+1 preconditioning steps + * and matrix-vector products. In this implementation the residual is updated + * and convergence is checked after each of these inner steps inside the outer + * iteration. If the user enables the history data, the residual at each of + * these steps is stored and therefore there will be multiple values per + * iteration. + * + * @author Conrad Clevenger, 2019 + */ +template > +class SolverIDR : public SolverBase +{ +public: + /** + * Structure for storing additional data needed by the solver. + */ + struct AdditionalData + { + /** + * Constructor. By default, an IDR(2) method is used. + */ + explicit AdditionalData(const unsigned int s = 2) + : s(s) + {} + + const unsigned int s; + }; + + /** + * Constructor. + */ + SolverIDR(SolverControl & cn, + VectorMemory &mem, + const AdditionalData & data = AdditionalData()); + + /** + * Constructor. Use an object of type GrowingVectorMemory as a default to + * allocate memory. + */ + explicit SolverIDR(SolverControl & cn, + const AdditionalData &data = AdditionalData()); + + /** + * Virtual destructor. + */ + virtual ~SolverIDR() override = default; + + /** + * Solve the linear system Ax=b for x. + */ + template + void + solve(const MatrixType & A, + VectorType & x, + const VectorType & b, + const PreconditionerType &preconditioner); + +protected: + /** + * Interface for derived class. This function gets the current iteration + * vector, the residual and the update vector in each step. It can be used + * for graphical output of the convergence history. + */ + virtual void + print_vectors(const unsigned int step, + const VectorType & x, + const VectorType & r, + const VectorType & d) const; + +private: + /** + * Additional solver parameters. + */ + AdditionalData additional_data; +}; + +/*@}*/ +/*------------------------- Implementation ----------------------------*/ + +#ifndef DOXYGEN + + +namespace internal +{ + namespace SolverIDRImplementation + { + template + inline TmpVectors::TmpVectors(const unsigned int s_param, + VectorMemory &vmem) + : mem(vmem) + , data(s_param) + {} + + + + template + inline VectorType &TmpVectors:: + operator[](const unsigned int i) const + { + Assert(i < data.size(), ExcIndexRange(i, 0, data.size())); + + Assert(data[i] != nullptr, ExcNotInitialized()); + return *data[i]; + } + + + + template + inline VectorType & + TmpVectors::operator()(const unsigned int i, + const VectorType & temp) + { + AssertIndexRange(i, data.size()); + if (data[i] == nullptr) + { + data[i] = std::move(typename VectorMemory::Pointer(mem)); + data[i]->reinit(temp); + } + return *data[i]; + } + } // namespace SolverIDRImplementation +} // namespace internal + + + +template +SolverIDR::SolverIDR(SolverControl & cn, + VectorMemory &mem, + const AdditionalData & data) + : SolverBase(cn, mem) + , additional_data(data) +{} + + + +template +SolverIDR::SolverIDR(SolverControl &cn, const AdditionalData &data) + : SolverBase(cn) + , additional_data(data) +{} + + + +template +void +SolverIDR::print_vectors(const unsigned int, + const VectorType &, + const VectorType &, + const VectorType &) const +{} + + + +template +template +void +SolverIDR::solve(const MatrixType & A, + VectorType & x, + const VectorType & b, + const PreconditionerType &preconditioner) +{ + LogStream::Prefix prefix("IDR(s)"); + + SolverControl::State iteration_state = SolverControl::iterate; + unsigned int step = 0; + + const unsigned int s = additional_data.s; + + // Define temporary vectors which do not do not depend on s + typename VectorMemory::Pointer r_pointer(this->memory); + typename VectorMemory::Pointer v_pointer(this->memory); + typename VectorMemory::Pointer vhat_pointer(this->memory); + typename VectorMemory::Pointer uhat_pointer(this->memory); + typename VectorMemory::Pointer ghat_pointer(this->memory); + + VectorType &r = *r_pointer; + VectorType &v = *v_pointer; + VectorType &vhat = *vhat_pointer; + VectorType &uhat = *uhat_pointer; + VectorType &ghat = *ghat_pointer; + + r.reinit(x, true); + v.reinit(x, true); + vhat.reinit(x, true); + uhat.reinit(x, true); + ghat.reinit(x, true); + + // Initial residual + A.vmult(r, x); + r.sadd(-1.0, 1.0, b); + + // Check for convergent initial guess + double res = r.l2_norm(); + iteration_state = this->iteration_status(step, res, x); + if (iteration_state == SolverControl::success) + return; + + // Initialize sets of vectors/matrices whose size dependent on s + internal::SolverIDRImplementation::TmpVectors G(s, this->memory); + internal::SolverIDRImplementation::TmpVectors U(s, this->memory); + internal::SolverIDRImplementation::TmpVectors Q(s, this->memory); + FullMatrix M(s, s); + + // Random number generator for vector entries of + // Q (normal distribution, mean=0 sigma=1) + std::mt19937 rng; + std::normal_distribution<> normal_distribution(0.0, 1.0); + for (unsigned int i = 0; i < s; ++i) + { + VectorType &tmp_g = G(i, x); + VectorType &tmp_u = U(i, x); + tmp_g = 0; + tmp_u = 0; + + // Compute random set of s orthonormalized vectors Q + // Note: the first vector is chosen to be the initial + // residual to match BiCGStab (as is done in comparisons + // with BiCGStab in the papers listed in the documentation + // of this function) + VectorType &tmp_q = Q(i, x); + if (i != 0) + { + for (auto indx : tmp_q.locally_owned_elements()) + tmp_q(indx) = normal_distribution(rng); + tmp_q.compress(VectorOperation::insert); + } + else + tmp_q = r; + + for (unsigned int j = 0; j < i; ++j) + { + v = Q[j]; + v *= (v * tmp_q) / (tmp_q * tmp_q); + tmp_q.add(-1.0, v); + } + + if (i != 0) + tmp_q *= 1.0 / tmp_q.l2_norm(); + + M(i, i) = 1.; + } + + double omega = 1.; + + bool early_exit = false; + + // Outer iteration + while (iteration_state == SolverControl::iterate) + { + ++step; + + // Compute phi + Vector phi(s); + for (unsigned int i = 0; i < s; ++i) + phi(i) = Q[i] * r; + + // Inner iteration over s + for (unsigned int k = 0; k < s; ++k) + { + // Solve M(k:s)*gamma = phi(k:s) + Vector gamma(s - k); + { + Vector phik(s - k); + FullMatrix Mk(s - k, s - k); + std::vector indices; + unsigned int j = 0; + for (unsigned int i = k; i < s; ++i, ++j) + { + indices.push_back(i); + phik(j) = phi(i); + } + Mk.extract_submatrix_from(M, indices, indices); + + FullMatrix Mk_inv(s - k, s - k); + Mk_inv.invert(Mk); + Mk_inv.vmult(gamma, phik); + } + + { + v = r; + + unsigned int j = 0; + for (unsigned int i = k; i < s; ++i, ++j) + v.add(-1.0 * gamma(j), G[i]); + preconditioner.vmult(vhat, v); + + uhat = vhat; + uhat *= omega; + j = 0; + for (unsigned int i = k; i < s; ++i, ++j) + uhat.add(gamma(j), U[i]); + A.vmult(ghat, uhat); + } + + // Update G and U + // Orthogonalize ghat to Q0,..,Q_{k-1} + // and update uhat + for (unsigned int i = 0; i < k; ++i) + { + double alpha = (Q[i] * ghat) / M(i, i); + ghat.add(-alpha, G[i]); + uhat.add(-alpha, U[i]); + } + G[k] = ghat; + U[k] = uhat; + + // Update kth column of M + for (unsigned int i = k; i < s; ++i) + M(i, k) = Q[i] * G[k]; + + // Orthogonalize r to Q0,...,Qk, + // update x + { + double beta = phi(k) / M(k, k); + r.add(-1.0 * beta, G[k]); + x.add(beta, U[k]); + + print_vectors(step, x, r, U[k]); + + // Check for early convergence. If so, store + // information in early_exit so that outer iteration + // is broken before recomputing the residual + res = r.l2_norm(); + iteration_state = this->iteration_status(step, res, x); + if (iteration_state != SolverControl::iterate) + { + early_exit = true; + break; + } + + // Update phi + if (k + 1 < s) + { + for (unsigned int i = 0; i < k + 1; ++i) + phi(i) = 0.0; + for (unsigned int i = k + 1; i < s; ++i) + phi(i) -= beta * M(i, k); + } + } + } + if (early_exit == true) + break; + + // Update r and x + preconditioner.vmult(vhat, r); + A.vmult(v, vhat); + + omega = (v * r) / (v * v); + + r.add(-1.0 * omega, v); + x.add(omega, vhat); + + print_vectors(step, x, r, vhat); + + // Check for convergence + res = r.l2_norm(); + iteration_state = this->iteration_status(step, res, x); + if (iteration_state != SolverControl::iterate) + break; + } + + if (iteration_state != SolverControl::success) + AssertThrow(false, SolverControl::NoConvergence(step, res)); +} + + +#endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/lac/solver_02.cc b/tests/lac/solver_02.cc index 46fce2fc6d..a046cc8136 100644 --- a/tests/lac/solver_02.cc +++ b/tests/lac/solver_02.cc @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -66,4 +67,5 @@ main() test>>(); test>>(); test>>(); + test>>(); } diff --git a/tests/lac/solver_02.output b/tests/lac/solver_02.output index c8977da68f..d097b92c72 100644 --- a/tests/lac/solver_02.output +++ b/tests/lac/solver_02.output @@ -1,16 +1,21 @@ DEAL:GMRES::Check 0 1.000 DEAL:GMRES::Starting value 1.000 -DEAL:GMRES::Check 1 0 -DEAL:GMRES::Convergence step 1 value 0 -DEAL::0 0 1.000 +DEAL:GMRES::Check 1 0.000 +DEAL:GMRES::Convergence step 1 value 0.000 +DEAL::0.000 0.000 1.000 DEAL:cg::Check 0 1.000 DEAL:cg::Starting value 1.000 -DEAL:cg::Check 1 0 -DEAL:cg::Convergence step 1 value 0 -DEAL::0 0 1.000 +DEAL:cg::Check 1 0.000 +DEAL:cg::Convergence step 1 value 0.000 +DEAL::0.000 0.000 1.000 DEAL:FGMRES::Check 0 1.000 DEAL:FGMRES::Starting value 1.000 -DEAL:FGMRES::Check 1 0 -DEAL:FGMRES::Convergence step 1 value 0 -DEAL::0 0 1.000 +DEAL:FGMRES::Check 1 0.000 +DEAL:FGMRES::Convergence step 1 value 0.000 +DEAL::0.000 0.000 1.000 +DEAL:IDR(s)::Check 0 1.000 +DEAL:IDR(s)::Starting value 1.000 +DEAL:IDR(s)::Check 1 0.000 +DEAL:IDR(s)::Convergence step 1 value 0.000 +DEAL::0.000 0.000 1.000