Program Listing for File sparseSolver.hpp

↰ Return to documentation for file (utils/sparseSolver.hpp)

#pragma once
#include "functions.hpp"
#include "qmatrix.hpp"
#include "timer.hpp"
#include <algorithm>
#include <cstddef>
#define EIGEN_USE_MKL_ALL
#include <eigen3/Eigen/Sparse>
#include <iostream>
#include <vector>
// Spectra - header-only sparse eigenvalue solver
#include <Spectra/GenEigsRealShiftSolver.h>
#include <Spectra/MatOp/SparseGenRealShiftSolve.h>

class exactSolver {
  qmatrix<double> qmat;
  size_t          column;

public:
  explicit exactSolver(size_t n = 0) : column(n) { qmat.resize(n, n, 0); }

  void set(size_t i, size_t j, double val) {
    qmat(i, j) = qmat(i, j) + val;
  }

  std::vector<double> solve() {
    auto                diagSolver = qmat.nonsys_diag_real();
    auto                leftVec    = std::get<0>(diagSolver);
    auto                rightVec   = std::get<1>(diagSolver);
    auto                eig        = std::get<2>(diagSolver);
    std::vector<double> X(column, {0});
    double              trace = 0;
    size_t              idx   = minIndex(eig); // index of minimum absolute value
#pragma omp parallel for reduction(+ : trace)
    for (size_t i = 0; i < column; i++) {
      X[i] = rightVec(i, idx).real();
      trace += X[i];
    }
    std::cout << "EigenValue" << eig[idx] << " trace :" << trace << std::endl;
    std::transform(X.begin(), X.end(), X.begin(), [trace](double x) { return x / trace; });
    return X;
  }

  auto vectorDot(std::vector<double> &a) {
    auto tmMat = qmat.dot(qmatrix<double>(a, column, 1));
    return std::vector<double>(tmMat.begin(), tmMat.end());
  }
};

class sparseEigen {
  using datatype = double;
  std::vector<Eigen::Triplet<datatype>> coefficients;
  size_t                        column;
  Eigen::SparseMatrix<datatype> SparseA;

public:
  explicit sparseEigen(size_t n = 0)
      : column(n), SparseA(Eigen::SparseMatrix<datatype>(n, n)) {}

  void set(size_t i, size_t j, double value) {
    coefficients.emplace_back(i, j, value);
  }

  auto getMatrix() {
    SparseA.setFromTriplets(coefficients.begin(), coefficients.end());
    return SparseA;
  }

  auto solve() {
    SparseA.setFromTriplets(coefficients.begin(), coefficients.end());

    Spectra::SparseGenRealShiftSolve<double> op(SparseA);
    Spectra::GenEigsRealShiftSolver          eigs(op, 1, column / 2, 1e-6);

    eigs.init();
    eigs.compute(Spectra::SortRule::LargestMagn);

    std::vector<double> X(column, {0});
    if (eigs.info() == Spectra::CompInfo::Successful) {
      double           rhoTrace = 0.0;
      Eigen::VectorXcd evalues  = eigs.eigenvalues();
      auto             eigVec   = eigs.eigenvectors();

      for (size_t i = 0; i < column; i++) {
        X[i]     = eigVec(i, 0).real();
        rhoTrace = rhoTrace + X[i];
      }

      for (auto &aa : X) {
        aa = aa / rhoTrace;
      }
      std::cout << "Lowest Eigenvector : " << evalues << "Trace:  " << rhoTrace << std::endl;
    } else {
      throw std::runtime_error("Failed to evaluate Eigenvalues\n");
    }
    return X;
  }

  auto vectorDot(const std::vector<double> &X) {
    SparseA.setFromTriplets(coefficients.begin(), coefficients.end());
    std::vector<double> result(column, 0);

    for (int k = 0; k < SparseA.outerSize(); ++k) {
      for (Eigen::SparseMatrix<datatype>::InnerIterator it(SparseA, k); it; ++it) {
        result[it.row()] += it.value() * X[it.col()];
      }
    }
    return result;
  }
};