Green's Functions

include/macis/gf/bandlan.hpp

@@ -30,36 +30,91 @@
 #include <utility>
 
 #include "macis/gf/inn_prods.hpp"
+#include "macis/solvers/davidson.hpp"
 
 typedef std::numeric_limits<double> dbl;
 
 namespace macis {
 
-extern "C" {
-extern int dgeqrf_(int *, int *, double *, int *, double *, double *, int *,
-                   int *);
-extern int dorgqr_(int *, int *, int *, double *, int *, double *, double *,
-                   int *, int *);
-extern int dsbtrd_(char *, char *, int *, int *, double *, int *, double *,
-                   double *, double *, int *, double *, int *);
-extern int dsytrd_(char *, int *, double *, int *, double *, double *, double *,
-                   double *, int *, int *);
-extern int dorgtr_(char *, int *, double *, int *, double *, double *, int *,
-                   int *);
-extern int dsteqr_(char *, int *, double *, double *, double *, int *, double *,
-                   int *);
-}
-
+/**
+ * @ brief Wrapper for QR decomposition in LAPACK, basically
+ *         calling dgeqrf and dorgqr to evaluate the R and
+ *         Q matrices, which are returned. Here, the input
+ *         matrix has more rows than columns.
+ *
+ * @param[inout] std::vector<std::vector<double> > &Q: On input,
+ *               matrix for which to evaluate the QR decomposition.
+ *               On output, Q-matrix.
+ * @param[out] std::vector<std::vector<double> > &R: On output, R
+ *             matrix in the QR decomposition.
+ *
+ * @returns bool: Error code from LAPACK routines.
+ *
+ * @author Carlos Mejuto-Zaera
+ * @date 25/04/2022
+ */
 bool QRdecomp(std::vector<std::vector<double> > &Q,
               std::vector<std::vector<double> > &R);
 
+/**
+ * @ brief Wrapper for QR decomposition in LAPACK, basically
+ *         calling dgeqrf and dorgqr to evaluate the R and
+ *         Q matrices, which are returned. Here, the input
+ *         matrix has more columns than rows.
+ *
+ * @param[inout] std::vector<std::vector<double> > &Q: On input,
+ *               matrix for which to evaluate the QR decomposition.
+ *               On output, Q-matrix.
+ * @param[out] std::vector<std::vector<double> > &R: On output, R
+ *             matrix in the QR decomposition.
+ *
+ * @returns bool: Error code from LAPACK routines.
+ *
+ * @author Carlos Mejuto-Zaera
+ * @date 25/04/2022
+ */
 bool QRdecomp_tr(std::vector<std::vector<double> > &Q,
                  std::vector<std::vector<double> > &R);
 
+/**
+ * @brief Wrapper to LAPACK routine to evaluate the eigenvectors
+ *        and eigenvalues of the symmetric matrix mat.
+ *
+ * @param[inout] std::vector<std::vector<double> > &mat: Matrix for
+ *               which to compute the eigenvalues/vectors. Erased
+ *               during computation.
+ * @param[out] std::vector<double> &eigvals: Eigenvalues, sorted from smallest
+ *             to largest.
+ * @param[out] std::vector<std::vector<double> > &eigvecs: Eigenvectors,
+ *             stored as row vectors.
+ *
+ * @returns bool: Error code from LAPACK.
+ *
+ * @author Carlos Mejuto Zaera
+ * @date 25/04/2022
+ */
 bool GetEigsys(std::vector<std::vector<double> > &mat,
                std::vector<double> &eigvals,
                std::vector<std::vector<double> > &eigvecs);
 
+/**
+ * @brief Wrapper to LAPACK routine to evaluate the eigenvectors
+ *        and eigenvalues of the symmetric band matrix mat.
+ *
+ * @param[inout] std::vector<std::vector<double> > &mat: Matrix for
+ *               which to compute the eigenvalues/vectors. Erased
+ *               during computation.
+ * @param[in] int nSupDiag: Nr. of bands.
+ * @param[out] std::vector<double> &eigvals: Eigenvalues, sorted from smallest
+ *             to largest.
+ * @param[out] std::vector<std::vector<double> > &eigvecs: Eigenvectors,
+ *             stored as row vectors.
+ *
+ * @returns bool: Error code from LAPACK.
+ *
+ * @author Carlos Mejuto Zaera
+ * @date 25/04/2022
+ */
 bool GetEigsysBand(std::vector<std::vector<double> > &mat, int nSupDiag,
                    std::vector<double> &eigvals,
                    std::vector<std::vector<double> > &eigvecs);
@@ -85,32 +140,35 @@ bool GetEigsysBand(std::vector<std::vector<double> > &mat, int nSupDiag,
  * @author Carlos Mejuto Zaera
  * @date 25/04/2022
  */
-template <class Cont>
-void MyBandLan(
-    const sparsexx::dist_sparse_matrix<sparsexx::csr_matrix<double, int32_t> >
-        &H,
-    std::vector<std::vector<Cont> > &qs, std::vector<std::vector<Cont> > &bandH,
-    int &nLanIts, double thres = 1.E-6, bool print = false) {
+template <class Cont, class Functor>
+void MyBandLan(const Functor &H,
+               // std::vector<std::vector<Cont> > &qs,
+               // std::vector<std::vector<Cont> > &bandH,
+               std::vector<Cont> &qs, std::vector<std::vector<Cont> > &bandH,
+               int &nLanIts, int nbands, int N, double thres = 1.E-6,
+               bool print = false) {
   // BAND LANCZOS ROUTINE. TAKES AS INPUT THE HAMILTONIAN H, INITIAL VECTORS qs
   // AND RETURNS THE BAND HAMILTONIAN bandH. IT PERFORMS nLanIts ITERATIONS,
   // STOPPING IF THE NORM OF ANY NEW KRYLOV VECTOR IS BELOW thres. IF LANCZOS IS
   // STOPPED PREMATURELY , nLanIts IS OVERWRITTEN WITH THE ACTUAL NUMBER OF
-  //ITERATIONS! THE qs VECTOR IS ERASED AT THE END OF THE CALCULATION
+  // ITERATIONS! THE qs VECTOR IS ERASED AT THE END OF THE CALCULATION
   bandH.clear();
   bandH.resize(nLanIts, std::vector<Cont>(nLanIts, 0.));
 
-  int nbands = qs.size();
-  auto spmv_info = sparsexx::spblas::generate_spmv_comm_info(H);
   // MAKE SPACE FOR 2 * nbands VECTORS
-  qs.resize(2 * nbands, std::vector<Cont>(qs[0].size(), 0.));
-  std::vector<Cont> temp(qs[0].size(), 0.);
+  // qs.resize(2 * nbands, std::vector<Cont>(qs[0].size(), 0.));
+  qs.resize(2 * nbands * N);
+  // std::vector<Cont> temp(qs[0].size(), 0.);
+  std::vector<Cont> temp(N, 0.);
   if(print) {
     for(int i = 0; i < nbands; i++) {
       std::ofstream ofile("lanvec_" + std::to_string(i + 1) + ".dat",
                           std::ios::out);
       ofile.precision(dbl::max_digits10);
-      for(size_t el = 0; el < qs[i].size(); el++)
-        ofile << std::scientific << qs[i][el] << std::endl;
+      // for(size_t el = 0; el < qs[i].size(); el++)
+      for(size_t el = 0; el < N; el++)
+        // ofile << std::scientific << qs[i][el] << std::endl;
+        ofile << std::scientific << qs[el + i * N] << std::endl;
       ofile.close();
     }
   }
@@ -125,8 +183,10 @@ void MyBandLan(
   for(int it = 1; it <= nLanIts; it++) {
     int band_indx_i =
         true_indx[it];  // TO WHAT ELEMENT OF THE VECTOR SET DO WE APPLY THIS
-    sparsexx::spblas::pgespmv(1., H, qs[band_indx_i].data(), 0., temp.data(),
-                              spmv_info);
+    // H.operator_action( 1, 1., qs[band_indx_i].data(), temp.size(), 0.,
+    // temp.data(), temp.size() );
+    H.operator_action(1, 1., qs.data() + band_indx_i * N, N, 0., temp.data(),
+                      N);
     if(print) {
       std::ofstream ofile("Htimes_lanvec_" + std::to_string(it) + ".dat",
                           std::ios::out);
@@ -140,15 +200,19 @@ void MyBandLan(
       int band_indx_j = true_indx[jt];
 #pragma omp parallel for
       for(size_t coeff = 0; coeff < temp.size(); coeff++)
-        temp[coeff] -= bandH[it - 1][jt - 1] * qs[band_indx_j][coeff];
+        // temp[coeff] -= bandH[it - 1][jt - 1] * qs[band_indx_j][coeff];
+        temp[coeff] -= bandH[it - 1][jt - 1] * qs[N * band_indx_j + coeff];
     }
     for(int jt = it; jt <= std::min(it + nbands - 1, nLanIts); jt++) {
       int band_indx_j = true_indx[jt];
-      bandH[it - 1][jt - 1] = MyInnProd(temp, qs[band_indx_j]);
+      // bandH[it - 1][jt - 1] = MyInnProd(temp, qs[band_indx_j]);
+      bandH[it - 1][jt - 1] =
+          blas::dot(N, temp.data(), 1, qs.data() + band_indx_j * N, 1);
       bandH[jt - 1][it - 1] = bandH[it - 1][jt - 1];
 #pragma omp parallel for
       for(size_t coeff = 0; coeff < temp.size(); coeff++)
-        temp[coeff] -= bandH[it - 1][jt - 1] * qs[band_indx_j][coeff];
+        // temp[coeff] -= bandH[it - 1][jt - 1] * qs[band_indx_j][coeff];
+        temp[coeff] -= bandH[it - 1][jt - 1] * qs[N * band_indx_j + coeff];
     }
     if(it + nbands <= nLanIts) {
       bandH[it - 1][it + nbands - 1] =
@@ -166,20 +230,24 @@ void MyBandLan(
         break;
 #pragma omp parallel for
         for(size_t coeff = 0; coeff < temp.size(); coeff++)
-          qs[true_indx[it + nbands]][coeff] = 0.;
+          // qs[true_indx[it + nbands]][coeff] = 0.;
+          qs[true_indx[it + nbands] * N + coeff] = 0.;
         std::cout << "FOUND A ZERO VECTOR AT POSITION " << next_indx
                   << std::endl;
       } else {
 #pragma omp parallel for
         for(size_t coeff = 0; coeff < temp.size(); coeff++)
-          qs[true_indx[it + nbands]][coeff] =
+          // qs[true_indx[it + nbands]][coeff] =
+          qs[true_indx[it + nbands] * N + coeff] =
               temp[coeff] / bandH[it - 1][it + nbands - 1];
         if(print) {
           std::ofstream ofile("lanvec_" + std::to_string(it + nbands) + ".dat",
                               std::ios::out);
           ofile.precision(dbl::max_digits10);
-          for(size_t el = 0; el < qs[true_indx[it + nbands]].size(); el++)
-            ofile << std::scientific << qs[true_indx[it + nbands]][el]
+          // for(size_t el = 0; el < qs[true_indx[it + nbands]].size(); el++)
+          for(size_t el = 0; el < N; el++)
+            // ofile << std::scientific << qs[true_indx[it + nbands]][el]
+            ofile << std::scientific << qs[true_indx[it + nbands] * N + el]
                   << std::endl;
           ofile.close();
         }

include/macis/gf/eigsolver.hpp

@@ -17,20 +17,13 @@
 #pragma once
 #include <Eigen/Core>
 #include <Eigen/Sparse>
+#include <lapack.hh>
 
 namespace macis {
 
-using namespace Eigen;
-
-typedef MatrixXd eigMatD;
-typedef SparseMatrix<double, RowMajor> SpMatD;
-
-extern "C" {
-extern int dsteqr_(char *, int *, double *, double *, double *, int *, double *,
-                   int *);
-extern int dsyev_(char *, char *, int *, double *, int *, double *, double *,
-                  int *, int *);
-}
+typedef Eigen::VectorXd VectorXd;
+typedef Eigen::MatrixXd eigMatD;
+typedef Eigen::SparseMatrix<double, Eigen::RowMajor> SpMatD;
 
 /**
  * @brief Computes the eigenvalues and eigenvectors of a tridiagonal, symmetric
@@ -45,7 +38,7 @@ extern int dsyev_(char *, char *, int *, double *, int *, double *, double *,
  * @date 05/04/2021
  */
 void Hste_v(const std::vector<double> &alphas, const std::vector<double> &betas,
-            VectorXd &eigvals, eigMatD &eigvecs);
+            Eigen::VectorXd &eigvals, eigMatD &eigvecs);
 
 /**
  * @brief Computes the eigenvalues and eigenvectors of a tridiagonal, symmetric
@@ -59,6 +52,6 @@ void Hste_v(const std::vector<double> &alphas, const std::vector<double> &betas,
  * @author Carlos Mejuto Zaera
  * @date 05/04/2021
  */
-void Hsyev(const eigMatD &H, VectorXd &eigvals, eigMatD &eigvecs);
+void Hsyev(const eigMatD &H, Eigen::VectorXd &eigvals, eigMatD &eigvecs);
 
 }  // namespace macis

include/macis/gf/gf.hpp

@@ -37,7 +37,7 @@
 
 namespace macis {
 
-using namespace Eigen;
+typedef Eigen::VectorXd VectorXd;
 
 typedef std::numeric_limits<double> dbl;
 typedef std::chrono::high_resolution_clock Clock;
@@ -165,21 +165,6 @@ void GF_Diag(const VectorXd &state_0, const MatOp &H,
   }
 }
 
-/**
- * @brief Class for comparing two determinants, to use in sorting routines.
- *
- * @author Carlos Mejuto Zaera
- * @date 28/01/2022
- */
-template <size_t nbits>
-class BitSetComparator {
- public:
-  bool operator()(const std::bitset<nbits> &c1,
-                  const std::bitset<nbits> &c2) const {
-    return bitset_less(c1, c2);
-  }
-};
-
 /**
  * @brief Routine to build basis to compute the basis for Green's function
  *        calculation on the state |wfn> considering orbital orb. Templated
@@ -244,10 +229,10 @@ void get_GF_basis_AS_1El(int orb, bool sp_up, bool is_part, const VectorXd &wfn,
   std::bitset<nbits> uni_string;
   std::vector<std::bitset<nbits>> founddets;
   founddets.reserve(trunc_size);
-  std::map<std::bitset<nbits>, size_t, BitSetComparator<nbits>> founddet_pos,
-      basedet_pos;
+  std::map<std::bitset<nbits>, size_t, bitset_less_comparator<nbits>>
+      founddet_pos, basedet_pos;
   typename std::map<std::bitset<nbits>, size_t,
-                    BitSetComparator<nbits>>::iterator it;
+                    bitset_less_comparator<nbits>>::iterator it;
   // ACTIVE SPACE
   std::vector<uint32_t> as_orbs;
   for(size_t i = 0; i < occs.size(); i++) {
@@ -401,7 +386,8 @@ std::vector<std::vector<double>> BuildWfn4Lanczos(
     const std::vector<std::bitset<nbits>> &GF_dets, bool is_part,
     std::vector<int> &todelete, double zero_thresh = 1.E-7) {
   // INITIALIZE THE DICTIONARY OF BASE DETERMINANTS
-  std::map<std::bitset<nbits>, size_t, BitSetComparator<nbits>> base_dets_pos;
+  std::map<std::bitset<nbits>, size_t, bitset_less_comparator<nbits>>
+      base_dets_pos;
   for(size_t iii = 0; iii < base_dets.size(); iii++)
     base_dets_pos[base_dets[iii]] = iii;
 
@@ -428,7 +414,7 @@ std::vector<std::vector<double>> BuildWfn4Lanczos(
         std::bitset<nbits> temp(GF_dets[ndet]);
         temp.flip(sporb);
         typename std::map<std::bitset<nbits>, size_t,
-                          BitSetComparator<nbits>>::const_iterator it =
+                          bitset_less_comparator<nbits>>::const_iterator it =
             base_dets_pos.find(temp);
         if(it != base_dets_pos.end()) {
           // IT DOES COME INDEED FROM A DETERMINANT IN THE ORIGINAL GROUND STATE

include/macis/gf/inn_prods.hpp

@@ -32,14 +32,7 @@ namespace macis {
  */
 inline double MyInnProd(const std::vector<double>& vecR,
                         const std::vector<double>& vecL) {
-  // SIMPLE INNER PRODUCT ROUTINE
-  double res = 0.;
-#pragma omp declare reduction(Vsum:double                   \
-                              : omp_out = omp_out + omp_in) \
-    initializer(omp_priv = 0.)
-#pragma omp parallel for reduction(Vsum : res)
-  for(size_t i = 0; i < vecR.size(); i++) res += vecR[i] * vecL[i];
-  return res;
+  return blas::dot(vecR.size(), vecR.data(), 1, vecL.data(), 1);
 }
 
 /**
@@ -56,14 +49,7 @@ inline double MyInnProd(const std::vector<double>& vecR,
 inline std::complex<double> MyInnProd(
     const std::vector<std::complex<double> >& vecR,
     const std::vector<std::complex<double> >& vecL) {
-  // SIMPLE INNER PRODUCT ROUTINE
-  std::complex<double> res(0., 0.);
-#pragma omp declare reduction \
-        (Vsum:std::complex<double>:omp_out=omp_out+omp_in)\
-        initializer(omp_priv=std::complex<double>(0.,0.))
-#pragma omp parallel for reduction(Vsum : res)
-  for(size_t i = 0; i < vecR.size(); i++) res += conj(vecR[i]) * vecL[i];
-  return res;
+  return blas::dot(vecR.size(), vecR.data(), 1, vecL.data(), 1);
 }
 
 /**
@@ -80,15 +66,7 @@ inline std::complex<double> MyInnProd(
 inline std::complex<double> MyInnProd(
     const std::vector<std::complex<double> >& vecR,
     const std::vector<double>& vecL) {
-  // SIMPLE INNER PRODUCT ROUTINE
-  std::complex<double> res(0., 0.);
-#pragma omp declare reduction \
-        (Vsum:std::complex<double>:omp_out=omp_out+omp_in)\
-        initializer(omp_priv=std::complex<double>(0.,0.))
-#pragma omp parallel for reduction(Vsum : res)
-  for(size_t i = 0; i < vecR.size(); i++)
-    res += conj(vecR[i]) * std::complex<double>(vecL[i], 0.);
-  return res;
+  return blas::dot(vecR.size(), vecR.data(), 1, vecL.data(), 1);
 }
 
 /**
@@ -103,8 +81,7 @@ inline std::complex<double> MyInnProd(
  */
 inline std::complex<double> MyInnProd(const Eigen::VectorXcd& vecR,
                                       const Eigen::VectorXcd& vecL) {
-  // SIMPLE INNER PRODUCT ROUTINE
-  return vecR.dot(vecL);
+  return blas::dot(vecR.size(), vecR.data(), 1, vecL.data(), 1);
 }
 
 /**
@@ -119,7 +96,6 @@ inline std::complex<double> MyInnProd(const Eigen::VectorXcd& vecR,
  */
 inline double MyInnProd(const Eigen::VectorXd& vecR,
                         const Eigen::VectorXd& vecL) {
-  // SIMPLE INNER PRODUCT ROUTINE
-  return vecR.dot(vecL);
+  return blas::dot(vecR.size(), vecR.data(), 1, vecL.data(), 1);
 }
 }  // namespace macis