// ---------------------------------------------------------------------
//
// Copyright (c) 2017-2018 The Regents of the University of Michigan and DFT-FE
// authors.
//
// This file is part of the DFT-FE code.
//
// The DFT-FE code 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 at
// the top level of the DFT-FE distribution.
//
// ---------------------------------------------------------------------
//
// @author Bikash Kanungo, Vishal Subramanian
//

#include <BFGSInverseDFTSolver.h>

namespace invDFT {
namespace {
template <typename T>
std::string to_string_with_precision(const T a_value, const int n = 6) {
  std::ostringstream out;
  out.precision(n);
  out << std::fixed << a_value;
  return out.str();
}
//
// y = a*x + b*y
//
void vecAdd(const dftfe::distributedCPUVec<double> &x,
            dftfe::distributedCPUVec<double> &y, const double a,
            const double b) {
  const unsigned int N = x.locally_owned_size();
  for (unsigned int i = 0; i < N; ++i)
    y.local_element(i) = a * x.local_element(i) + b * y.local_element(i);
}

void vecScale(dftfe::distributedCPUVec<double> &x, const double a) {
  const unsigned int N = x.locally_owned_size();
  for (unsigned int i = 0; i < N; ++i)
    x.local_element(i) *= a;
}

} // namespace

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::
    BFGSInverseDFTSolver(int numComponents, double tol, double lineSearchTol,
                         unsigned int maxNumIter, int historySize,
                         int numLineSearch, const MPI_Comm &mpi_comm_parent)
    : d_numComponents(numComponents), d_tol(tol),
      d_lineSearchTol(lineSearchTol), d_maxNumIter(maxNumIter),
      d_historySize(historySize), d_numLineSearch(numLineSearch),
      d_k(numComponents, 0),
      d_y(numComponents, std::list<dftfe::distributedCPUVec<double>>(0)),
      d_s(numComponents, std::list<dftfe::distributedCPUVec<double>>(0)),
      d_rho(numComponents, std::list<double>(0)),
      pcout(std::cout,
            (dealii::Utilities::MPI::this_mpi_process(mpi_comm_parent) == 0)) {}

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
void BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::
    inverseJacobianTimesVec(
        const dftfe::distributedCPUVec<double> &g,
        dftfe::distributedCPUVec<double> &z, const unsigned int component,
        InverseDFTSolverFunction<FEOrder, FEOrderElectro, memorySpace>
            &iDFTSolverFunction) {
  int N = d_k[component];
  dftfe::distributedCPUVec<double> q = g;
  auto itReverseY = d_y[component].rbegin();
  auto itReverseS = d_s[component].rbegin();
  auto itReverseRho = d_rho[component].rbegin();
  double gamma = 1.0;
  if (itReverseY != d_y[component].rend()) {
    // gamma = (s_k^T*y_k)/(y_k^Ty_K)
    std::vector<double> dotProd1(1, 0.0);
    std::vector<double> dotProd2(1, 0.0);
    iDFTSolverFunction.dotProduct(*itReverseS, *itReverseY, 1, dotProd1);
    iDFTSolverFunction.dotProduct(*itReverseY, *itReverseY, 1, dotProd2);
    gamma = dotProd1[0] / dotProd2[0];
  }

  std::vector<double> alpha(N, 0.0);
  std::vector<double> beta(N, 0.0);
  for (int i = N - 1; i >= 0; --i) {
    std::vector<double> dotProd(1, 0.0);
    iDFTSolverFunction.dotProduct(*itReverseS, q, 1, dotProd);
    alpha[i] = (*itReverseRho) * dotProd[0];
    // q = q - alpha*y
    vecAdd(*itReverseY, q, -alpha[i], 1.0);
    ++itReverseY;
    ++itReverseS;
    ++itReverseRho;
  }

  z = q;
  vecScale(z, gamma);

  auto itY = d_y[component].begin();
  auto itS = d_s[component].begin();
  auto itRho = d_rho[component].begin();
  for (int i = 0; i < N; ++i) {
    std::vector<double> dotProd(1, 0.0);
    iDFTSolverFunction.dotProduct(*itY, z, 1, dotProd);
    beta[i] = (*itRho) * dotProd[0];
    // z = z + (alpha-beta)*s
    vecAdd(*itS, z, alpha[i] - beta[i], 1.0);
    ++itY;
    ++itS;
    ++itRho;
  }
}

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
void BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::fnormLoss(
    const std::vector<dftfe::distributedCPUVec<double>> &x,
    const std::vector<dftfe::distributedCPUVec<double>> &p,
    const std::vector<double> &alpha, std::vector<double> &fnorms,
    InverseDFTSolverFunction<FEOrder, FEOrderElectro, memorySpace>
        &iDFTSolverFunction) {
  std::vector<dftfe::distributedCPUVec<double>> xnew(d_numComponents);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    xnew[iComp].reinit(x[iComp], false);
    // xnew = x + alpha*p
    xnew[iComp] = x[iComp];
    vecAdd(p[iComp], xnew[iComp], alpha[iComp], 1.0);
  }

  std::vector<dftfe::distributedCPUVec<double>> g(d_numComponents);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    g[iComp].reinit(x[iComp], false);
    g[iComp] = 0.0;
  }

  std::vector<double> L(d_numComponents);
  iDFTSolverFunction.getForceVector(xnew, g, L);
  fnorms.resize(d_numComponents, 0.0);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    fnorms[iComp] = L[iComp];
  }
}

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
void BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::fnormGrad(
    const std::vector<dftfe::distributedCPUVec<double>> &x,
    const std::vector<dftfe::distributedCPUVec<double>> &p,
    const std::vector<double> &alpha, std::vector<double> &fnorms,
    InverseDFTSolverFunction<FEOrder, FEOrderElectro, memorySpace>
        &iDFTSolverFunction) {
  std::vector<dftfe::distributedCPUVec<double>> xnew(d_numComponents);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    xnew[iComp].reinit(x[iComp], false);
    // xnew = x + alpha*p
    xnew[iComp] = x[iComp];
    vecAdd(p[iComp], xnew[iComp], alpha[iComp], 1.0);
  }

  std::vector<dftfe::distributedCPUVec<double>> g(d_numComponents);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    g[iComp].reinit(x[iComp], false);
    g[iComp] = 0.0;
  }

  std::vector<double> L(d_numComponents);
  double constraint;
  iDFTSolverFunction.getForceVector(xnew, g, L);
  fnorms.resize(d_numComponents, 0.0);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    std::vector<double> dotProd(1, 0.0);
    iDFTSolverFunction.dotProduct(g[iComp], g[iComp], 1, dotProd);
    fnorms[iComp] = std::sqrt(dotProd[0]);
  }
}

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
void BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::fnormCP(
    const std::vector<dftfe::distributedCPUVec<double>> &x,
    const std::vector<dftfe::distributedCPUVec<double>> &p,
    const std::vector<double> &alpha, std::vector<double> &fnorms,
    InverseDFTSolverFunction<FEOrder, FEOrderElectro, memorySpace>
        &iDFTSolverFunction) {
  std::vector<dftfe::distributedCPUVec<double>> xnew(d_numComponents);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    xnew[iComp].reinit(x[iComp], false);
    // xnew = x + alpha*p
    xnew[iComp] = x[iComp];
    vecAdd(p[iComp], xnew[iComp], alpha[iComp], 1.0);
  }

  std::vector<dftfe::distributedCPUVec<double>> g(d_numComponents);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    g[iComp].reinit(x[iComp], false);
    g[iComp] = 0.0;
  }

  std::vector<double> L(d_numComponents);
  iDFTSolverFunction.getForceVector(xnew, g, L);
  fnorms.resize(d_numComponents, 0.0);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    std::vector<double> dotProd(1, 0.0);
    iDFTSolverFunction.dotProduct(g[iComp], d_p[iComp], 1, dotProd);
    fnorms[iComp] = dotProd[0];
  }
}

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
void BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::
    solveLineSearchCP(
        std::vector<std::vector<double>> &lambda,
        std::vector<std::vector<double>> &f, const int maxIter,
        const double tolerance,
        InverseDFTSolverFunction<FEOrder, FEOrderElectro, memorySpace>
            &iDFTSolverFunction) {
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    dftfe::utils::throwException(
        lambda[iComp].size() >= 2,
        "At least two initial values are need for a secant method.");
  }

  std::vector<double> lambda0(d_numComponents);
  std::vector<double> lambda1(d_numComponents);
  std::vector<double> f0(d_numComponents);
  std::vector<double> f1(d_numComponents);
  for (unsigned int i = 0; i < maxIter; ++i) {
    int N = lambda[0].size();
    int M = f[0].size();
    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      lambda0[iComp] = lambda[iComp][N - 2];
      lambda1[iComp] = lambda[iComp][N - 1];
    }
    if (M > 0) {
      for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp)
        f0[iComp] = f[iComp][M - 1];
    } else {
      fnormCP(d_x, d_p, lambda0, f0, iDFTSolverFunction);
      for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp)
        f[iComp].push_back(f0[iComp]);
    }

    if (std::all_of(f0.begin(), f0.end(), [tolerance](double x) {
          return std::sqrt(std::fabs(x)) < tolerance;
        })) {
      // remove the last element (i.e., lambda1)
      for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp)
        lambda[iComp].pop_back();

      std::cout << "f0 in secantCP below tolerance" << std::endl;
      break;
    }

    this->fnormCP(d_x, d_p, lambda1, f1, iDFTSolverFunction);
    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp)
      f[iComp].push_back(f1[iComp]);

    if (std::all_of(f1.begin(), f1.end(), [tolerance](double x) {
          return std::sqrt(std::fabs(x)) < tolerance;
        })) {
      pcout << "f1 in secantCP below tolerance" << std::endl;
      break;
    }

    //
    // TODO Fetch the tolerance from dftParams
    //
    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      if (fabs((f1[iComp] - f0[iComp]) / f1[iComp]) < 1e-15) {
        std::string message =
            "Secant line search failed, possibly because f'(x) = 0 in the "
            "interval"
            "The last two alphas and their function values are:\n"
            " Alphas: (" +
            to_string_with_precision(lambda0[iComp], 18) + "," +
            to_string_with_precision(lambda1[iComp], 18) + ")\t" + "fs: (" +
            to_string_with_precision(f0[iComp], 18) + "," +
            to_string_with_precision(f1[iComp], 18) + ").";
        dftfe::utils::throwException(false, message);
      }
    }

    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      double s = (f1[iComp] - f0[iComp]) / (lambda1[iComp] - lambda0[iComp]);
      /* if the solve is going in the wrong direction, reverse it */
      if (s > 0.0)
        s = -s;
      double lambdaNext = lambda1[iComp] - f1[iComp] / s;
      // switch directions if we stepped out of bounds
      if (lambdaNext < 0.0)
        lambdaNext = lambda1[iComp] + f1[iComp] / s;

      lambda[iComp].push_back(lambdaNext);
    }
  }
}

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
void BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::
    solveLineSearchSecantLoss(
        std::vector<std::vector<double>> &lambda,
        std::vector<std::vector<double>> &f, const int maxIter,
        const double tolerance,
        InverseDFTSolverFunction<FEOrder, FEOrderElectro, memorySpace>
            &iDFTSolverFunction) {
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    if (lambda[iComp].size() < 2)
      dftfe::utils::throwException(
          false, "At least two initial values are need for a secant method.");
  }

  std::vector<double> lambdaOld(d_numComponents);
  std::vector<double> lambdaMid(d_numComponents);
  std::vector<double> lambdaNew(d_numComponents);
  std::vector<double> fOld(d_numComponents);
  std::vector<double> fMid(d_numComponents);
  std::vector<double> fNew(d_numComponents);

  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    lambdaOld[iComp] = lambda[iComp][0];
    lambdaNew[iComp] = lambda[iComp][1];
    lambdaMid[iComp] = 0.5 * (lambdaNew[iComp] + lambdaOld[iComp]);
    fOld[iComp] = f[iComp][0];
  }

  for (int i = 0; i < maxIter; i++) {
    /* compute the objective at the midpoint */
    this->fnormLoss(d_x, d_p, lambdaMid, fMid, iDFTSolverFunction);
    /* compute the objective at the new endpoint */
    this->fnormLoss(d_x, d_p, lambdaNew, fNew, iDFTSolverFunction);

    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      pcout << "LineSearch L2: " << i << " for component: " << iComp
            << std::endl;
      pcout << lambdaOld[iComp] << " " << fOld[iComp] << std::endl;
      pcout << lambdaMid[iComp] << " " << fMid[iComp] << std::endl;
      pcout << lambdaNew[iComp] << " " << fNew[iComp] << std::endl;
    }
    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      lambdaNew[iComp] = .5 * (lambdaNew[iComp] + lambdaOld[iComp]);
      lambdaMid[iComp] = .5 * (lambdaNew[iComp] + lambdaOld[iComp]);
      const double delLambda = lambdaNew[iComp] - lambdaOld[iComp];
      /* compute f'() at the end points using second order one sided
       * differencing */
      const double delF =
          (3. * fNew[iComp] - 4. * fMid[iComp] + 1. * fOld[iComp]) / delLambda;
      const double delFOld =
          (-3. * fOld[iComp] + 4. * fMid[iComp] - 1. * fNew[iComp]) / delLambda;
      /* compute f''() at the midpoint using centered differencing */
      const double del2F = (delF - delFOld) / delLambda;

      double lambdaUpdate = 0.0;
      /* compute the secant (Newton) update -- always go downhill */
      if (del2F > 0.)
        lambdaUpdate = lambdaNew[iComp] - delF / del2F;
      else if (del2F < 0.)
        lambdaUpdate = lambdaNew[iComp] + delF / del2F;
      else
        break;

      if (lambdaUpdate < 0.0)
        lambdaUpdate = 0.5 * (lambdaNew[iComp] + lambdaOld[iComp]);

      /* update the endpoints and the midpoint of the bracketed secant region */
      lambdaOld[iComp] = lambdaNew[iComp];
      lambdaNew[iComp] = lambdaUpdate;
      fOld[iComp] = fNew[iComp];
      lambdaMid[iComp] = 0.5 * (lambdaNew[iComp] + lambdaOld[iComp]);
      lambda[iComp].push_back(lambdaNew[iComp]);
      f[iComp].push_back(fNew[iComp]);
      pcout << "Updated lambda for component: " << iComp << ": "
            << lambdaNew[iComp] << std::endl;
    }
  }
}

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
void BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::
    solveLineSearchSecantForceNorm(
        std::vector<std::vector<double>> &lambda,
        std::vector<std::vector<double>> &f, const int maxIter,
        const double tolerance,
        InverseDFTSolverFunction<FEOrder, FEOrderElectro, memorySpace>
            &iDFTSolverFunction) {
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    if (lambda[iComp].size() < 2)
      dftfe::utils::throwException(
          false, "At least two initial values are need for a secant method.");
  }

  std::vector<double> lambdaOld(d_numComponents);
  std::vector<double> lambdaMid(d_numComponents);
  std::vector<double> lambdaNew(d_numComponents);
  std::vector<double> fOld(d_numComponents);
  std::vector<double> fMid(d_numComponents);
  std::vector<double> fNew(d_numComponents);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    lambdaOld[iComp] = lambda[iComp][0];
    lambdaNew[iComp] = lambda[iComp][1];
    lambdaMid[iComp] = 0.5 * (lambdaNew[iComp] + lambdaOld[iComp]);
    fOld[iComp] = f[iComp][0];
  }
  for (int i = 0; i < maxIter; i++) {

    /* compute the objective at the midpoint */
    this->fnormGrad(d_x, d_p, lambdaMid, fMid, iDFTSolverFunction);
    /* compute the objective at the new endpoint */
    this->fnormGrad(d_x, d_p, lambdaNew, fNew, iDFTSolverFunction);

    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      pcout << "LineSearch L2: " << i << " for component: " << iComp
            << std::endl;
      pcout << lambdaOld[iComp] << " " << fOld[iComp] << std::endl;
      pcout << lambdaMid[iComp] << " " << fMid[iComp] << std::endl;
      pcout << lambdaNew[iComp] << " " << fNew[iComp] << std::endl;
    }

    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      lambdaNew[iComp] = .5 * (lambdaNew[iComp] + lambdaOld[iComp]);
      lambdaMid[iComp] = .5 * (lambdaNew[iComp] + lambdaOld[iComp]);
      const double delLambda = lambdaNew[iComp] - lambdaOld[iComp];
      /* compute f'() at the end points using second order one sided
       * differencing */
      const double delF =
          (3. * fNew[iComp] - 4. * fMid[iComp] + 1. * fOld[iComp]) / delLambda;
      const double delFOld =
          (-3. * fOld[iComp] + 4. * fMid[iComp] - 1. * fNew[iComp]) / delLambda;
      /* compute f''() at the midpoint using centered differencing */
      const double del2F = (delF - delFOld) / delLambda;

      double lambdaUpdate = 0.0;
      /* compute the secant (Newton) update -- always go downhill */
      if (del2F > 0.)
        lambdaUpdate = lambdaNew[iComp] - delF / del2F;
      else if (del2F < 0.)
        lambdaUpdate = lambdaNew[iComp] + delF / del2F;
      else
        break;

      if (lambdaUpdate < 0.0)
        lambdaUpdate = 0.5 * (lambdaNew[iComp] + lambdaOld[iComp]);

      /* update the endpoints and the midpoint of the bracketed secant region */
      lambdaOld[iComp] = lambdaNew[iComp];
      lambdaNew[iComp] = lambdaUpdate;
      fOld[iComp] = fNew[iComp];
      lambdaMid[iComp] = 0.5 * (lambdaNew[iComp] + lambdaOld[iComp]);
      lambda[iComp].push_back(lambdaNew[iComp]);
      f[iComp].push_back(fNew[iComp]);
      pcout << "Updated lambda for component: " << iComp << ": "
            << lambdaNew[iComp] << std::endl;
    }
  }
}

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
void BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::solve(
    InverseDFTSolverFunction<FEOrder, FEOrderElectro, memorySpace>
        &iDFTSolverFunction,
    const BFGSInverseDFTSolver::LSType lsType) {
  std::vector<dftfe::distributedCPUVec<double>> y(d_numComponents);
  std::vector<dftfe::distributedCPUVec<double>> s(d_numComponents);

  d_x = iDFTSolverFunction.getInitialGuess();

  //
  // allocate d_g, d_p, y, and s to be of the same size as d_x
  // and initialize them to zero
  //
  d_g.resize(d_numComponents);
  d_p.resize(d_numComponents);
  for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
    d_g[iComp].reinit(d_x[iComp], false);
    d_p[iComp].reinit(d_x[iComp], false);
    y[iComp].reinit(d_x[iComp], false);
    s[iComp].reinit(d_x[iComp], false);
    d_g[iComp] = 0.0;
    d_p[iComp] = 0.0;
    y[iComp] = 0.0;
    s[iComp] = 0.0;
  }

  std::vector<double> L(d_numComponents);
  for (unsigned int iter = 0; iter < d_maxNumIter; ++iter) {
    pcout << " bfgs iter = " << iter << "\n";
    if (iter == 0) {
      iDFTSolverFunction.getForceVector(d_x, d_g, L);
    }

    std::vector<double> gnorm(d_numComponents);
    std::vector<std::vector<double>> lsNorm(d_numComponents);
    std::vector<std::vector<double>> lambdas(d_numComponents);
    std::vector<bool> hasConverged(d_numComponents, false);
    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      std::vector<double> dotProd(1, 0.0);
      iDFTSolverFunction.dotProduct(d_g[iComp], d_g[iComp], 1, dotProd);
      const double forceNorm = std::sqrt(dotProd[0]);
      pcout << "BFGS-Inverse " << iter << " forceNorm[" << iComp
            << "]: " << forceNorm << std::endl;

      if (forceNorm < d_tol)
        hasConverged[iComp] = true;

      inverseJacobianTimesVec(d_g[iComp], d_p[iComp], iComp,
                              iDFTSolverFunction);
      vecScale(d_p[iComp], -1.0);

      double f0 = 0.0;
      if (lsType == BFGSInverseDFTSolver::LSType::SECANT_LOSS) {
        f0 = L[iComp];
      } else if (lsType == BFGSInverseDFTSolver::LSType::SECANT_FORCE_NORM) {
        f0 = forceNorm;
      } else if (lsType == BFGSInverseDFTSolver::LSType::CP) {
        std::vector<double> gDotP(1, 0.0);
        iDFTSolverFunction.dotProduct(d_g[iComp], d_p[iComp], 1, gDotP);
        f0 = gDotP[0];
      } else {
        dftfe::utils::throwException(
            false,
            "Invalid line search type passed to BFGSInverseDFTSolver.solve()");
      }

      lambdas[iComp].push_back(0.0);
      lambdas[iComp].push_back(1.0);
      lsNorm[iComp].push_back(f0);
    }

    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      pcout << " printing the loss val for iComp = " << iComp
            << " loss = " << L[iComp] << "\n";
    }

    if (std::all_of(hasConverged.begin(), hasConverged.end(),
                    [](bool x) { return x; })) {
      break;
    }

    if (lsType == BFGSInverseDFTSolver::LSType::SECANT_LOSS) {
      this->solveLineSearchSecantLoss(lambdas, lsNorm, d_numLineSearch,
                                      d_lineSearchTol, iDFTSolverFunction);
    } else if (lsType == BFGSInverseDFTSolver::LSType::SECANT_FORCE_NORM) {
      this->solveLineSearchSecantForceNorm(lambdas, lsNorm, d_numLineSearch,
                                           d_lineSearchTol, iDFTSolverFunction);
    } else {
      this->solveLineSearchCP(lambdas, lsNorm, d_numLineSearch, d_lineSearchTol,
                              iDFTSolverFunction);
    }

    double optimalLambda = 1.0;
    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      const unsigned int numLambdas = lambdas[iComp].size();
      pcout << "Line search for component" << iComp << std::endl;
      for (unsigned int i = 0; i < numLambdas - 1; ++i) {
        pcout << "Lambda: " << lambdas[iComp][i]
              << " Norm: " << lsNorm[iComp][i] << std::endl;
      }

      //
      // If the lambda is negative set lambda to 1.0, else it will lead to
      // ascent instead of descent along the gradient
      //
      optimalLambda = (lambdas[iComp][numLambdas - 1] > 0.0)
                          ? lambdas[iComp][numLambdas - 1]
                          : 1.0;
      pcout << "Optimal lambda for iComp " << iComp << " is " << optimalLambda
            << std::endl;

      double alpha = optimalLambda;

      // s = alpha*d_p
      s[iComp] = d_p[iComp];
      vecScale(s[iComp], alpha);

      // d_x = s + d_x
      vecAdd(s[iComp], d_x[iComp], 1.0, 1.0);

      // y = g(x_{k+1}) - g(x_k)
      // first set y = -d_g (i.e. y = -g(x_k))
      y[iComp] = d_g[iComp];
      vecScale(y[iComp], -1.0);
    }

    // evaluate g(x_{k+1})
    iDFTSolverFunction.getForceVector(d_x, d_g, L);

    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      // add g(x_{k+1}) to y
      // to obtain y = g(x_{k+1}) - g(x_k)
      vecAdd(d_g[iComp], y[iComp], 1.0, 1.0);

      std::vector<double> curvature(1);
      iDFTSolverFunction.dotProduct(y[iComp], s[iComp], 1, curvature);
      const double rho = curvature[0];
      pcout << "Curvature condition (y^Ts) [" << iComp << "]: " << rho
            << std::endl;

      //
      // add to history, if curvature is positive
      // TODO Use a small finite value instead
      if (rho > 0.0) {
        d_y[iComp].push_back(y[iComp]);
        d_s[iComp].push_back(s[iComp]);
        d_rho[iComp].push_back(1.0 / rho);
        d_k[iComp]++;
      }
    }

    for (unsigned int iComp = 0; iComp < d_numComponents; ++iComp) {
      if (d_k[iComp] > d_historySize) {
        d_k[iComp] = d_historySize;
        d_y[iComp].pop_front();
        d_s[iComp].pop_front();
        d_rho[iComp].pop_front();
      }
    }
  }
}

template <unsigned int FEOrder, unsigned int FEOrderElectro,
          dftfe::utils::MemorySpace memorySpace>
std::vector<dftfe::distributedCPUVec<double>>
BFGSInverseDFTSolver<FEOrder, FEOrderElectro, memorySpace>::getSolution()
    const {
  return d_x;
}

template class BFGSInverseDFTSolver<2, 2, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<2, 3, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<2, 4, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<3, 3, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<3, 4, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<3, 5, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<3, 6, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<4, 4, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<4, 5, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<4, 6, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<4, 7, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<5, 5, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<5, 6, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<5, 7, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<5, 8, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<6, 6, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<6, 7, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<6, 8, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<6, 9, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<7, 7, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<7, 8, dftfe::utils::MemorySpace::HOST>;
template class BFGSInverseDFTSolver<7, 9, dftfe::utils::MemorySpace::HOST>;
#ifdef DFTFE_WITH_DEVICE
template class BFGSInverseDFTSolver<2, 2, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<2, 3, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<2, 4, dftfe::utils::MemorySpace::DEVICE>;

template class BFGSInverseDFTSolver<3, 3, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<3, 4, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<3, 5, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<3, 6, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<4, 4, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<4, 5, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<4, 6, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<4, 7, dftfe::utils::MemorySpace::DEVICE>;

template class BFGSInverseDFTSolver<5, 5, dftfe::utils::MemorySpace::DEVICE>;

template class BFGSInverseDFTSolver<5, 6, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<5, 7, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<5, 8, dftfe::utils::MemorySpace::DEVICE>;

template class BFGSInverseDFTSolver<6, 6, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<6, 7, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<6, 8, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<6, 9, dftfe::utils::MemorySpace::DEVICE>;

template class BFGSInverseDFTSolver<7, 7, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<7, 8, dftfe::utils::MemorySpace::DEVICE>;
template class BFGSInverseDFTSolver<7, 9, dftfe::utils::MemorySpace::DEVICE>;
#endif

} // namespace invDFT