// 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
	{
		/**     \internal Low-level Induced Dimension Reduction algoritm
		        \param A The matrix A
		        \param b The right hand side vector b
		        \param x On input and initial solution, on output the computed solution.
		        \param precond A preconditioner being able to efficiently solve for an
		                  approximation of Ax=b (regardless of b)
		        \param iter On input the max number of iteration, on output the number of performed iterations.
		        \param relres On input the tolerance error, on output an estimation of the relative error.
		        \param S On input Number of the dimension of the shadow space.
				\param smoothing switches residual smoothing on.
				\param angle small omega lead to faster convergence at the expense of numerical stability
				\param replacement switches on a residual replacement strategy to increase accuracy of residual at the expense of more Mat*vec products
		        \return false in the case of numerical issue, for example a break down of IDRS.
		*/
		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


/** \ingroup IterativeLinearSolvers_Module
  * \brief The Induced Dimension Reduction method (IDR(s)) is a short-recurrences Krylov method for sparse square problems.
  *
  * This class allows to solve for A.x = b sparse linear problems. The vectors x and b can be either dense or sparse.
  * he Induced Dimension Reduction method, IDR(), is a robust and efficient short-recurrence Krylov subspace method for
  * solving large nonsymmetric systems of linear equations.
  *
  * For indefinite systems IDR(S) outperforms both BiCGStab and BiCGStab(L). Additionally, IDR(S) can handle matrices
  * with complex eigenvalues more efficiently than BiCGStab.
  *
  * Many problems that do not converge for BiCGSTAB converge for IDR(s) (for larger values of s). And if both methods 
  * converge the convergence for IDR(s) is typically much faster for difficult systems (for example indefinite problems). 
  *
  * IDR(s) is a limited memory finite termination method. In exact arithmetic it converges in at most N+N/s iterations,
  * with N the system size.  It uses a fixed number of 4+3s vector. In comparison, BiCGSTAB terminates in 2N iterations 
  * and uses 7 vectors. GMRES terminates in at most N iterations, and uses I+3 vectors, with I the number of iterations. 
  * Restarting GMRES limits the memory consumption, but destroys the finite termination property.
  *
  * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
  * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
  *
  * \implsparsesolverconcept
  *
  * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
  * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
  * and NumTraits<Scalar>::epsilon() for the tolerance.
  *
  * The tolerance corresponds to the relative residual error: |Ax-b|/|b|
  *
  * \b Performance: when using sparse matrices, best performance is achied for a row-major sparse matrix format.
  * Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled.
  * See \ref TopicMultiThreading for details.
  *
  * By default the iterations start with x=0 as an initial guess of the solution.
  * One can control the start using the solveWithGuess() method.
  *
  * IDR(s) can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
  *
  * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
  */
	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:
			/** Default constructor. */
			IDRS(): m_S(4), m_smoothing(false), m_angle(RealScalar(0.7)), m_residual(false) {}

			/**     Initialize the solver with matrix \a A for further \c Ax=b solving.

			        This constructor is a shortcut for the default constructor followed
			        by a call to compute().

			        \warning this class stores a reference to the matrix A as well as some
			        precomputed values that depend on it. Therefore, if \a A is changed
			        this class becomes invalid. Call compute() to update it with the new
			        matrix A, or modify a copy of A.
			*/
			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) {}


			/** \internal */
			/**     Loops over the number of columns of b and does the following:
			                1. sets the tolerence and maxIterations
			                2. Calls the function that has the core solver routine
			*/
			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;
			}

			/** Sets the parameter S, indicating the dimension of the shadow space. Default is 4*/
			void setS(Index S)
			{
				if (S < 1)
				{
					S = 4;
				}

				m_S = S;
			}

			/** Switches off and on smoothing.
			Residual smoothing results in monotonically decreasing residual norms at
			the expense of two extra vectors of storage and a few extra vector
			operations. Although monotonic decrease of the residual norms is a
			desirable property, the rate of convergence of the unsmoothed process and
			the smoothed process is basically the same. Default is off */
			void setSmoothing(bool smoothing)
			{
				m_smoothing=smoothing;
			}

			/** The angle must be a real scalar. In IDR(s), a value for the
			iteration parameter omega must be chosen in every s+1th step. The most
			natural choice is to select a value to minimize the norm of the next residual.
			This corresponds to the parameter omega = 0. In practice, this may lead to
			values of omega that are so small that the other iteration parameters
			cannot be computed with sufficient accuracy. In such cases it is better to
			increase the value of omega sufficiently such that a compromise is reached
			between accurate computations and reduction of the residual norm. The
			parameter angle =0.7 (”maintaining the convergence strategy”)
			results in such a compromise. */
			void setAngle(RealScalar angle)
			{
				m_angle=angle;
			}

			/** The parameter replace is a logical that determines whether a
			residual replacement strategy is employed to increase the accuracy of the
			solution. */
			void setResidualUpdate(bool update)
			{
				m_residual=update;
			}

	};

}  // namespace Eigen

#endif /* EIGEN_IDRS_H */