non-hermitian davidson

pyblock2/driver/core.py

@@ -176,7 +176,7 @@ def dmrg(
  thrds = [1e-6] * 4 + [1e-7] * 1
  else:
  thrds = [1e-5] * 4 + [5e-6] * 1
- if dav_type == "ExactNonHermitian":
+ if dav_type is not None and "NonHermitian" in dav_type:
  bra = ket.deep_copy("BRA")
  else:
  bra = ket

src/core/matrix_functions.hpp

@@ -248,7 +248,8 @@ enum struct DavidsonTypes : uint16_t {
  NoPrecond = 64,
  NonHermitian = 128,
  Exact = 256,
- ExactNonHermitian = 128 | 256
+ ExactNonHermitian = 128 | 256,
+ NonHermitianDavidsonPrecond = 128 | 32
 };
 
 inline bool operator&(DavidsonTypes a, DavidsonTypes b) {

src/main.cpp

@@ -788,31 +788,61 @@ int main(int argc, char *argv[]) {
 
  if (single_prec && su2 && use_complex && !sgf) {
  cout << "SPIN-ADAPTED SINGLE-PREC COMPLEX" << endl;
+#ifdef _USE_COMPLEX
+#ifdef _USE_SINGLE_PREC
  run<SU2, complex<float>>(params);
+#endif
+#endif
  } else if (single_prec && !su2 && use_complex && !sgf) {
  cout << "NON-SPIN-ADAPTED SINGLE-PREC COMPLEX" << endl;
+#ifdef _USE_COMPLEX
+#ifdef _USE_SINGLE_PREC
  run<SZ, complex<float>>(params);
+#endif
+#endif
  } else if (single_prec && !su2 && use_complex && sgf) {
  cout << "GENERAL-SPIN SINGLE-PREC COMPLEX" << endl;
+#ifdef _USE_SG
+#ifdef _USE_COMPLEX
+#ifdef _USE_SINGLE_PREC
  run<SGF, complex<float>>(params);
+#endif
+#endif
+#endif
  } else if (su2 && use_complex && !sgf) {
  cout << "SPIN-ADAPTED COMPLEX" << endl;
+#ifdef _USE_COMPLEX
  run<SU2, complex<double>>(params);
+#endif
  } else if (use_complex && !sgf) {
  cout << "NON-SPIN-ADAPTED COMPLEX" << endl;
+#ifdef _USE_COMPLEX
  run<SZ, complex<double>>(params);
+#endif
  } else if (use_complex && sgf) {
  cout << "GENERAL-SPIN COMPLEX" << endl;
+#ifdef _USE_SG
+#ifdef _USE_COMPLEX
  run<SGF, complex<double>>(params);
+#endif
+#endif
  } else if (single_prec && su2 && !sgf) {
  cout << "SPIN-ADAPTED SINGLE-PREC" << endl;
+#ifdef _USE_SINGLE_PREC
  run<SU2, float>(params);
+#endif
  } else if (single_prec && !su2 && !sgf) {
  cout << "NON-SPIN-ADAPTED SINGLE-PREC" << endl;
+#ifdef _USE_SINGLE_PREC
  run<SZ, float>(params);
+#endif
  } else if (single_prec && sgf) {
  cout << "GENERAL-SPIN SINGLE-PREC" << endl;
+#ifdef _USE_SG
+#ifdef _USE_SINGLE_PREC
  run<SGF, float>(params);
+#endif
+#endif
  } else if (su2 && !sgf) {
  cout << "SPIN-ADAPTED" << endl;
  run<SU2, double>(params);
@@ -821,7 +851,9 @@ int main(int argc, char *argv[]) {
  run<SZ, double>(params);
  } else if (sgf) {
  cout << "GENERAL-SPIN" << endl;
+#ifdef _USE_SG
  run<SGF, double>(params);
+#endif
  }
 
  auto finish = chrono::system_clock::now();

src/pybind/pybind_core.hpp

@@ -1521,6 +1521,8 @@ template <typename S = void> void bind_types(py::module &m) {
  .value("NonHermitian", DavidsonTypes::NonHermitian)
  .value("Exact", DavidsonTypes::Exact)
  .value("ExactNonHermitian", DavidsonTypes::ExactNonHermitian)
+ .value("NonHermitianDavidsonPrecond",
+ DavidsonTypes::NonHermitianDavidsonPrecond)
  .def(py::self & py::self)
  .def(py::self | py::self);
 

unit_test/debug_test_dmrg_sa_n2_sto3g_nh.cpp

@@ -92,8 +92,8 @@ void TestDMRGN2STO3GSA<FL>::test_dmrg(
  // DMRG
  shared_ptr<DMRG<S, FL, FL>> dmrg =
  make_shared<DMRG<S, FL, FL>>(me, bdims, noises);
- dmrg->davidson_type = DavidsonTypes::ExactNonHermitian;
- dmrg->iprint = 3;
+ dmrg->davidson_type = DavidsonTypes::NonHermitianDavidsonPrecond;
+ dmrg->iprint = 2;
  dmrg->noise_type = NoiseTypes::ReducedPerturbativeCollected;
  dmrg->trunc_type = dmrg->trunc_type | TruncationTypes::RealDensityMatrix;
  dmrg->cutoff = 1E-20;
@@ -126,6 +126,7 @@ TYPED_TEST_CASE(TestDMRGN2STO3GSA, TestFL);
 
 TYPED_TEST(TestDMRGN2STO3GSA, TestSU2) {
  using FL = TypeParam;
+ using FLL = typename GMatrix<FL>::FL;
 
  shared_ptr<FCIDUMP<FL>> fcidump = make_shared<FCIDUMP<FL>>();
  PGTypes pg = PGTypes::D2H;
@@ -170,6 +171,7 @@ TYPED_TEST(TestDMRGN2STO3GSA, TestSU2) {
 
 TYPED_TEST(TestDMRGN2STO3GSA, TestSZ) {
  using FL = TypeParam;
+ using FLL = typename GMatrix<FL>::FL;
 
  shared_ptr<FCIDUMP<FL>> fcidump = make_shared<FCIDUMP<FL>>();
  PGTypes pg = PGTypes::D2H;

unit_test/test_matrix.cpp

@@ -1085,7 +1085,7 @@ TYPED_TEST(TestMatrix, TestEigRealNonSymmetric) {
  }
 }
 
-TYPED_TEST(TestMatrix, TestDavidsonRealNonSymmetric) {
+TYPED_TEST(TestMatrix, TestDavidsonRealNonSymmetricExact) {
  using FL = TypeParam;
  const int sz = is_same<FL, double>::value ? 200 : 75;
  const FL conv = is_same<FL, double>::value ? 1E-8 : 1E-7;
@@ -1167,6 +1167,93 @@ TYPED_TEST(TestMatrix, TestDavidsonRealNonSymmetric) {
  }
 }
 
+TYPED_TEST(TestMatrix, TestDavidsonRealNonSymmetric) {
+ using FL = TypeParam;
+ const int sz = is_same<FL, double>::value ? 120 : 50;
+ const FL conv = is_same<FL, double>::value ? 1E-14 : 1E-7;
+ const FL thrd = is_same<FL, double>::value ? 1E-3 : 1E+2;
+ const FL thrd2 = is_same<FL, double>::value ? 5E-1 : 1E+2;
+ const FL thrd3 = is_same<FL, double>::value ? 1E+0 : 1E+3;
+ using MatMul = typename TestMatrix<FL>::MatMul;
+ for (int i = 0; i < this->n_tests; i++) {
+ MKL_INT n = Random::rand_int(1, sz);
+ MKL_INT k = min(n, (MKL_INT)Random::rand_int(1, 10));
+ int ndav = 0;
+ GMatrix<FL> a(dalloc_<FL>()->allocate(n * n), n, n);
+ GMatrix<FL> ap(dalloc_<FL>()->allocate(n * n), n, n);
+ GMatrix<FL> ag(dalloc_<FL>()->allocate(n * n), n, n);
+ GMatrix<FL> ax(dalloc_<FL>()->allocate(n * n), n, n);
+ GDiagonalMatrix<FL> aa(dalloc_<FL>()->allocate(n), n);
+ GDiagonalMatrix<FL> ww(dalloc_<FL>()->allocate(n), n);
+ GDiagonalMatrix<FL> wi(dalloc_<FL>()->allocate(n), n);
+ vector<GMatrix<FL>> bs(k * 2, GMatrix<FL>(nullptr, n, 1));
+ Random::fill<FL>(a.data, a.size());
+ GMatrixFunctions<FL>::qr(a, ap, ag);
+ GMatrixFunctions<FL>::iadd(a, ap, (FL)0.95, false, (FL)0.05);
+ Random::fill<FL>(ww.data, ww.size());
+ GMatrixFunctions<FL>::copy(ag, a);
+ GMatrixFunctions<FL>::inverse(a);
+ for (MKL_INT ki = 0; ki < n; ki++)
+ for (MKL_INT kj = 0; kj < n; kj++)
+ ap(ki, kj) = ww(ki, kj);
+ ww.clear();
+ GMatrixFunctions<FL>::multiply(a, false, ap, false, ax, 1.0, 0.0);
+ GMatrixFunctions<FL>::multiply(ax, false, ag, false, a, 1.0, 0.0);
+ for (MKL_INT ki = 0; ki < n; ki++)
+ aa(ki, ki) = a(ki, ki);
+ for (int i = 0; i < k * 2; i++) {
+ bs[i].allocate();
+ bs[i].clear();
+ bs[i].data[i] = 1;
+ }
+ MatMul mop(a);
+ vector<FL> vw = IterativeMatrixFunctions<FL>::davidson(
+ mop, aa, bs, 0,
+ DavidsonTypes::NonHermitian | DavidsonTypes::DavidsonPrecond, ndav,
+ false, (shared_ptr<ParallelCommunicator<SZ>>)nullptr, conv,
+ n * k * 50, n * k * 40, k * 2, min(max((MKL_INT)5, k * 5 + 10), n));
+ ASSERT_EQ((int)vw.size(), k);
+ GDiagonalMatrix<FL> w(&vw[0], k);
+ GMatrixFunctions<FL>::eig(a, ww, wi, ap);
+ ASSERT_TRUE(GMatrixFunctions<FL>::all_close(
+ wi, GIdentityMatrix<FL>(n, (FL)0.0), thrd2, thrd2));
+ vector<int> idxs(n);
+ for (int i = 0; i < n; i++)
+ idxs[i] = i;
+ sort(idxs.begin(), idxs.begin() + n,
+ [&ww](int i, int j) { return ww.data[i] < ww.data[j]; });
+ GDiagonalMatrix<FL> w2(wi.data, k);
+ for (int i = 0; i < k; i++)
+ w2.data[i] = ww.data[idxs[i]];
+ ASSERT_TRUE(GMatrixFunctions<FL>::all_close(w, w2, thrd, thrd));
+ for (int i = 0; i < k; i++)
+ ASSERT_TRUE(GMatrixFunctions<FL>::all_close(
+ bs[i], GMatrix<FL>(a.data + a.n * idxs[i], a.n, 1),
+ thrd2, thrd2) ||
+ GMatrixFunctions<FL>::all_close(
+ bs[i], GMatrix<FL>(a.data + a.n * idxs[i], a.n, 1),
+ thrd2, thrd2, -1.0));
+ for (int i = 0; i < k; i++) {
+ ASSERT_TRUE(
+ GMatrixFunctions<FL>::all_close(
+ bs[i + k], GMatrix<FL>(ap.data + a.n * idxs[i], a.n, 1),
+ thrd3, thrd3) ||
+ GMatrixFunctions<FL>::all_close(
+ bs[i + k], GMatrix<FL>(ap.data + a.n * idxs[i], a.n, 1),
+ thrd3, thrd3, -1.0));
+ }
+ for (int i = k * 2 - 1; i >= 0; i--)
+ bs[i].deallocate();
+ wi.deallocate();
+ ww.deallocate();
+ aa.deallocate();
+ ax.deallocate();
+ ag.deallocate();
+ ap.deallocate();
+ a.deallocate();
+ }
+}
+
 TYPED_TEST(TestMatrix, TestSVD) {
  using FL = TypeParam;
  const int sz = is_same<FL, double>::value ? 200 : 75;