Implement wrapper of sbd eigensolver

.github/workflows/test_latest_versions.yml

@@ -71,7 +71,12 @@ jobs:
         if: runner.os == 'Linux'
         run: |
           sudo apt-get update
-          sudo apt-get install -y g++ clang cmake ninja-build
+          sudo apt-get install -y g++ clang cmake ninja-build \
+            libopenblas-dev liblapack-dev libopenmpi-dev openmpi-bin libomp-dev
+      - name: Install dependencies (macOS)
+        if: runner.os == 'macOS'
+        run: |
+          brew install open-mpi libomp
       - name: Configure CMake
         run: >-
           cmake -S . -B build

.gitmodules

@@ -43,3 +43,6 @@
 [submodule "deps/bitset2"]
 	path = deps/bitset2
 	url = https://github.com/ClaasBontus/bitset2.git
+[submodule "deps/sbd"]
+	path = deps/sbd
+	url = https://github.com/r-ccs-cms/sbd.git

CMakeLists.txt

@@ -16,6 +16,7 @@ endif()
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
+# Enable warnings
 if (MSVC)
     add_compile_options(/W4)
 else()
@@ -28,6 +29,21 @@ include_directories(include)
 add_subdirectory(deps/doctest)
 add_subdirectory(deps/nanobench)
 
+# MPI and OpenMP
+find_package(MPI REQUIRED)
+find_package(OpenMP REQUIRED)
+
+# BLAS
+find_library(ACCELERATE_LIBRARY Accelerate)
+if (ACCELERATE_LIBRARY)
+    message(STATUS "Using Accelerate Framework")
+else()
+    set(BLA_VENDOR OpenBLAS)
+    find_package(BLAS REQUIRED)
+    include_directories(${BLAS_INCLUDE_DIRS})
+    link_libraries(${BLAS_LIBRARIES})
+endif()
+
 # boost::dynamic_bitset
 add_library(boost_dynamic_bitset INTERFACE)
 target_include_directories(boost_dynamic_bitset INTERFACE
@@ -50,18 +66,30 @@ target_include_directories(bitset2 INTERFACE
     ${CMAKE_CURRENT_SOURCE_DIR}/deps/bitset2
 )
 
+add_library(sbd INTERFACE)
+target_include_directories(sbd SYSTEM INTERFACE
+    ${CMAKE_CURRENT_SOURCE_DIR}/deps/sbd/include
+)
+
 add_executable(sqd_tests
     test/doctest_main.cpp
     test/test_postselection.cpp
     test/test_subsampling.cpp
     test/test_configuration_recovery.cpp
     test/test_fermion.cpp
+    test/test_sbd_eigensolver.cpp
+)
+target_include_directories(sqd_tests PRIVATE
+    deps/doctest/doctest
+    ${MPI_C_INCLUDE_DIRS}
 )
-target_include_directories(sqd_tests PRIVATE deps/doctest/doctest)
 target_link_libraries(sqd_tests
     PRIVATE
     boost_dynamic_bitset
     bitset2
+    sbd
+    MPI::MPI_C
+    OpenMP::OpenMP_CXX
 )
 
 include(${doctest_SOURCE_DIR}/scripts/cmake/doctest.cmake)
@@ -78,4 +106,5 @@ target_link_libraries(sqd_benchmarks
     nanobench
     boost_dynamic_bitset
     bitset2
+    sbd
 )

data/fcidump_mock.txt

@@ -0,0 +1,69 @@
+ &FCI NORB=  4,NELEC=4,MS2=0,
+  ORBSYM=1,1,1,1,
+  ISYM=1,
+ &END
+ 0.2840843372397457    1    1    1    1
+ 0.0002854816936369414    2    1    1    1
+ 0.04417807442382433    2    1    2    1
+ 0.2622091838554816    2    2    1    1
+ 0.0006468616587143983    2    2    2    1
+ 0.2730936605028639    2    2    2    2
+ -6.264064000199201e-05    3    1    1    1
+ 0.002101357535092042    3    1    2    1
+ 0.0002255741757854556    3    1    2    2
+ 0.03573411911251266    3    1    3    1
+ 0.0001350220017667965    3    2    1    1
+ 0.0004833657197940787    3    2    2    1
+ 6.60155191565777e-05    3    2    2    2
+ -0.0003287992979621698    3    2    3    1
+ 0.008161769399712743    3    2    3    2
+ 0.2615787261186168    3    3    1    1
+ -4.738450152304652e-05    3    3    2    1
+ 0.2582635280856485    3    3    2    2
+ 0.0006070623455779637    3    3    3    1
+ -0.002443459518692521    3    3    3    2
+ 0.2799043482114074    3    3    3    3
+ 0.004948020321549084    4    1    1    1
+ -0.0001423669308772119    4    1    2    1
+ 0.01355099400828998    4    1    2    2
+ -9.862707170814855e-05    4    1    3    1
+ 0.00690828349990277    4    1    3    2
+ -0.005561548693420539    4    1    3    3
+ 0.02473119582336451    4    1    4    1
+ 0.0002187162198814015    4    2    1    1
+ 0.0222393679678707    4    2    2    1
+ 0.0007535913042595605    4    2    2    2
+ 0.01456013639219515    4    2    3    1
+ 4.134481095727774e-05    4    2    3    2
+ 0.0007998106625613357    4    2    3    3
+ -0.000631142382644568    4    2    4    1
+ 0.0274480584179987    4    2    4    2
+ 0.0002458927780511462    4    3    1    1
+ 0.01592029588904504    4    3    2    1
+ 0.0001645551990137675    4    3    2    2
+ -0.02154892525410992    4    3    3    1
+ 0.0006066874646118994    4    3    3    2
+ -0.0005508824155951026    4    3    3    3
+ 7.170799160342295e-05    4    3    4    1
+ -0.002178587382308951    4    3    4    2
+ 0.03709458071687864    4    3    4    3
+ 0.2682959348286622    4    4    1    1
+ -0.0003600848098100874    4    4    2    1
+ 0.2731565325485292    4    4    2    2
+ 0.0002213930846584909    4    4    3    1
+ -0.0009824570080940911    4    4    3    2
+ 0.2773340631231709    4    4    3    3
+ 0.004390876956748526    4    4    4    1
+ 0.0004732853856294346    4    4    4    2
+ -0.000517460610178891    4    4    4    3
+ 0.3038445820400585    4    4    4    4
+ -12.52718671684592    1    1  0  0
+ 0.001483911550546467    2    1  0  0
+ -12.70422837256603    2    2  0  0
+ 0.006259774312687916    3    1  0  0
+ 0.04567009757527558    3    2  0  0
+ -12.83732437964672    3    3  0  0
+ -0.1188168544704107    4    1  0  0
+ -0.008552528363446082    4    2  0  0
+ -0.005989008051549467    4    3  0  0
+ -13.28315426027296    4    4  0  0

include/qiskit/addon/sqd/eigensolver/sbd.hpp

@@ -0,0 +1,115 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2025.
+//
+// This code is licensed under the Apache License, Version 2.0. You may
+// obtain a copy of this license in the LICENSE.txt file in the root directory
+// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+//
+// Any modifications or derivative works of this code must retain this
+// copyright notice, and modified files need to carry a notice indicating
+// that they have been altered from the originals.
+
+#ifndef QISKIT_ADDON_SQD_SUBSAMPLING_HPP_
+#define QISKIT_ADDON_SQD_SUBSAMPLING_HPP_
+
+/// Support for SBD eigensolver
+
+#include <string>
+#include <vector>
+
+#include "qiskit/addon/sqd/bitset_full.hpp"
+
+#if __has_include(<sbd/sbd.h>)
+
+#include <mpi.h>
+#include <sbd/sbd.h>
+
+namespace Qiskit
+{
+
+namespace addon
+{
+
+namespace sqd
+{
+
+namespace internal
+{
+
+}
+
+template <typename BitstringVectorType>
+struct SBDResult {
+    const double energy;
+    const std::vector<double> density;
+    BitstringVectorType carryover_bitstrings;
+    const std::vector<std::vector<double>> one_p_rdm, two_p_rdm;
+};
+
+template <typename BitstringVectorType>
+SBDResult<BitstringVectorType> sbd_solve(
+    const MPI_Comm &comm, const sbd::tpb::SBD &sbd_data, const sbd::FCIDump &fcidump,
+    const BitstringVectorType &alpha_determinants,
+    const BitstringVectorType &beta_determinants, const std::string &loadname = "",
+    const std::string &savename = ""
+)
+{
+    // Populate adet and bdet in storage format that sbd expects
+    std::vector<std::vector<size_t>> adet, bdet;
+    adet.reserve(alpha_determinants.size());
+    for (const auto &bitstring : alpha_determinants) {
+        adet.push_back(
+            sbd::from_string(
+                internal::bitset_to_string(bitstring), sbd_data.bit_length,
+                bitstring.size()
+            )
+        );
+    }
+    bdet.reserve(beta_determinants.size());
+    for (const auto &bitstring : beta_determinants) {
+        bdet.push_back(
+            sbd::from_string(
+                internal::bitset_to_string(bitstring), sbd_data.bit_length,
+                bitstring.size()
+            )
+        );
+    }
+
+    // Call sbd
+    double energy;
+    std::vector<double> density;
+    std::vector<std::vector<size_t>> carryover_sbd_bitstrings;
+    std::vector<std::vector<double>> one_p_rdm, two_p_rdm;
+    sbd::tpb::diag(
+        comm, sbd_data, fcidump, adet, bdet, loadname, savename, energy, density,
+        carryover_sbd_bitstrings, one_p_rdm, two_p_rdm
+    );
+
+    // Convert carryover bitstrings to desired bitstring format
+    BitstringVectorType carryover_bitstrings;
+    carryover_bitstrings.reserve(carryover_sbd_bitstrings.size());
+    for (const auto &sbd_bitstring : carryover_sbd_bitstrings) {
+        carryover_bitstrings.emplace_back(
+            sbd::makestring(
+                sbd_bitstring, sbd_data.bit_length, alpha_determinants[0].size()
+            )
+        );
+    }
+
+    // Construct and return result object
+    return {
+        energy, std::move(density), std::move(carryover_bitstrings),
+        std::move(one_p_rdm), std::move(two_p_rdm)
+    };
+}
+
+} // namespace sqd
+
+} // namespace addon
+
+} // namespace Qiskit
+
+#endif // __has_include(<sbd/sbd.h>)
+
+#endif // QISKIT_ADDON_SQD_SUBSAMPLING_HPP_

include/qiskit/addon/sqd/internal/bitset_common.hpp

@@ -18,6 +18,7 @@
 #include <array>
 #include <bitset>
 #include <cstddef>
+#include <string>
 
 namespace Qiskit
 {
@@ -55,6 +56,12 @@ std::array<HalfSize<std::bitset<N>>, 2> split_bitstring(const std::bitset<N> &bi
     return {right, left};
 }
 
+template <std::size_t N>
+std::string bitset_to_string(const std::bitset<N> &bitset)
+{
+    return bitset.to_string();
+}
+
 } // namespace internal
 
 } // namespace sqd

include/qiskit/addon/sqd/support/boost_dynamic_bitset.hpp

@@ -55,6 +55,14 @@ split_bitstring(const boost::dynamic_bitset<Block, Allocator> &bitset)
     return {right, left};
 }
 
+template <typename Block, typename Allocator>
+std::string bitset_to_string(const boost::dynamic_bitset<Block, Allocator> &bitstring)
+{
+    std::string str;
+    boost::to_string(bitstring, str);
+    return str;
+}
+
 } // namespace internal
 
 } // namespace sqd

test/doctest_main.cpp

@@ -10,9 +10,29 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
+#define DOCTEST_CONFIG_IMPLEMENT
 #include "doctest.h"
+#include <mpi.h>
 
+int main(int argc, char *argv[])
+{
+    // Initialize MPI
+    MPI_Init(&argc, &argv);
+
+    // Run the doctest tests without exiting
+    doctest::Context context(argc, argv);
+    context.setOption("no-run", false);
+    context.setOption("no-exitcode", true);
+    int result = context.run();
+
+    // Finalize MPI
+    MPI_Finalize();
+
+    // Now we can exit
+    return result;
+}
+
+// Special functions for when exceptions are disabled
 #if QKA_SQD_DISABLE_EXCEPTIONS
 #include <exception>
 #include <iostream>

test/test_sbd_eigensolver.cpp

@@ -0,0 +1,30 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2025.
+//
+// This code is licensed under the Apache License, Version 2.0. You may
+// obtain a copy of this license in the LICENSE.txt file in the root directory
+// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+//
+// Any modifications or derivative works of this code must retain this
+// copyright notice, and modified files need to carry a notice indicating
+// that they have been altered from the originals.
+
+#include "qiskit/addon/sqd/eigensolver/sbd.hpp"
+
+#include <bitset>
+#include <unordered_set>
+#include <vector>
+
+#include "bitset_compat.hpp"
+
+#include "doctest.h"
+
+TEST_CASE("Basic SBD test")
+{
+    auto sbd_data = sbd::tpb::generate_sbd_data(0, nullptr);
+    auto fcidump = sbd::LoadFCIDump("data/fcidump_mock.txt");
+    std::vector<std::bitset<4>> adets{3, 5, 6};
+    auto result =
+        Qiskit::addon::sqd::sbd_solve(MPI_COMM_WORLD, sbd_data, fcidump, adets, adets);
+}