Merge pull request #339 from SpM-lab/338-add-c-api-for-retrieving-default-sampling-for-spir_funcs

Add spir_basis_get_default_taus_ext

Author: shinaoka

include/sparseir/sparseir.h

@@ -598,6 +598,21 @@ int spir_basis_get_n_default_ws(const spir_basis *b, int *num_points);
  */
 int spir_basis_get_default_ws(const spir_basis *b, double *points);
 
+
+/***
+ * @brief Gets the default tau sampling points for ann IR basis.
+ *
+ * This function returns default tau sampling points for an IR basis object.
+ *
+ * @param b Pointer to the basis object
+ * @param n_points Number of requested sampling points.
+ * @param points Pre-allocated array to store the sampling points. The size of the array must be at least n_points.
+ * @param n_points_returned Number of sampling points returned.
+ * @return An integer status code:
+ *         - 0 (SPIR_COMPUTATION_SUCCESS) on success
+ */
+int spir_basis_get_default_taus_ext(const spir_basis *b, int n_points, double *points, int *n_points_returned);
+
 /**
  * @brief Gets the number of default Matsubara sampling points for an IR basis.
  *

src/cinterface.cpp

@@ -900,6 +900,56 @@ int spir_basis_get_default_taus(const spir_basis *b, double *points)
     }
 }
 
+int spir_basis_get_default_taus_ext(
+    const spir_basis *b, int n_points, double *points, int *n_points_returned)
+{
+    if (!b || !points) {
+        return SPIR_INVALID_ARGUMENT;
+    }
+
+    auto impl = get_impl_basis(b);
+    if (!impl) {
+        return SPIR_GET_IMPL_FAILED;
+    }
+
+    if (!is_ir_basis(b)) {
+        DEBUG_LOG("Error: The basis is not an IR basis");
+        return SPIR_INVALID_ARGUMENT;
+    }
+
+    double beta = impl->get_beta();
+
+    try {
+        Eigen::VectorXd tau_points;
+        if (impl->get_statistics() == SPIR_STATISTICS_FERMIONIC) {
+            auto ir_basis = _safe_static_pointer_cast<_IRBasis<sparseir::Fermionic>>(impl);
+            tau_points = sparseir::default_sampling_points(
+                *(ir_basis->get_impl()->sve_result->u), n_points
+            );
+            *n_points_returned = tau_points.size();
+        } else {
+            auto ir_basis = _safe_static_pointer_cast<_IRBasis<sparseir::Bosonic>>(impl);
+            tau_points = sparseir::default_sampling_points(
+                *(ir_basis->get_impl()->sve_result->u), n_points
+            );
+            *n_points_returned = tau_points.size();
+        }
+
+        // Copy the requested number of points
+        // rescale the points to the original domain
+        for (int i = 0; i < *n_points_returned; ++i) {
+            tau_points(i) = (tau_points(i) + 1) / 2 * beta;
+            if (tau_points(i) > 0.5 * beta) {
+                tau_points(i) -= beta;
+            }
+        }
+        std::copy(tau_points.data(), tau_points.data() + *n_points_returned, points);
+        return SPIR_COMPUTATION_SUCCESS;
+    } catch (const std::exception &e) {
+        return SPIR_GET_IMPL_FAILED;
+    }
+}
+
 int spir_basis_get_n_default_matsus(const spir_basis *b, bool positive_only, int *num_points)
 {
     if (!b || !num_points) {

test/cpp/cinterface_sampling.cxx

@@ -82,6 +82,18 @@ void test_tau_sampling()
     status = spir_basis_get_default_taus(basis, tau_points_org);
     REQUIRE(status == SPIR_COMPUTATION_SUCCESS);
 
+    int n_tau_points_ext = n_tau_points + 1;
+    double *tau_points_ext = (double *)malloc(n_tau_points_ext * sizeof(double));
+    int n_tau_points_returned;
+    status = spir_basis_get_default_taus_ext(
+        basis, n_tau_points_ext, tau_points_ext, &n_tau_points_returned);
+    REQUIRE(status == SPIR_COMPUTATION_SUCCESS);
+    REQUIRE(n_tau_points_returned == n_tau_points_ext);
+    for (int i = 0; i < n_tau_points_returned; i++) {
+        REQUIRE(tau_points_ext[i] >= -0.5 * beta);
+        REQUIRE(tau_points_ext[i] <=  0.5 * beta);
+    }
+
     int sampling_status;
     spir_sampling *sampling = spir_tau_sampling_new(basis, n_tau_points, tau_points_org, &sampling_status);
     REQUIRE(sampling_status == SPIR_COMPUTATION_SUCCESS);