diff --git a/unsupported/Eigen/IterativeSolvers b/unsupported/Eigen/IterativeSolvers
index bea96af2c..ad7a0be27 100644
--- a/unsupported/Eigen/IterativeSolvers
+++ b/unsupported/Eigen/IterativeSolvers
@@ -22,6 +22,7 @@
   *  - a constrained conjugate gradient
   *  - a Householder GMRES implementation
   *  - an IDR(s) implementation
+  *  - a BiCGSTAB(L) implementation
   *  - a DGMRES implementation
   *  - a MINRES implementation
   *
@@ -87,6 +88,7 @@
 #include "src/IterativeSolvers/DGMRES.h"
 #include "src/IterativeSolvers/MINRES.h"
 #include "src/IterativeSolvers/IDRS.h"
+#include "src/IterativeSolvers/BiCGSTABL.h"
 
 #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
 
diff --git a/unsupported/Eigen/src/IterativeSolvers/BiCGSTABL.h b/unsupported/Eigen/src/IterativeSolvers/BiCGSTABL.h
new file mode 100755
index 000000000..141d7056d
--- /dev/null
+++ b/unsupported/Eigen/src/IterativeSolvers/BiCGSTABL.h
@@ -0,0 +1,339 @@
+// 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>
+// Copyright (C) 2020 Adithya Vijaykumar
+//
+// 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/.
+
+/*
+
+  This implementation of BiCGStab(L) is based on the papers
+      General algorithm:
+      1. G.L.G. Sleijpen, D.R. Fokkema. (1993). BiCGstab(l) for linear equations
+  involving unsymmetric matrices with complex spectrum. Electronic Transactions
+  on Numerical Analysis. Polynomial step update:
+      2. G.L.G. Sleijpen, M.B. Van Gijzen. (2010) Exploiting BiCGstab(l)
+  strategies to induce dimension reduction SIAM Journal on Scientific Computing.
+      3. Fokkema, Diederik R. Enhanced implementation of BiCGstab (l) for
+  solving linear systems of equations. Universiteit Utrecht. Mathematisch
+  Instituut, 1996
+      4. Sleijpen, G. L., & van der Vorst, H. A. (1996). Reliable updated
+  residuals in hybrid Bi-CG methods. Computing, 56(2), 141-163.
+*/
+
+#ifndef EIGEN_BICGSTABL_H
+#define EIGEN_BICGSTABL_H
+
+namespace Eigen {
+
+namespace internal {
+/**     \internal Low-level bi conjugate gradient stabilized algorithm with L
+   additional residual minimization steps \param mat The matrix A \param rhs 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 iters
+   On input the max number of iteration, on output the number of performed
+   iterations. \param tol_error On input the tolerance error, on output an
+   estimation of the relative error. \param L On input Number of additional
+   GMRES steps to take. If L is too large (~20) instabilities occur. \return
+   false in the case of numerical issue, for example a break down of BiCGSTABL.
+*/
+template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
+bool bicgstabl(const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond, Index &iters,
+               typename Dest::RealScalar &tol_error, Index L) {
+  using numext::abs;
+  using numext::sqrt;
+  typedef typename Dest::RealScalar RealScalar;
+  typedef typename Dest::Scalar Scalar;
+  const Index N = rhs.size();
+  L = L < x.rows() ? L : x.rows();
+
+  Index k = 0;
+
+  const RealScalar tol = tol_error;
+  const Index maxIters = iters;
+
+  typedef Matrix<Scalar, Dynamic, 1> VectorType;
+  typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> DenseMatrixType;
+
+  DenseMatrixType rHat(N, L + 1);
+  DenseMatrixType uHat(N, L + 1);
+
+  // We start with an initial guess x_0 and let us set r_0 as (residual
+  // calculated from x_0)
+  VectorType x0 = x;
+  rHat.col(0) = rhs - mat * x0;  // r_0
+
+  x.setZero();  // This will contain the updates to the solution.
+  // rShadow is arbritary, but must never be orthogonal to any residual.
+  VectorType rShadow = VectorType::Random(N);
+
+  VectorType x_prime = x;
+
+  // Redundant: x is already set to 0
+  // x.setZero();
+  VectorType b_prime = rHat.col(0);
+
+  // Other vectors and scalars initialization
+  Scalar rho0 = 1.0;
+  Scalar alpha = 0.0;
+  Scalar omega = 1.0;
+
+  uHat.col(0).setZero();
+
+  bool bicg_convergence = false;
+
+  const RealScalar normb = rhs.stableNorm();
+  if (internal::isApprox(normb, RealScalar(0))) {
+    x.setZero();
+    iters = 0;
+    return true;
+  }
+  RealScalar normr = rHat.col(0).stableNorm();
+  RealScalar Mx = normr;
+  RealScalar Mr = normr;
+
+  // Keep track of the solution with the lowest residual
+  RealScalar normr_min = normr;
+  VectorType x_min = x_prime + x;
+
+  // Criterion for when to apply the group-wise update, conform ref 3.
+  const RealScalar delta = 0.01;
+
+  bool compute_res = false;
+  bool update_app = false;
+
+  while (normr > tol * normb && k < maxIters) {
+    rho0 *= -omega;
+
+    for (Index j = 0; j < L; ++j) {
+      const Scalar rho1 = rShadow.dot(rHat.col(j));
+
+      if (!(numext::isfinite)(rho1) || rho0 == RealScalar(0.0)) {
+        // We cannot continue computing, return the best solution found.
+        x += x_prime;
+
+        // Check if x is better than the best stored solution thus far.
+        normr = (rhs - mat * (precond.solve(x) + x0)).stableNorm();
+
+        if (normr > normr_min || !(numext::isfinite)(normr)) {
+          // x_min is a better solution than x, return x_min
+          x = x_min;
+          normr = normr_min;
+        }
+        tol_error = normr / normb;
+        iters = k;
+        // x contains the updates to x0, add those back to obtain the solution
+        x = precond.solve(x);
+        x += x0;
+        return (normr < tol * normb);
+      }
+
+      const Scalar beta = alpha * (rho1 / rho0);
+      rho0 = rho1;
+      // Update search directions
+      uHat.leftCols(j + 1) = rHat.leftCols(j + 1) - beta * uHat.leftCols(j + 1);
+      uHat.col(j + 1) = mat * precond.solve(uHat.col(j));
+      const Scalar sigma = rShadow.dot(uHat.col(j + 1));
+      alpha = rho1 / sigma;
+      // Update residuals
+      rHat.leftCols(j + 1) -= alpha * uHat.middleCols(1, j + 1);
+      rHat.col(j + 1) = mat * precond.solve(rHat.col(j));
+      // Complete BiCG iteration by updating x
+      x += alpha * uHat.col(0);
+      normr = rHat.col(0).stableNorm();
+      // Check for early exit
+      if (normr < tol * normb) {
+        /*
+          Convergence was achieved during BiCG step.
+          Without this check BiCGStab(L) fails for trivial matrices, such as
+          when the preconditioner already is the inverse, or the input matrix is
+          identity.
+        */
+        bicg_convergence = true;
+        break;
+      } else if (normr < normr_min) {
+        // We found an x with lower residual, keep this one.
+        x_min = x + x_prime;
+        normr_min = normr;
+      }
+    }
+    if (!bicg_convergence) {
+      /*
+        The polynomial/minimize residual step.
+
+        QR Householder method for argmin is more stable than (modified)
+        Gram-Schmidt, in the sense that there is less loss of orthogonality. It
+        is more accurate than solving the normal equations, since the normal
+        equations scale with condition number squared.
+      */
+      const VectorType gamma = rHat.rightCols(L).householderQr().solve(rHat.col(0));
+      x += rHat.leftCols(L) * gamma;
+      rHat.col(0) -= rHat.rightCols(L) * gamma;
+      uHat.col(0) -= uHat.rightCols(L) * gamma;
+      normr = rHat.col(0).stableNorm();
+      omega = gamma(L - 1);
+    }
+    if (normr < normr_min) {
+      // We found an x with lower residual, keep this one.
+      x_min = x + x_prime;
+      normr_min = normr;
+    }
+
+    k++;
+
+    /*
+      Reliable update part
+
+      The recursively computed residual can deviate from the actual residual
+      after several iterations. However, computing the residual from the
+      definition costs extra MVs and should not be done at each iteration. The
+      reliable update strategy computes the true residual from the definition:
+      r=b-A*x at strategic intervals. Furthermore a "group wise update" strategy
+      is used to combine updates, which improves accuracy.
+    */
+
+    // Maximum norm of residuals since last update of x.
+    Mx = numext::maxi(Mx, normr);
+    // Maximum norm of residuals since last computation of the true residual.
+    Mr = numext::maxi(Mr, normr);
+
+    if (normr < delta * normb && normb <= Mx) {
+      update_app = true;
+    }
+
+    if (update_app || (normr < delta * Mr && normb <= Mr)) {
+      compute_res = true;
+    }
+
+    if (bicg_convergence) {
+      update_app = true;
+      compute_res = true;
+      bicg_convergence = false;
+    }
+
+    if (compute_res) {
+      // Explicitly compute residual from the definition
+
+      // This is equivalent to the shifted version of rhs - mat *
+      // (precond.solve(x)+x0)
+      rHat.col(0) = b_prime - mat * precond.solve(x);
+      normr = rHat.col(0).stableNorm();
+      Mr = normr;
+
+      if (update_app) {
+        // After the group wise update, the original problem is translated to a
+        // shifted one.
+        x_prime += x;
+        x.setZero();
+        b_prime = rHat.col(0);
+        Mx = normr;
+      }
+    }
+    if (normr < normr_min) {
+      // We found an x with lower residual, keep this one.
+      x_min = x + x_prime;
+      normr_min = normr;
+    }
+
+    compute_res = false;
+    update_app = false;
+  }
+
+  // Convert internal variable to the true solution vector x
+  x += x_prime;
+
+  normr = (rhs - mat * (precond.solve(x) + x0)).stableNorm();
+  if (normr > normr_min || !(numext::isfinite)(normr)) {
+    // x_min is a better solution than x, return x_min
+    x = x_min;
+    normr = normr_min;
+  }
+  tol_error = normr / normb;
+  iters = k;
+
+  // x contains the updates to x0, add those back to obtain the solution
+  x = precond.solve(x);
+  x += x0;
+  return true;
+}
+
+}  // namespace internal
+
+template <typename MatrixType_, typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar>>
+class BiCGSTABL;
+
+namespace internal {
+
+template <typename MatrixType_, typename Preconditioner_>
+struct traits<Eigen::BiCGSTABL<MatrixType_, Preconditioner_>> {
+  typedef MatrixType_ MatrixType;
+  typedef Preconditioner_ Preconditioner;
+};
+
+}  // namespace internal
+
+template <typename MatrixType_, typename Preconditioner_>
+class BiCGSTABL : public IterativeSolverBase<BiCGSTABL<MatrixType_, Preconditioner_>> {
+  typedef IterativeSolverBase<BiCGSTABL> Base;
+  using Base::m_error;
+  using Base::m_info;
+  using Base::m_isInitialized;
+  using Base::m_iterations;
+  using Base::matrix;
+  Index m_L;
+
+ public:
+  typedef MatrixType_ MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef Preconditioner_ Preconditioner;
+
+  /** Default constructor. */
+  BiCGSTABL() : m_L(2) {}
+
+  /**
+  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 BiCGSTABL(const EigenBase<MatrixDerived> &A) : Base(A.derived()), m_L(2) {}
+
+  /** \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::bicgstabl(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, m_L);
+    m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence;
+  }
+
+  /** Sets the parameter L, indicating how many minimize residual steps are
+   * used. Default: 2 */
+  void setL(Index L) {
+    eigen_assert(L >= 1 && "L needs to be positive");
+    m_L = L;
+  }
+};
+
+}  // namespace Eigen
+
+#endif /* EIGEN_BICGSTABL_H */
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index f87bacd48..6bd8430cb 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -99,6 +99,7 @@ ei_add_test(gmres)
 ei_add_test(dgmres)
 ei_add_test(minres)
 ei_add_test(idrs)
+ei_add_test(bicgstabl)
 ei_add_test(levenberg_marquardt)
 ei_add_test(kronecker_product)
 ei_add_test(bessel_functions)
diff --git a/unsupported/test/bicgstabl.cpp b/unsupported/test/bicgstabl.cpp
new file mode 100644
index 000000000..302848c98
--- /dev/null
+++ b/unsupported/test/bicgstabl.cpp
@@ -0,0 +1,31 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2012 Kolja Brix <brix@igpm.rwth-aaachen.de>
+//
+// 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/.
+
+#include "../../test/sparse_solver.h"
+#include <Eigen/IterativeSolvers>
+
+template<typename T> void test_bicgstabl_T()
+{
+  BiCGSTABL<SparseMatrix<T>, DiagonalPreconditioner<T> > bicgstabl_colmajor_diag;
+  BiCGSTABL<SparseMatrix<T>, IncompleteLUT<T> >           bicgstabl_colmajor_ilut;
+
+  //This does not change the tolerance of the test, only the tolerance of the solver.
+  bicgstabl_colmajor_diag.setTolerance(NumTraits<T>::epsilon()*20);
+  bicgstabl_colmajor_ilut.setTolerance(NumTraits<T>::epsilon()*20);
+
+  CALL_SUBTEST( check_sparse_square_solving(bicgstabl_colmajor_diag)  );
+  CALL_SUBTEST( check_sparse_square_solving(bicgstabl_colmajor_ilut)     );
+}
+
+EIGEN_DECLARE_TEST(bicgstabl)
+{
+  CALL_SUBTEST_1(test_bicgstabl_T<double>());
+  CALL_SUBTEST_2(test_bicgstabl_T<std::complex<double> >());
+}