IDRS.h

Go to the documentation of this file.

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2020 Chris Schoutrop <c.e.m.schoutrop@tue.nl>
// Copyright (C) 2020 Jens Wehner <j.wehner@esciencecenter.nl>
// Copyright (C) 2020 Jan van Dijk <j.v.dijk@tue.nl>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.


#ifndef EIGEN_IDRS_H
#define EIGEN_IDRS_H

namespace Eigen
{

        namespace internal
        {
                template<typename Vector, typename RealScalar>
                typename Vector::Scalar omega(const Vector& t, const Vector& s, RealScalar angle)
                {
                        using numext::abs;
                        typedef typename Vector::Scalar Scalar;
                        const RealScalar ns = s.norm();
                        const RealScalar nt = t.norm();
                        const Scalar ts = t.dot(s);
                        const RealScalar rho = abs(ts / (nt * ns));

                        if (rho < angle) {
                                if (ts == Scalar(0)) {
                                        return Scalar(0);
                                }
                                // Original relation for om is given by
                                // om = om * angle / rho;
                                // To alleviate potential (near) division by zero this can be rewritten as
                                // om = angle * (ns / nt) * (ts / abs(ts)) = angle * (ns / nt) * sgn(ts)
                                return angle * (ns / nt) * (ts / abs(ts));
                        }
                        return ts / (nt * nt);
                }

                template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
                bool idrs(const MatrixType& A, const Rhs& b, Dest& x, const Preconditioner& precond,
                        Index& iter,
                        typename Dest::RealScalar& relres, Index S, bool smoothing, typename Dest::RealScalar angle, bool replacement)
                {
                        typedef typename Dest::RealScalar RealScalar;
                        typedef typename Dest::Scalar Scalar;
                        typedef Matrix<Scalar, Dynamic, 1> VectorType;
                        typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> DenseMatrixType;
                        const Index N = b.size();
                        S = S < x.rows() ? S : x.rows();
                        const RealScalar tol = relres;
                        const Index maxit = iter;

                        Index replacements = 0;
                        bool trueres = false;

                        FullPivLU<DenseMatrixType> lu_solver;

                        DenseMatrixType P;
                        {
                                HouseholderQR<DenseMatrixType> qr(DenseMatrixType::Random(N, S));
                            P = (qr.householderQ() * DenseMatrixType::Identity(N, S));
                        }

                        const RealScalar normb = b.norm();

                        if (internal::isApprox(normb, RealScalar(0)))
                        {
                                //Solution is the zero vector
                                x.setZero();
                                iter = 0;
                                relres = 0;
                                return true;
                        }
                         // from http://homepage.tudelft.nl/1w5b5/IDRS/manual.pdf
                         // A peak in the residual is considered dangerously high if‖ri‖/‖b‖> C(tol/epsilon).
                         // With epsilon the
             // relative machine precision. The factor tol/epsilon corresponds to the size of a
             // finite precision number that is so large that the absolute round-off error in
             // this number, when propagated through the process, makes it impossible to
             // achieve the required accuracy.The factor C accounts for the accumulation of
             // round-off errors. This parameter has beenset to 10−3.
                         // mp is epsilon/C
                         // 10^3 * eps is very conservative, so normally no residual replacements will take place. 
                         // It only happens if things go very wrong. Too many restarts may ruin the convergence.
                        const RealScalar mp = RealScalar(1e3) * NumTraits<Scalar>::epsilon();



                        //Compute initial residual
                        const RealScalar tolb = tol * normb; //Relative tolerance
                        VectorType r = b - A * x;

                        VectorType x_s, r_s;

                        if (smoothing)
                        {
                                x_s = x;
                                r_s = r;
                        }

                        RealScalar normr = r.norm();

                        if (normr <= tolb)
                        {
                                //Initial guess is a good enough solution
                                iter = 0;
                                relres = normr / normb;
                                return true;
                        }

                        DenseMatrixType G = DenseMatrixType::Zero(N, S);
                        DenseMatrixType U = DenseMatrixType::Zero(N, S);
                        DenseMatrixType M = DenseMatrixType::Identity(S, S);
                        VectorType t(N), v(N);
                        Scalar om = 1.;

                        //Main iteration loop, guild G-spaces:
                        iter = 0;

                        while (normr > tolb && iter < maxit)
                        {
                                //New right hand size for small system:
                                VectorType f = (r.adjoint() * P).adjoint();

                                for (Index k = 0; k < S; ++k)
                                {
                                        //Solve small system and make v orthogonal to P:
                                        //c = M(k:s,k:s)\f(k:s);
                                        lu_solver.compute(M.block(k , k , S -k, S - k ));
                                        VectorType c = lu_solver.solve(f.segment(k , S - k ));
                                        //v = r - G(:,k:s)*c;
                                        v = r - G.rightCols(S - k ) * c;
                                        //Preconditioning
                                        v = precond.solve(v);

                                        //Compute new U(:,k) and G(:,k), G(:,k) is in space G_j
                                        U.col(k) = U.rightCols(S - k ) * c + om * v;
                                        G.col(k) = A * U.col(k );

                                        //Bi-Orthogonalise the new basis vectors:
                                        for (Index i = 0; i < k-1 ; ++i)
                                        {
                                                //alpha =  ( P(:,i)'*G(:,k) )/M(i,i);
                                                Scalar alpha = P.col(i ).dot(G.col(k )) / M(i, i );
                                                G.col(k ) = G.col(k ) - alpha * G.col(i );
                                                U.col(k ) = U.col(k ) - alpha * U.col(i );
                                        }

                                        //New column of M = P'*G  (first k-1 entries are zero)
                                        //M(k:s,k) = (G(:,k)'*P(:,k:s))';
                                        M.block(k , k , S - k , 1) = (G.col(k ).adjoint() * P.rightCols(S - k )).adjoint();

                                        if (internal::isApprox(M(k,k), Scalar(0)))
                                        {
                                                return false;
                                        }

                                        //Make r orthogonal to q_i, i = 0..k-1
                                        Scalar beta = f(k ) / M(k , k );
                                        r = r - beta * G.col(k );
                                        x = x + beta * U.col(k );
                                        normr = r.norm();

                                        if (replacement && normr > tolb / mp)
                                        {
                                                trueres = true;
                                        }

                                        //Smoothing:
                                        if (smoothing)
                                        {
                                                t = r_s - r;
                                                //gamma is a Scalar, but the conversion is not allowed
                                                Scalar gamma = t.dot(r_s) / t.norm();
                                                r_s = r_s - gamma * t;
                                                x_s = x_s - gamma * (x_s - x);
                                                normr = r_s.norm();
                                        }

                                        if (normr < tolb || iter == maxit)
                                        {
                                                break;
                                        }

                                        //New f = P'*r (first k  components are zero)
                                        if (k < S-1)
                                        {
                                                f.segment(k + 1, S - (k + 1) ) = f.segment(k + 1 , S - (k + 1)) - beta * M.block(k + 1 , k , S - (k + 1), 1);
                                        }
                                }//end for

                                if (normr < tolb || iter == maxit)
                                {
                                        break;
                                }

                                //Now we have sufficient vectors in G_j to compute residual in G_j+1
                                //Note: r is already perpendicular to P so v = r
                                //Preconditioning
                                v = r;
                                v = precond.solve(v);

                                //Matrix-vector multiplication:
                                t = A * v;

                                //Computation of a new omega
                                om = internal::omega(t, r, angle);

                                if (om == RealScalar(0.0))
                                {
                                        return false;
                                }

                                r = r - om * t;
                                x = x + om * v;
                                normr = r.norm();

                                if (replacement && normr > tolb / mp)
                                {
                                        trueres = true;
                                }

                                //Residual replacement?
                                if (trueres && normr < normb)
                                {
                                        r = b - A * x;
                                        trueres = false;
                                        replacements++;
                                }

                                //Smoothing:
                                if (smoothing)
                                {
                                        t = r_s - r;
                                        Scalar gamma = t.dot(r_s) /t.norm();
                                        r_s = r_s - gamma * t;
                                        x_s = x_s - gamma * (x_s - x);
                                        normr = r_s.norm();
                                }

                                iter++;

                        }//end while

                        if (smoothing)
                        {
                                x = x_s;
                        }
                        relres=normr/normb;
                        return true;
                }

        }  // namespace internal

        template <typename _MatrixType, typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
        class IDRS;

        namespace internal
        {

                template <typename _MatrixType, typename _Preconditioner>
                struct traits<Eigen::IDRS<_MatrixType, _Preconditioner> >
                {
                        typedef _MatrixType MatrixType;
                        typedef _Preconditioner Preconditioner;
                };

        }  // namespace internal


        template <typename _MatrixType, typename _Preconditioner>
        class IDRS : public IterativeSolverBase<IDRS<_MatrixType, _Preconditioner> >
        {

                public:
                        typedef _MatrixType MatrixType;
                        typedef typename MatrixType::Scalar Scalar;
                        typedef typename MatrixType::RealScalar RealScalar;
                        typedef _Preconditioner Preconditioner;

                private:
                        typedef IterativeSolverBase<IDRS> Base;
                        using Base::m_error;
                        using Base::m_info;
                        using Base::m_isInitialized;
                        using Base::m_iterations;
                        using Base::matrix;
                        Index m_S;
                        bool m_smoothing;
                        RealScalar m_angle;
                        bool m_residual;

                public:
                        IDRS(): m_S(4), m_smoothing(false), m_angle(RealScalar(0.7)), m_residual(false) {}

                        template <typename MatrixDerived>
                        explicit IDRS(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_S(4), m_smoothing(false),
                                                                                                                           m_angle(RealScalar(0.7)), m_residual(false) {}


                        template <typename Rhs, typename Dest>
                        void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
                        {
                                m_iterations = Base::maxIterations();
                                m_error = Base::m_tolerance;

                                bool ret = internal::idrs(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, m_S,m_smoothing,m_angle,m_residual);

                                m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence;
                        }

                        void setS(Index S)
                        {
                                if (S < 1)
                                {
                                        S = 4;
                                }

                                m_S = S;
                        }

                        void setSmoothing(bool smoothing)
                        {
                                m_smoothing=smoothing;
                        }

                        void setAngle(RealScalar angle)
                        {
                                m_angle=angle;
                        }

                        void setResidualUpdate(bool update)
                        {
                                m_residual=update;
                        }

        };

}  // namespace Eigen

#endif /* EIGEN_IDRS_H */