Commit e9949e08 authored by Nguyen, Thien's avatar Nguyen, Thien

Fleshing out QLanczos algorithm

Following descriptions in the paper to regularize/stabilize the matrices.
Signed-off-by: Nguyen, Thien's avatarThien Nguyen <nguyentm@ornl.gov>
parent b3c0e6b8
......@@ -685,7 +685,76 @@ double QLanczos::calcQlanczosEnergy(const std::vector<double>& normVec) const
j++;
}
// TODO
const auto matFieldSampling = [](const arma::mat& in_H, const arma::mat& in_S, const std::vector<size_t> in_indexVec) {
arma::mat H(in_indexVec.size(), in_indexVec.size(), arma::fill::zeros);
arma::mat S(in_indexVec.size(), in_indexVec.size(), arma::fill::zeros);
assert(in_S.n_cols == in_H.n_cols && in_S.n_rows == in_H.n_rows);
for (size_t i = 0; i < in_indexVec.size(); ++i)
{
for (size_t j = 0; j < in_indexVec.size(); ++j)
{
if (xacc::container::contains(in_indexVec, i) && xacc::container::contains(in_indexVec, j))
{
S(i, j) = in_S(in_indexVec[i], in_indexVec[j]);
H(i, j) = in_H(in_indexVec[i], in_indexVec[j]);
}
}
}
return std::make_pair(H, S);
};
// Regularize/stabilize H and S matrices with s and epsilon stabilization parameters.
// Returns the stabilized (H, S) matrices.
const auto regularizeMatrices = [&](double s, double eps) -> std::pair<arma::mat, arma::mat> {
std::vector<size_t> indexVec { 0 };
size_t ii = 0;
size_t jj = 0;
while (ii < H.n_rows && jj < (H.n_rows - 1))
{
for(jj = ii + 1; jj < H.n_rows; ++jj)
{
if(S(ii, jj) < s)
{
indexVec.emplace_back(jj);
break;
}
}
ii = indexVec.back();
}
indexVec.emplace_back(H.n_rows);
auto [Hnew, Snew] = matFieldSampling(H, S, indexVec);
arma::vec sigma;
arma::mat V;
arma::eig_sym(sigma, V, Snew);
std::vector<size_t> indexVecEig;
for (size_t i = 0; i < sigma.n_elem; ++i)
{
if (sigma[i] > eps)
{
indexVecEig.emplace_back(i);
}
}
auto Snew2 = V.t() * Snew * V;
auto Hnew2 = V.t() * Hnew * V ;
auto [Hnew3, Snew3] = matFieldSampling(Hnew2, Snew2, indexVecEig);
return std::make_pair(Hnew3, Snew3);
};
// Stabilization parameters:
// TODO: enables as heterogenous params
const double s = 0.75;
const double epsilon = 1e-12;
auto [Hreg, Sreg] = regularizeMatrices(s, epsilon);
// TODO: solving the *generalized* eigenvalues: Hx = ESx
// to compute the Lanczos energy E
// TODO: returns the smallest eigenvalue from the previous step.
return 0.0;
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment