// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {

namespace internal {

template <typename Scalar>
void dogleg(const Matrix<Scalar, Dynamic, Dynamic> &qrfac, const Matrix<Scalar, Dynamic, 1> &diag,
            const Matrix<Scalar, Dynamic, 1> &qtb, Scalar delta, Matrix<Scalar, Dynamic, 1> &x) {
  using std::abs;
  using std::sqrt;

  typedef DenseIndex Index;

  /* Local variables */
  Index i, j;
  Scalar sum, temp, alpha, bnorm;
  Scalar gnorm, qnorm;
  Scalar sgnorm;

  /* Function Body */
  const Scalar epsmch = NumTraits<Scalar>::epsilon();
  const Index n = qrfac.cols();
  eigen_assert(n == qtb.size());
  eigen_assert(n == x.size());
  eigen_assert(n == diag.size());
  Matrix<Scalar, Dynamic, 1> wa1(n), wa2(n);

  /* first, calculate the gauss-newton direction. */
  for (j = n - 1; j >= 0; --j) {
    temp = qrfac(j, j);
    if (temp == 0.) {
      temp = epsmch * qrfac.col(j).head(j + 1).maxCoeff();
      if (temp == 0.) temp = epsmch;
    }
    if (j == n - 1)
      x[j] = qtb[j] / temp;
    else
      x[j] = (qtb[j] - qrfac.row(j).tail(n - j - 1).dot(x.tail(n - j - 1))) / temp;
  }

  /* test whether the gauss-newton direction is acceptable. */
  qnorm = diag.cwiseProduct(x).stableNorm();
  if (qnorm <= delta) return;

  // TODO : this path is not tested by Eigen unit tests

  /* the gauss-newton direction is not acceptable. */
  /* next, calculate the scaled gradient direction. */

  wa1.fill(0.);
  for (j = 0; j < n; ++j) {
    wa1.tail(n - j) += qrfac.row(j).tail(n - j) * qtb[j];
    wa1[j] /= diag[j];
  }

  /* calculate the norm of the scaled gradient and test for */
  /* the special case in which the scaled gradient is zero. */
  gnorm = wa1.stableNorm();
  sgnorm = 0.;
  alpha = delta / qnorm;
  if (gnorm == 0.) goto algo_end;

  /* calculate the point along the scaled gradient */
  /* at which the quadratic is minimized. */
  wa1.array() /= (diag * gnorm).array();
  // TODO : once unit tests cover this part,:
  // wa2 = qrfac.template triangularView<Upper>() * wa1;
  for (j = 0; j < n; ++j) {
    sum = 0.;
    for (i = j; i < n; ++i) {
      sum += qrfac(j, i) * wa1[i];
    }
    wa2[j] = sum;
  }
  temp = wa2.stableNorm();
  sgnorm = gnorm / temp / temp;

  /* test whether the scaled gradient direction is acceptable. */
  alpha = 0.;
  if (sgnorm >= delta) goto algo_end;

  /* the scaled gradient direction is not acceptable. */
  /* finally, calculate the point along the dogleg */
  /* at which the quadratic is minimized. */
  bnorm = qtb.stableNorm();
  temp = bnorm / gnorm * (bnorm / qnorm) * (sgnorm / delta);
  temp = temp - delta / qnorm * numext::abs2(sgnorm / delta) +
         sqrt(numext::abs2(temp - delta / qnorm) +
              (1. - numext::abs2(delta / qnorm)) * (1. - numext::abs2(sgnorm / delta)));
  alpha = delta / qnorm * (1. - numext::abs2(sgnorm / delta)) / temp;
algo_end:

  /* form appropriate convex combination of the gauss-newton */
  /* direction and the scaled gradient direction. */
  temp = (1. - alpha) * (std::min)(sgnorm, delta);
  x = temp * wa1 + alpha * x;
}

}  // end namespace internal

}  // end namespace Eigen