diff --git a/c++/triqs/det_manip/det_manip.hpp b/c++/triqs/det_manip/det_manip.hpp index 2bf8c79d62..6c6295a5b3 100644 --- a/c++/triqs/det_manip/det_manip.hpp +++ b/c++/triqs/det_manip/det_manip.hpp @@ -32,10 +32,19 @@ #include #include +#include #include +#include +#include #include #include +#include +#include #include +#include +#include +#include +#include #include namespace triqs::det_manip { @@ -114,11 +123,9 @@ namespace triqs::det_manip { * storages. * * @param f triqs::det_manip::MatrixBuilder object. - * @param ncap Initial capacity for the size of the matrix, i.e. the maximum number of rows and columns. - * @param kcap Initial capacity for the maximum number of rows and columns that can be added or removed in a single - * operation. + * @param cap Initial capacity for the size of the matrix, i.e. the maximum number of rows and columns. */ - det_manip(F f, long ncap, long kcap = 1) : f_(std::move(f)) { reserve(ncap, kcap); } + det_manip(F f, long cap) : f_(std::move(f)) { reserve(cap); } /** * @brief Construct a det_manip object with a triqs::det_manip::MatrixBuilder and two ranges containing the @@ -132,22 +139,23 @@ namespace triqs::det_manip { */ template requires(MatrixBuilderXRange && MatrixBuilderYRange) - det_manip(F f, X &&x_rg, Y &&y_rg) : f_(std::move(f)), n_(std::ranges::size(x_rg)) { // NOLINT (ranges need not be forwarded) + det_manip(F f, X &&x_rg, Y &&y_rg) : f_(std::move(f)) { // NOLINT (ranges need not be forwarded) // check input sizes - if (n_ != std::ranges::size(y_rg)) TRIQS_RUNTIME_ERROR << "Error in det_manip::det_manip: Argument ranges have different sizes"; + auto const sz = std::ranges::size(x_rg); + if (sz != std::ranges::size(y_rg)) TRIQS_RUNTIME_ERROR << "Error in det_manip::det_manip: Argument ranges have different sizes"; // early return if the argument ranges are empty - if (n_ == 0) { + if (sz == 0) { reserve(30); return; } // reserve memory and fill the data storages - reserve(n_ * 2); + reserve(sz * 2); set_xy(x_rg, y_rg); // determinant and inverse matrix - auto M_v = M_(nda::range(size()), nda::range(size())); + auto M_v = M_(nda::range(sz), nda::range(sz)); nda::for_each(M_v.shape(), [this, &M_v](auto i, auto j) { M_v(i, j) = f_(x_[i], y_[j]); }); det_ = nda::determinant(M_v); M_v = nda::inverse(M_v); @@ -156,38 +164,27 @@ namespace triqs::det_manip { /** * @brief Reserve memory and resize the data storages. * - * @details It only reserves/resizes if the current capacities are smaller than the new capacities. + * @details It only reserves/resizes if the current capacities are smaller than the new capacities. It does not + * reserve memory for the working data. * - * @param new_ncap New capacity for the size of the matrix, i.e. the maximum number of rows and columns. - * @param new_kcap New capacity for the maximum number of rows and columns that can be added or removed in a single - * operation. + * @param cap New capacity for the size of the matrix, i.e. the maximum number of rows and columns. */ - void reserve(long new_ncap, long new_kcap = 1) { - if (new_kcap > kcap_) { - kcap_ = new_kcap; - if (new_ncap <= ncap_) wk_.resize(ncap_, kcap_); - } - if (new_ncap > ncap_) { - ncap_ = 2 * new_ncap; - + void reserve(long cap) { + if (cap > capacity()) { matrix_type M_copy(M_); - M_.resize(ncap_, ncap_); + M_.resize(cap, cap); auto rg = nda::range(M_copy.extent(0)); M_(rg, rg) = M_copy; - row_perm_.reserve(ncap_); - col_perm_.reserve(ncap_); - x_.reserve(ncap_); - y_.reserve(ncap_); - - w1_.resize(ncap_); - wk_.resize(ncap_, kcap_); + row_perm_.reserve(cap); + col_perm_.reserve(cap); + x_.reserve(cap); + y_.reserve(cap); } } - /// Clear the data storages and set the size to zero. + /// Clear the data storages and reset the matrix to size zero. void clear() { - n_ = 0; sign_ = 1; det_ = 1; last_try_ = try_tag::NoTry; @@ -236,7 +233,7 @@ namespace triqs::det_manip { void set_precision_error(double eps) { precision_error_ = eps; } /// Get the current size of the matrix. - [[nodiscard]] auto size() const { return n_; } + [[nodiscard]] auto size() const { return static_cast(x_.size()); } /// Get current capacity of the data storages. [[nodiscard]] auto capacity() const { return M_.shape()[0]; } @@ -268,16 +265,16 @@ namespace triqs::det_manip { /// Get a vector with all matrix builder arguments \f$ \mathbf{x} \f$. [[nodiscard]] auto get_x() const { std::vector res; - res.reserve(n_); - for (auto i : nda::range(n_)) res.emplace_back(x_[row_perm_[i]]); + res.reserve(size()); + for (auto i : nda::range{size()}) res.emplace_back(x_[row_perm_[i]]); return res; } /// Get a vector with all matrix builder arguments \f$ \mathbf{y} \f$. [[nodiscard]] auto get_y() const { std::vector res; - res.reserve(n_); - for (auto i : nda::range(n_)) res.emplace_back(y_[col_perm_[i]]); + res.reserve(size()); + for (auto i : nda::range{size()}) res.emplace_back(y_[col_perm_[i]]); return res; } @@ -467,22 +464,22 @@ namespace triqs::det_manip { value_type try_insert(long i, long j, x_type const &x, y_type const &y) { // check input arguments and copy them to the working data EXPECTS(last_try_ == try_tag::NoTry); - EXPECTS(0 <= i and i <= n_); - EXPECTS(0 <= j and j <= n_); + EXPECTS(0 <= i and i <= size()); + EXPECTS(0 <= j and j <= size()); std::tie(wins_.i, wins_.j, wins_.x, wins_.y) = std::make_tuple(i, j, x, y); // set the try tag last_try_ = try_tag::Insert; // early return if the current matrix is empty - if (n_ == 0) { + if (size() == 0) { newdet_ = f_(x, y); newsign_ = 1; return newdet_; } // reserve memory for the working data - if (n_ + 1 > wins_.capacity()) wins_.reserve(2 * (n_ + 1)); + if (size() + 1 > wins_.capacity()) wins_.reserve(2 * (size() + 1)); // calculate the new column B and the new row C of the matrix G (except for the element D) for (long k = 0; k < size(); ++k) { @@ -491,7 +488,7 @@ namespace triqs::det_manip { } // calculate S^{-1} = D - C * M * B - auto rg_n = nda::range(n_); + auto rg_n = nda::range(size()); nda::blas::gemv(1.0, M_(rg_n, rg_n), wins_.B(rg_n), 0.0, wins_.MB(rg_n)); wins_.S_inv = f_(x, y) - nda::blas::dot(wins_.C(rg_n), wins_.MB(rg_n)); @@ -505,9 +502,8 @@ namespace triqs::det_manip { private: // Complete the insert operation. void complete_insert() { - auto const new_size = n_ + 1; - auto const old_size = n_; - ++n_; + auto const new_size = size() + 1; + auto const old_size = size(); // reserve data storages if (new_size > capacity()) reserve(2 * new_size); @@ -517,7 +513,7 @@ namespace triqs::det_manip { y_.push_back(wins_.y); // early return if the new matrix is size 1 - if (n_ == 1) { + if (new_size == 1) { M_(0, 0) = 1 / newdet_; row_perm_.push_back(0); col_perm_.push_back(0); @@ -588,12 +584,12 @@ namespace triqs::det_manip { */ value_type try_insert_k(std::vector i, std::vector j, std::vector x, std::vector y) { // check input argument sizes - k_ = static_cast(i.size()); + auto const k = static_cast(i.size()); EXPECTS(last_try_ == try_tag::NoTry); - EXPECTS(k_ > 0); - EXPECTS(j.size() == k_); - EXPECTS(x.size() == k_); - EXPECTS(y.size() == k_); + EXPECTS(k > 0); + EXPECTS(j.size() == k); + EXPECTS(x.size() == k); + EXPECTS(y.size() == k); // move the input arguments to the working data winsk_.i = std::move(i); @@ -606,41 +602,41 @@ namespace triqs::det_manip { std::ranges::sort(std::ranges::zip_view(winsk_.i, winsk_.x), comp); std::ranges::sort(std::ranges::zip_view(winsk_.j, winsk_.y), comp); EXPECTS(std::ranges::adjacent_find(winsk_.i) == winsk_.i.end()); - EXPECTS(winsk_.i.front() >= 0 && winsk_.i.back() < n_ + k_); + EXPECTS(winsk_.i.front() >= 0 && winsk_.i.back() < size() + k); EXPECTS(std::ranges::adjacent_find(winsk_.j) == winsk_.j.end()); - EXPECTS(winsk_.j.front() >= 0 && winsk_.j.back() < n_ + k_); + EXPECTS(winsk_.j.front() >= 0 && winsk_.j.back() < size() + k); // set the try tag last_try_ = try_tag::InsertK; // reserve memory for the working data auto const [n_cap, k_cap] = winsk_.capacity(); - if (n_ + k_ > n_cap || k_ > k_cap) winsk_.reserve(2 * (n_ + k_), k_); + if (size() + k > n_cap || k > k_cap) winsk_.reserve(2 * (size() + k), k); // build the matrix D as part of S^{-1} = D - C M B - nda::for_each(std::array{k_, k_}, [this](auto l, auto m) { winsk_.S_inv(l, m) = f_(winsk_.x[l], winsk_.y[m]); }); + nda::for_each(std::array{k, k}, [this](auto l, auto m) { winsk_.S_inv(l, m) = f_(winsk_.x[l], winsk_.y[m]); }); // early return if the current matrix is empty - if (n_ == 0) { - newdet_ = detail::determinant(winsk_.S_inv, k_); + if (size() == 0) { + newdet_ = detail::determinant(winsk_.S_inv, k); newsign_ = 1; return newdet_; } // calculate the new columns B and the new rows C of the matrix G (except for the block matrix D) - for (long l = 0; l < n_; ++l) { - for (long m = 0; m < k_; ++m) { + for (long l = 0; l < size(); ++l) { + for (long m = 0; m < k; ++m) { winsk_.B(l, m) = f_(x_[l], winsk_.y[m]); winsk_.C(m, l) = f_(winsk_.x[m], y_[l]); } } // calculate S^{-1} = D - C M B and its determinant - auto rg_n = nda::range(n_); - auto rg_k = nda::range(k_); + auto rg_n = nda::range(size()); + auto rg_k = nda::range(k); nda::blas::gemm(1.0, M_(rg_n, rg_n), winsk_.B(rg_n, rg_k), 0.0, winsk_.MB(rg_n, rg_k)); nda::blas::gemm(-1.0, winsk_.C(rg_k, rg_n), winsk_.MB(rg_n, rg_k), 1.0, winsk_.S_inv(rg_k, rg_k)); - auto const det_S_inv = detail::determinant(winsk_.S_inv, k_); + auto const det_S_inv = detail::determinant(winsk_.S_inv, k); // calculate the new determinant = det(G^{(n)}) * det(S^{-1}) and sign = old sign * (-1)^{\sum_l i_l + j_l} newdet_ = det_ * det_S_inv; @@ -679,8 +675,9 @@ namespace triqs::det_manip { private: // Complete the insert_k operation. void complete_insert_k() { - auto const new_size = n_ + k_; - auto const old_size = n_; + auto const k = static_cast(winsk_.i.size()); + auto const new_size = size() + k; + auto const old_size = size(); // reserve data storages if (new_size > capacity()) reserve(2 * new_size); @@ -692,18 +689,17 @@ namespace triqs::det_manip { std::ranges::copy(std::ranges::iota_view(old_size, new_size), std::back_inserter(col_perm_)); // early return if the old matrix was empty - auto rg_k = nda::range(k_); + auto rg_k = nda::range(k); if (old_size == 0) { - n_ = new_size; M_(rg_k, rg_k) = nda::inverse(winsk_.S_inv(rg_k, rg_k)); return; } // update the permutation vectors - for (auto l : rg_k) { - ++n_; - std::rotate(row_perm_.begin() + winsk_.i[l], row_perm_.begin() + n_ - 1, row_perm_.begin() + n_); - std::rotate(col_perm_.begin() + winsk_.j[l], col_perm_.begin() + n_ - 1, col_perm_.begin() + n_); + for (auto tmp_sz = old_size + 1; auto l : rg_k) { + std::rotate(row_perm_.begin() + winsk_.i[l], row_perm_.begin() + tmp_sz - 1, row_perm_.begin() + tmp_sz); + std::rotate(col_perm_.begin() + winsk_.j[l], col_perm_.begin() + tmp_sz - 1, col_perm_.begin() + tmp_sz); + ++tmp_sz; } // calculate the matrix product C M and the matrix S @@ -742,16 +738,16 @@ namespace triqs::det_manip { value_type try_remove(long i, long j) { // check input arguments and copy them to the working data EXPECTS(last_try_ == try_tag::NoTry); - EXPECTS(i >= 0 and i < n_); - EXPECTS(j >= 0 and j < n_); + EXPECTS(i >= 0 and i < size()); + EXPECTS(j >= 0 and j < size()); std::tie(wrem_.i, wrem_.j, wrem_.ip, wrem_.jp) = std::make_tuple(i, j, row_perm_[i], col_perm_[j]); // set the try tag last_try_ = try_tag::Remove; // calculate the signs associated with P1, P2, P3 and P4 - int s_p1p2 = (wrem_.ip == n_ - 1 ? 1 : -1); - s_p1p2 = (wrem_.jp == n_ - 1 ? s_p1p2 : -s_p1p2); + int s_p1p2 = (wrem_.ip == size() - 1 ? 1 : -1); + s_p1p2 = (wrem_.jp == size() - 1 ? s_p1p2 : -s_p1p2); int s_p3p4 = ((i + j) % 2 == 0 ? 1 : -1); // set the diagonal element S @@ -768,39 +764,40 @@ namespace triqs::det_manip { // Complete the remove operation. void complete_remove() { // early return if the resulting matrix is empty - if (n_ == 1) { + if (size() == 1) { clear(); return; } // perform the P1 and P2 permutations by swapping the row and column to be removed with the last row and column - auto rg_n = nda::range{n_}; - if (wrem_.ip != n_ - 1) { + auto const old_size = size(); + auto const new_size = old_size - 1; + auto rg_n = nda::range{old_size}; + if (wrem_.ip != new_size) { // for M, we have to apply P1^T to the columns - deep_swap(M_(rg_n, wrem_.ip), M_(rg_n, n_ - 1)); + deep_swap(M_(rg_n, wrem_.ip), M_(rg_n, new_size)); // update the x arguments and the row permutation vector - x_[wrem_.ip] = x_[n_ - 1]; + x_[wrem_.ip] = x_[new_size]; auto it1 = std::ranges::find(row_perm_, wrem_.ip); - auto it2 = std::ranges::find(row_perm_, n_ - 1); + auto it2 = std::ranges::find(row_perm_, new_size); std::swap(*it1, *it2); } - if (wrem_.jp != n_ - 1) { + if (wrem_.jp != new_size) { // for M, we have to apply P2^T to the rows - deep_swap(M_(wrem_.jp, rg_n), M_(n_ - 1, rg_n)); + deep_swap(M_(wrem_.jp, rg_n), M_(new_size, rg_n)); // update the y arguments and the column permutation vector - y_[wrem_.jp] = y_[n_ - 1]; + y_[wrem_.jp] = y_[new_size]; auto it1 = std::ranges::find(col_perm_, wrem_.jp); - auto it2 = std::ranges::find(col_perm_, n_ - 1); + auto it2 = std::ranges::find(col_perm_, new_size); std::swap(*it1, *it2); } // update the size of the matrix - --n_; - rg_n = nda::range{n_}; + rg_n = nda::range{new_size}; // remove elements from the row and column permutation vectors and from the x and y arguments - std::ignore = std::ranges::remove(row_perm_, n_); - std::ignore = std::ranges::remove(col_perm_, n_); + std::ignore = std::ranges::remove(row_perm_, new_size); + std::ignore = std::ranges::remove(col_perm_, new_size); row_perm_.pop_back(); col_perm_.pop_back(); x_.pop_back(); @@ -808,11 +805,11 @@ namespace triqs::det_manip { // calculate -S^{-1} auto mS_inv = -1 / wrem_.S; - ASSERT(std::isfinite(std::abs(w1_.ksi))); + ASSERT(std::isfinite(std::abs(mS_inv))); // solve P = \widetilde{M}^{(n-1)} + \widetilde{M}^{(n-1)} B S C \widetilde{M}^{(n-1)} for \widetilde{M}^{(n-1)} // by using the fact that we know -\widetilde{M}^{(n-1)} B S, -S C \widetilde{M}^{(n-1)} and S^{-1} - blas::ger(mS_inv, M_(rg_n, n_), M_(n_, rg_n), M_(rg_n, rg_n)); + blas::ger(mS_inv, M_(rg_n, new_size), M_(new_size, rg_n), M_(rg_n, rg_n)); } public: @@ -881,27 +878,27 @@ namespace triqs::det_manip { */ value_type try_remove_k(std::vector i, std::vector j) { // check input argument sizes - k_ = static_cast(i.size()); + auto const k = static_cast(i.size()); EXPECTS(last_try_ == try_tag::NoTry); - EXPECTS(k_ > 0 && k_ <= n_); - EXPECTS(j.size() == k_); + EXPECTS(k > 0 && k <= size()); + EXPECTS(j.size() == k); // sort and check input arguments std::ranges::sort(i); std::ranges::sort(j); - EXPECTS(std::ranges::adjacent_find(i) == i.end() && i.front() >= 0 && i.back() < n_); - EXPECTS(std::ranges::adjacent_find(j) == j.end() && j.front() >= 0 && j.back() < n_); + EXPECTS(std::ranges::adjacent_find(i) == i.end() && i.front() >= 0 && i.back() < size()); + EXPECTS(std::ranges::adjacent_find(j) == j.end() && j.front() >= 0 && j.back() < size()); // set the try tag last_try_ = try_tag::RemoveK; // reserve memory for the working data - wremk_.reserve(k_); + wremk_.reserve(k); // move input arguments to the working data and get the corresponding row/column positions in the matrix G wremk_.i = std::move(i); wremk_.j = std::move(j); - for (long l = 0; l < k_; ++l) { + for (long l = 0; l < k; ++l) { wremk_.ip[l] = row_perm_[wremk_.i[l]]; wremk_.jp[l] = col_perm_[wremk_.j[l]]; } @@ -909,8 +906,8 @@ namespace triqs::det_manip { // compute the signs of the permutations P1, P2, P3, P4 and set the matrix S int s_p1p2 = 1; long idx_sum = 0; - long target = n_ - k_; - for (long l = 0; l < k_; ++l) { + long target = size() - k; + for (long l = 0; l < k; ++l) { // the combined sign of P3 and P4 is simply (-1)^{\sum i_k + j_k} idx_sum += wremk_.i[l] + wremk_.j[l]; @@ -919,8 +916,8 @@ namespace triqs::det_manip { // if not, P1 has to swap it with the corresponding row s_p1p2 = -s_p1p2; // we have to take care of the case where the row is swapped with another row that we want to remove - auto it = std::find(wremk_.ip.begin() + l + 1, wremk_.ip.begin() + k_, target); - if (it != wremk_.ip.begin() + k_) { + auto it = std::find(wremk_.ip.begin() + l + 1, wremk_.ip.begin() + k, target); + if (it != wremk_.ip.begin() + k) { std::swap(wremk_.ip[l], *it); } else { wremk_.ip[l] = target; @@ -932,8 +929,8 @@ namespace triqs::det_manip { // if not, P2 has to swap it with the corresponding column s_p1p2 = -s_p1p2; // we have to take care of the case where the column is swapped with another column that we want to remove - auto it = std::find(wremk_.jp.begin() + l + 1, wremk_.jp.begin() + k_, target); - if (it != wremk_.jp.begin() + k_) { + auto it = std::find(wremk_.jp.begin() + l + 1, wremk_.jp.begin() + k, target); + if (it != wremk_.jp.begin() + k) { std::swap(wremk_.jp[l], *it); } else { wremk_.jp[l] = target; @@ -942,12 +939,12 @@ namespace triqs::det_manip { ++target; // set the elements of the matrix S - for (long m = 0; m < k_; ++m) { wremk_.S(l, m) = M_(col_perm_[wremk_.j[l]], row_perm_[wremk_.i[m]]); } + for (long m = 0; m < k; ++m) { wremk_.S(l, m) = M_(col_perm_[wremk_.j[l]], row_perm_[wremk_.i[m]]); } } int s_p3p4 = (idx_sum % 2 == 0 ? 1 : -1); // compute the new determinant and sign - auto det_S = detail::determinant(wremk_.S, k_); + auto det_S = detail::determinant(wremk_.S, k); newdet_ = det_ * det_S * s_p1p2; newsign_ = sign_ * s_p1p2 * s_p3p4; @@ -973,15 +970,19 @@ namespace triqs::det_manip { private: // Complete the remove_k operation. void complete_remove_k() { + auto const k = static_cast(wremk_.i.size()); + auto const new_size = size() - k; + auto const old_size = size(); + // early return if the resulting matrix is empty - if (n_ == k_) { + if (new_size == 0) { clear(); return; } // perform the P1 and P2 permutations by swapping the rows and columns accordingly - auto rg_n = nda::range{n_}; - for (long m = 0, target = n_ - k_; m < k_; ++m, ++target) { + auto rg_n = nda::range{old_size}; + for (long m = 0, target = new_size; m < k; ++m, ++target) { if (row_perm_[wremk_.i[m]] != target) { // for M, we have to apply P1^T to the columns deep_swap(M_(rg_n, row_perm_[wremk_.i[m]]), M_(rg_n, target)); @@ -1002,27 +1003,24 @@ namespace triqs::det_manip { } } - // update the size of the matrix - n_ -= k_; - rg_n = nda::range{n_}; - // remove elements from the row and column permutation vectors and from the x and y arguments - auto ge_n = [this](auto i) { return i >= n_; }; + auto ge_n = [new_size](auto i) { return i >= new_size; }; std::ignore = std::ranges::remove_if(row_perm_, ge_n); std::ignore = std::ranges::remove_if(col_perm_, ge_n); - row_perm_.resize(n_); - col_perm_.resize(n_); - x_.resize(n_); - y_.resize(n_); + row_perm_.resize(new_size); + col_perm_.resize(new_size); + x_.resize(new_size); + y_.resize(new_size); // calculate S^{-1} - auto rg_k = nda::range{k_}; - auto rg_n_nk = nda::range{n_, n_ + k_}; + auto rg_nk = nda::range{new_size}; + auto rg_k = nda::range{k}; + auto rg_nk_n = nda::range{new_size, old_size}; nda::inverse_in_place(wremk_.S(rg_k, rg_k)); // solve P = \widetilde{M}^{(n-k)} + \widetilde{M}^{(n-k)} B S C \widetilde{M}^{(n-k)} for \widetilde{M}^{(n-k)} // by using the fact that we know -\widetilde{M}^{(n-k)} B S, -S C \widetilde{M}^{(n-k)} and S^{-1} - blas::gemm(-1.0, M_(rg_n, rg_n_nk), wremk_.S(rg_k, rg_k) * M_(rg_n_nk, rg_n), 1.0, M_(rg_n, rg_n)); + blas::gemm(-1.0, M_(rg_nk, rg_nk_n), wremk_.S(rg_k, rg_k) * M_(rg_nk_n, rg_nk), 1.0, M_(rg_nk, rg_nk)); } // Complete the remove2 operation. @@ -1156,83 +1154,6 @@ namespace triqs::det_manip { nda::blas::ger(-1 / wrow_.xi, wrow_.Mu(rg), wrow_.vTM(rg), M_(rg, rg)); } - public: - /** - * Consider the change the row i and column j and the corresponding x and y - * - * Returns the ratio of det Minv_new / det Minv. - * This routine does NOT make any modification. It has to be completed with complete_operation(). - */ - value_type try_change_col_row(long i, long j, x_type const &x, y_type const &y) { - TRIQS_ASSERT(last_try_ == try_tag::NoTry); - TRIQS_ASSERT(0 <= i and i < n_); - TRIQS_ASSERT(0 <= j and j < n_); - - last_try_ = try_tag::ChangeRowCol; - w1_.i = i; - w1_.j = j; - w1_.ireal = row_perm_[i]; - w1_.jreal = col_perm_[j]; - w1_.x = x; - w1_.y = y; - - // Compute the col B. - for (long idx = 0; idx < n_; idx++) { // MC : delta_x, MB : delta_y - w1_.MC(idx) = f_(x_[idx], y) - f_(x_[idx], y_[w1_.jreal]); - w1_.MB(idx) = f_(x, y_[idx]) - f_(x_[w1_.ireal], y_[idx]); - } - w1_.MC(w1_.ireal) = f_(x, y) - f_(x_[w1_.ireal], y_[w1_.jreal]); - w1_.MB(w1_.jreal) = 0; - - nda::range RN(n_); - // C : X, B : Y - //w1.C(R) = mat_inv(R,R) * w1.MC(R);// OPTIMIZE BELOW - blas::gemv(1.0, M_(RN, RN), w1_.MC(RN), 0.0, w1_.C(RN)); - //w1.B(R) = transpose(mat_inv(R,R)) * w1.MB(R); // OPTIMIZE BELOW - blas::gemv(1.0, transpose(M_(RN, RN)), w1_.MB(RN), 0.0, w1_.B(RN)); - - // compute the det_ratio - auto Xn = w1_.C(w1_.jreal); - auto Yn = w1_.B(w1_.ireal); - auto Z = nda::blas::dot(w1_.MB(RN), w1_.C(RN)); - auto Mnn = M_(w1_.jreal, w1_.ireal); - auto det_ratio = (1 + Xn) * (1 + Yn) - Mnn * Z; - w1_.ksi = det_ratio; - newdet_ = det_ * det_ratio; - newsign_ = sign_; - return det_ratio; // newsign/sign is unity - } - //------------------------------------------------------------------------------------------ - private: - void complete_change_col_row() { - nda::range RN(n_); - x_[w1_.ireal] = w1_.x; - y_[w1_.jreal] = w1_.y; - - // FIXME : Use blas for this ? Is it better - auto Xn = w1_.C(w1_.jreal); - auto Yn = w1_.B(w1_.ireal); - auto Mnn = M_(w1_.jreal, w1_.ireal); - - auto D = w1_.ksi; // get back - auto a = -(1 + Yn) / D; // D in the notes - auto b = -(1 + Xn) / D; - auto Z = nda::blas::dot(w1_.MB(RN), w1_.C(RN)); - Z = Z / D; - Mnn = Mnn / D; - w1_.MB(RN) = M_(w1_.jreal, RN); // Mnj - w1_.MC(RN) = M_(RN, w1_.ireal); // Min - - for (long i = 0; i < n_; ++i) - for (long j = 0; j < n_; ++j) { - auto Xi = w1_.C(i); - auto Yj = w1_.B(j); - auto Mnj = w1_.MB(j); - auto Min = w1_.MC(i); - M_(i, j) += a * Xi * Mnj + b * Min * Yj + Mnn * Xi * Yj + Z * Min * Mnj; - } - } - public: /** * @brief Try to fill the original matrix \f$ F^{(n)} \f$ with new elements. @@ -1299,34 +1220,40 @@ namespace triqs::det_manip { } // reserve memory and reset the data storages and permutation vectors - n_ = wref_.size(); if (wref_.size() > capacity()) reserve(2 * wref_.size()); set_xy(wref_.x, wref_.y); // set the new inverse matrix M - auto rg = nda::range(n_); + auto rg = nda::range(size()); M_(rg, rg) = nda::inverse(wref_.G(rg, rg)); } public: /** - * Finish the move of the last try_xxx called. - * Throws if no try_xxx has been done or if the last operation was complete_operation. - */ + * @brief Complete the last try-operation. + * + * @details It completes the last try-operation by calling the correct completion function depending on the try tag + * set in the last try function call. + * + * If the number of operations exceeds the threshold set in the triqs::det_manip::param_type object, the inverse + * matrix \f$ M^{(n)} \f$, the determinant \f$ \det(G^{(n)}) \f$ and the sign \f$ s^{(n)} \f$ are regenerated using + * the matrix builder (see regenerate_and_check()). + * + * A possible warning is emitted or an exception is thrown if the current objects are not consistent with the + * regenerated ones. + */ void complete_operation() { switch (last_try_) { case (try_tag::Insert): complete_insert(); break; case (try_tag::Remove): complete_remove(); break; case (try_tag::ChangeCol): complete_change_col(); break; case (try_tag::ChangeRow): complete_change_row(); break; - case (try_tag::ChangeRowCol): complete_change_col_row(); break; case (try_tag::InsertK): complete_insert_k(); break; case (try_tag::RemoveK): complete_remove_k(); break; case (try_tag::Refill): complete_refill(); break; case (try_tag::NoTry): return; break; - default: TRIQS_RUNTIME_ERROR << "Misuing det_manip"; // Never used? + default: return; } - det_ = newdet_; sign_ = newsign_; ++nops_; @@ -1334,77 +1261,9 @@ namespace triqs::det_manip { last_try_ = try_tag::NoTry; } - /** - * Reject the previous try_xxx called. - * All try_xxx have to be either accepted (complete_operation) or rejected. - */ + /// Reject the last try-operation. void reject_last_try() { last_try_ = try_tag::NoTry; } - // ----------------- A few short cuts ----------------- - - public: - /// Insert (try_insert + complete) - value_type insert(long i, long j, x_type const &x, y_type const &y) { - auto r = try_insert(i, j, x, y); - complete_operation(); - return r; - } - - /// Insert_at_end (try_insert + complete) - value_type insert_at_end(x_type const &x, y_type const &y) { return insert(n_, n_, x, y); } - - /// Insert2 (try_insert2 + complete) - value_type insert2(long i0, long i1, long j0, long j1, x_type const &x0, x_type const &x1, y_type const &y0, y_type const &y1) { - auto r = try_insert2(i0, i1, j0, j1, x0, x1, y0, y1); - complete_operation(); - return r; - } - - /// Insert2_at_end (try_insert2 + complete) - value_type insert2_at_end(x_type const &x0, x_type const &x1, y_type const &y0, y_type const &y1) { - return insert2(n_, n_ + 1, n_, n_ + 1, x0, x1, y0, y1); - } - - /// Remove (try_remove + complete) - value_type remove(long i, long j) { - auto r = try_remove(i, j); - complete_operation(); - return r; - } - - /// Remove_at_end (try_remove + complete) - value_type remove_at_end() { return remove(n_ - 1, n_ - 1); } - - /// Remove2 (try_remove2 + complete) - value_type remove2(long i0, long i1, long j0, long j1) { - auto r = try_remove2(i0, i1, j0, j1); - complete_operation(); - return r; - } - - /// Remove2_at_end (try_remove2 + complete) - value_type remove2_at_end() { return remove2(n_ - 1, n_ - 2, n_ - 1, n_ - 2); } - - /// change_col (try_change_col + complete) - value_type change_col(long j, y_type const &y) { - auto r = try_change_col(j, y); - complete_operation(); - return r; - } - - /// change_row (try_change_row + complete) - value_type change_row(long i, x_type const &x) { - auto r = try_change_row(i, x); - complete_operation(); - return r; - } - - value_type change_one_row_and_one_col(long i, long j, x_type const &x, y_type const &y) { - auto r = try_change_col_row(i, j, x, y); - complete_operation(); - return r; - } - /** * @brief Regenerate the inverse matrix, determinant and sign and check the consistency of the current * values/objects. @@ -1476,43 +1335,60 @@ namespace triqs::det_manip { sign_ = (exp_sign > 0 ? 1 : -1); } - /// Write into HDF5 - friend void h5_write(h5::group fg, std::string subgroup_name, det_manip const &g) { - auto gr = fg.create_group(subgroup_name); - h5_write(gr, "N", g.n_); - h5_write(gr, "mat_inv", g.M_); - h5_write(gr, "det", g.det_); - h5_write(gr, "sign", g.sign_); - h5_write(gr, "row_num", g.row_perm_); - h5_write(gr, "col_num", g.col_perm_); - h5_write(gr, "x_values", g.x_); - h5_write(gr, "y_values", g.y_); - h5_write(gr, "n_opts", g.nops_); - h5_write(gr, "n_opts_max_before_check", g.nops_before_check_); - h5_write(gr, "singular_threshold", g.singular_threshold_); + /** + * @brief Write a triqs::det_manip::det_manip object to HDF5. + * + * @param g h5::group containing the subgroup to be written to. + * @param name Name of the subgroup. + * @param dm Manipulator object to be written. + */ + friend void h5_write(h5::group g, std::string name, det_manip const &dm) { + auto gr = g.create_group(name); + h5::write(gr, "x", dm.x_); + h5::write(gr, "y", dm.y_); + h5::write(gr, "M", dm.M_); + h5::write(gr, "det", dm.det_); + h5::write(gr, "row_perm", dm.row_perm_); + h5::write(gr, "col_perm", dm.col_perm_); + h5::write(gr, "sign", dm.sign_); + h5::write(gr, "nops_before_check", dm.nops_before_check_); + h5::write(gr, "singular_threshold", dm.singular_threshold_); + h5::write(gr, "precision_warning", dm.precision_warning_); + h5::write(gr, "precision_error", dm.precision_error_); + h5::write(gr, "nops", dm.nops_); } - /// Read from HDF5 - friend void h5_read(h5::group fg, std::string subgroup_name, det_manip &g) { - auto gr = fg.open_group(subgroup_name); - h5_read(gr, "N", g.n_); - h5_read(gr, "mat_inv", g.M_); - g.ncap_ = first_dim(g.M_); // restore Nmax - g.last_try_ = try_tag::NoTry; - h5_read(gr, "det", g.det_); - h5_read(gr, "sign", g.sign_); - h5_read(gr, "row_num", g.row_perm_); - h5_read(gr, "col_num", g.col_perm_); - h5_read(gr, "x_values", g.x_); - h5_read(gr, "y_values", g.y_); - h5_read(gr, "n_opts", g.nops_); - h5_read(gr, "n_opts_max_before_check", g.nops_before_check_); - h5_read(gr, "singular_threshold", g.singular_threshold_); + /** + * @brief Read a triqs::det_manip::det_manip object from HDF5. + * + * @param g h5::group containing the subgroup to be read from. + * @param name Name of the subgroup. + * @param dm Manipulator object to be read into. + */ + friend void h5_read(h5::group g, std::string name, det_manip &dm) { + auto gr = g.open_group(name); + h5::read(gr, "x", dm.x_); + h5::read(gr, "y", dm.y_); + h5::read(gr, "M", dm.M_); + h5::read(gr, "det", dm.det_); + h5::read(gr, "row_perm", dm.row_perm_); + h5::read(gr, "col_perm", dm.col_perm_); + h5::read(gr, "sign", dm.sign_); + h5::read(gr, "nops_before_check", dm.nops_before_check_); + h5::read(gr, "singular_threshold", dm.singular_threshold_); + h5::read(gr, "precision_warning", dm.precision_warning_); + h5::read(gr, "precision_error", dm.precision_error_); + h5::read(gr, "nops", dm.nops_); + dm.x_.reserve(dm.capacity()); + dm.y_.reserve(dm.capacity()); + dm.row_perm_.reserve(dm.capacity()); + dm.col_perm_.reserve(dm.capacity()); + dm.last_try_ = try_tag::NoTry; } private: // Enumerate the different operations supported by the det_manip class that have a try - complete step. - enum class try_tag { NoTry, Insert, Remove, ChangeCol, ChangeRow, ChangeRowCol, InsertK, RemoveK, Refill }; + enum class try_tag { NoTry, Insert, Remove, ChangeCol, ChangeRow, InsertK, RemoveK, Refill }; // Set the matrix builder arguments to the given ranges and reset the permutation vectors. template @@ -1560,8 +1436,6 @@ namespace triqs::det_manip { detail::work_data_change_col wcol_; detail::work_data_change_row wrow_; detail::work_data_refill wref_; - detail::work_data_type1 w1_; - detail::work_data_typek wk_; value_type newdet_{1}; int newsign_{1}; @@ -1574,11 +1448,6 @@ namespace triqs::det_manip { // tag and operation counter try_tag last_try_{try_tag::NoTry}; std::uint64_t nops_{0}; - - // sizes of matrices and capacities of their data storages - long n_{0}; - long ncap_{0}; - long k_{0}; - long kcap_{1}; }; + } // namespace triqs::det_manip diff --git a/c++/triqs/det_manip/work_data.hpp b/c++/triqs/det_manip/work_data.hpp index 5d6bc7551d..6d2361f87d 100644 --- a/c++/triqs/det_manip/work_data.hpp +++ b/c++/triqs/det_manip/work_data.hpp @@ -26,78 +26,6 @@ namespace triqs::det_manip::detail { - // ================ Work Data Types ===================== - - // For single-row/column operations - template struct work_data_type1 { - x_type x; - y_type y; - long i, j, ireal, jreal; - // MB = A^(-1)*B, - // MC = C*A^(-1) - nda::vector MB, MC, B, C; - // ksi = newdet/det - value_type ksi; - void resize(long N) { - MB.resize(N); - MC.resize(N); - B.resize(N); - C.resize(N); - } - }; - - // For multiple-row/column operations - template struct work_data_typek { - std::vector x; - std::vector y; - std::vector i, j, ireal, jreal; - // MB = A^(-1)*B, - // MC = C*A^(-1) - nda::matrix MB, MC, B, C, ksi; - void resize(long N, long k) { - if (k < 2) return; - x.resize(k); - y.resize(k); - i.resize(k); - j.resize(k); - ireal.resize(k); - jreal.resize(k); - MB.resize(N, k); - MC.resize(k, N); - B.resize(N, k); - C.resize(k, N); - ksi.resize(k, k); - } - value_type det_ksi(long k) const { - if (k == 2) { - return ksi(0, 0) * ksi(1, 1) - ksi(1, 0) * ksi(0, 1); - } else if (k == 3) { - return // Rule of Sarrus - ksi(0, 0) * ksi(1, 1) * ksi(2, 2) + // - ksi(0, 1) * ksi(1, 2) * ksi(2, 0) + // - ksi(0, 2) * ksi(1, 0) * ksi(2, 1) - // - ksi(2, 0) * ksi(1, 1) * ksi(0, 2) - // - ksi(2, 1) * ksi(1, 2) * ksi(0, 0) - // - ksi(2, 2) * ksi(1, 0) * ksi(0, 1); // - } else { - auto Rk = nda::range(k); - return nda::determinant(ksi(Rk, Rk)); - }; - } - }; - - // For refill operations - template struct work_data_type_refill { - std::vector x_values; - std::vector y_values; - nda::matrix M; - void reserve(long N) { - x_values.reserve(N); - y_values.reserve(N); - M.resize(N, N); - } - }; - // Data storage for temporary data used in the det_manip class when inserting a new row and column. // // - x and y: MatrixBuilder arguments for the new row and column. @@ -108,13 +36,9 @@ namespace triqs::det_manip::detail { template struct work_data_insert { X x; Y y; - long i; - long j; + long i, j; T S_inv; - nda::vector B; - nda::vector C; - nda::vector MB; - nda::vector CM; + nda::vector B, C, MB, CM; // Get current capacity of the data storages. auto capacity() const { return B.size(); } @@ -143,13 +67,8 @@ namespace triqs::det_manip::detail { template struct work_data_insert_k { std::vector x; std::vector y; - std::vector i; - std::vector j; - nda::matrix S_inv; - nda::matrix B; - nda::matrix C; - nda::matrix MB; - nda::matrix CM; + std::vector i, j; + nda::matrix S_inv, B, C, MB, CM; // Get current capacity of the data storages. auto capacity() const { return std::make_pair(B.shape()[0], B.shape()[1]); } @@ -173,10 +92,7 @@ namespace triqs::det_manip::detail { // - ip and jp: Positions of the row and column in the matrix G. // - S: Diagonal element of \widetilde{M}^{(n)}. template struct work_data_remove { - long i; - long j; - long ip; - long jp; + long i, j, ip, jp; T S; }; @@ -186,10 +102,7 @@ namespace triqs::det_manip::detail { // - ip and jp: Positions of the rows and columns in the matrix G. // - S: Block matrix of \widetilde{M}^{(n)}. template struct work_data_remove_k { - std::vector i; - std::vector j; - std::vector ip; - std::vector jp; + std::vector i, j, ip, jp; nda::matrix S; // Get current capacity of the data storages. @@ -216,11 +129,8 @@ namespace triqs::det_manip::detail { // - xi: Factor appearing in the matrix determinant lemma, i.e. xi = (1 + v^T M u). template struct work_data_change_col { Y y; - long j; - long jp; - nda::vector u; - nda::vector Mu; - nda::vector vTM; + long j, jp; + nda::vector u, Mu, vTM; T xi; // Get current capacity of the data storages. @@ -247,11 +157,8 @@ namespace triqs::det_manip::detail { // - xi: Factor appearing in the matrix determinant lemma, i.e. xi = (1 + v^T M u). template struct work_data_change_row { X x; - long i; - long ip; - nda::vector vT; - nda::vector vTM; - nda::vector Mu; + long i, ip; + nda::vector vT, vTM, Mu; T xi; // Get current capacity of the data storages. diff --git a/test/c++/det_manip/det_manip1.cpp b/test/c++/det_manip/det_manip1.cpp deleted file mode 100644 index 11ec77c2fd..0000000000 --- a/test/c++/det_manip/det_manip1.cpp +++ /dev/null @@ -1,190 +0,0 @@ -// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) -// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS) -// Copyright (c) 2018-2023 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Henri Menke, Olivier Parcollet, Nils Wentzell - -#include -#include -#include -#include -#include "./old_test_tool.hpp" - -struct fun { - - typedef double result_type; - typedef double argument_type; - - /*double f(double eps, double tau, double beta) const { - bool s = (tau>0); - tau = (s ? tau : beta + tau); - return (s ? 1 : -1 ) * exp( - eps* (tau)) / (1 + exp(- beta *eps)); - }*/ - - double operator()(double x, double y) const { - const double pi = acos(-1); - const double beta = 10.0; - const double epsi = 0.1; - double tau = x - y; - bool s = (tau > 0); - tau = (s ? tau : beta + tau); - double r = epsi + tau / beta * (1 - 2 * epsi); - return -2 * (pi / beta) / std::sin(pi * r); - - //return - 0.5*(f(0.12, x-y, beta) + f(0.43,x-y,beta) ); //+ f(-0.7,x-y,beta)); - } -}; - -template void assert_close(T1 const &A, T2 const &B, double precision) { - if (std::abs(A - B) > precision) TRIQS_RUNTIME_ERROR << "assert_close error : " << A << "\n" << B; -} -const double PRECISION = 1.e-6; - -struct test { - - fun f; - triqs::det_manip::det_manip D; - double det_old, detratio; - - test() : f(), D(f, 100) {} - - //#define PRINT_ALL - void check() { - std::cerr << "---- check --- " << std::endl; -#ifndef PRINT_ALL - std::cerr << "det = " << D.determinant() << " == " << double(determinant(D.matrix())) << std::endl; -#else - std::cerr << "det = " << D.determinant() << " == " << double(determinant(D.matrix())) << std::endl - << D.inverse_matrix() << D.matrix() << nda::matrix(inverse(D.matrix())) << std::endl; - std::cerr << "det_old = " << det_old << "detratio = " << detratio << " determin " << D.determinant() << std::endl; -#endif - auto diff = nda::matrix(inverse(D.matrix()) - D.inverse_matrix()); - //std::cerr << diff < 0) detratio = D.try_remove(RNG(s), RNG(s)); - break; - case 2: { - long k = 2 + RNG(6); - std::cerr << " Insert" << k << std::endl; - std::vector x_vec(k), y_vec(k); - std::vector i_vec(k), j_vec(k); - for (auto m : itertools::range(k)) { - x_vec.at(m) = RNG(10.0); - y_vec.at(m) = RNG(10.0); - i_vec.at(m) = RNG(s + m); - j_vec.at(m) = RNG(s + m); - } - std::sort(i_vec.begin(), i_vec.end()); - std::sort(j_vec.begin(), j_vec.end()); - if (std::unique(i_vec.begin(), i_vec.end()) == i_vec.end() && std::unique(j_vec.begin(), j_vec.end()) == j_vec.end()) { - detratio = D.try_insert_k(i_vec, j_vec, x_vec, y_vec); - } else - do_something = false; - break; - } - case 3: { - long k = 2 + RNG(6); - std::cerr << " Remove" << k << std::endl; - if (D.size() >= k) { - std::vector i_vec(k), j_vec(k); - for (auto m : itertools::range(k)) { - i_vec.at(m) = RNG(s); - j_vec.at(m) = RNG(s); - } - std::sort(i_vec.begin(), i_vec.end()); - std::sort(j_vec.begin(), j_vec.end()); - if (std::unique(i_vec.begin(), i_vec.end()) == i_vec.end() && std::unique(j_vec.begin(), j_vec.end()) == j_vec.end()) { - detratio = D.try_remove_k(i_vec, j_vec); - } else - do_something = false; - } - break; - } - case 4: - if (D.size() == 0) break; - y = RNG(10.0); - i0 = RNG(s); - std::cerr << " try_change_col" << i0 << std::endl; - detratio = D.try_change_col(i0, y); - break; - case 5: - if (D.size() == 0) break; - y = RNG(10.0); - i0 = RNG(s); - std::cerr << " try_change_row" << i0 << std::endl; - detratio = D.try_change_row(i0, y); - break; - case 6: - if (D.size() == 0) break; - x = RNG(10.0); - y = RNG(10.0); - i0 = RNG(s); - j0 = RNG(s); - - std::cerr << " try_change_col_row" << i0 << " " << j0 << std::endl; - - //D.try_change_col(j0,y); - //D.try_change_row(i0,x); - //D.reject_last_try(); - - detratio = D.try_change_col_row(i0, j0, x, y); - break; - default: TRIQS_RUNTIME_ERROR << " TEST INTERNAL ERROR"; - }; - - if (do_something) { - if (std::abs(detratio * det_old) > 1.e-3) { - D.complete_operation(); - if (D.size() > 0) check(); - } else { - std::cerr << " reject since new det is = " << std::abs(detratio * det_old) << std::endl; - D.reject_last_try(); - } - } - } - } -}; - -int main(int, char **) { test().run(); } diff --git a/test/c++/det_manip/det_manip2.cpp b/test/c++/det_manip/det_manip2.cpp deleted file mode 100644 index 2f7f7994df..0000000000 --- a/test/c++/det_manip/det_manip2.cpp +++ /dev/null @@ -1,125 +0,0 @@ -// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) -// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS) -// Copyright (c) 2018-2023 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Henri Menke, Olivier Parcollet, Nils Wentzell - -#include -#include -#include -#include "./old_test_tool.hpp" - -#include - -struct fun { - typedef double result_type; - typedef double argument_type; - double operator()(double x, double y) const { return 0.5 * (y >= x ? 1 : -1); } -}; - -template void assert_close(T1 const &A, T2 const &B, double precision) { - if (std::abs(A - B) > precision) TRIQS_RUNTIME_ERROR << "assert_close error : " << A << "\n" << B; -} -const double PRECISION = 1.e-12; - -struct test { - - fun f; - triqs::det_manip::det_manip D; - double det_old, detratio; - - test() : f(), D(f, 100) {} - - void check() { -#define PRINT_ALLL -#ifndef PRINT_ALL - std::cerr << "det = " << D.determinant() << " == " << double(determinant(D.matrix())) << std::endl; -#else - std::cerr << "det = " << D.determinant() << " == " << double(determinant(D.matrix())) << std::endl; - std::cerr << "inverse_matrix = " << D.inverse_matrix() << std::endl; - //std::cerr << "inverse inverse_matrix = " << nda::matrix(inverse(D.inverse_matrix()) ) << std::endl; - std::cerr << "matrix = " << D.matrix() << std::endl; - std::cerr << nda::matrix(inverse(D.matrix())) << std::endl; -#endif - assert_close(D.determinant(), 0.5, PRECISION); - assert_close(D.determinant(), 1 / determinant(D.inverse_matrix()), PRECISION); - nda::assert_all_close(inverse(D.matrix()), D.inverse_matrix(), PRECISION); - assert_close(det_old * detratio, D.determinant(), PRECISION); - } - - void run() { - triqs::mc_tools::random_generator RNG("mt19937", 23432); - for (size_t i = 0; i < 100; ++i) { - std::cerr << " ------------------------------------------------" << std::endl; - std::cerr << " i = " << i << " size = " << D.size() << std::endl; - // choose a move - size_t s = D.size(); - size_t w; - det_old = D.determinant(); - detratio = 1; - double x; - - switch (RNG((i > 10 ? 4 : 1))) { - case 0: - std::cout << "Insert" << std::endl; - x = RNG(10.0); - w = RNG(s); - detratio = D.try_insert(w, w, x, x); - break; - case 1: - std::cout << "Remove" << std::endl; - w = RNG(s); - if (s > 0) detratio = D.try_remove(w, w); - break; - case 2: { - long k = 2 + RNG(6); - std::vector x_vec(k); - std::vector w_vec(k); - for (auto m : itertools::range(k)) { - x_vec.at(m) = RNG(10.0); - w_vec.at(m) = RNG(s + m); - } - std::sort(w_vec.begin(), w_vec.end()); - if (std::unique(w_vec.begin(), w_vec.end()) == w_vec.end()) { - std::cout << " Insert" << k << std::endl; - detratio = D.try_insert_k(w_vec, w_vec, x_vec, x_vec); - } - break; - } - case 3: { - long k = 2 + RNG(6); - std::cout << " Remove" << k << std::endl; - if (D.size() >= 2) { - std::vector w_vec(k); - for (auto m : itertools::range(k)) { w_vec.at(m) = RNG(s); } - std::sort(w_vec.begin(), w_vec.end()); - if (std::unique(w_vec.begin(), w_vec.end()) == w_vec.end()) detratio = D.try_remove_k(w_vec, w_vec); - } - break; - } - default: TRIQS_RUNTIME_ERROR << " TEST INTERNAL ERROR"; - }; - - std::cout << "about to complete ..." << std::endl; - D.complete_operation(); - std::cout << " .. done" << std::endl; - if (D.size() > 0) check(); - std::cout << " check done" << std::endl; - } - } -}; - -int main(int, char **) { test().run(); } diff --git a/test/c++/det_manip/det_manip3.cpp b/test/c++/det_manip/det_manip3.cpp deleted file mode 100644 index 90bb8d8559..0000000000 --- a/test/c++/det_manip/det_manip3.cpp +++ /dev/null @@ -1,97 +0,0 @@ -// Copyright (c) 2016-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) -// Copyright (c) 2016-2018 Centre national de la recherche scientifique (CNRS) -// Copyright (c) 2018-2021 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Olivier Parcollet, Nils Wentzell - -#include -#include -#include -#include -namespace arrays = nda; -using triqs::det_manip::det_manip; - -struct func { - // gives the coefficients of the matrix (function F of the documentation) - double operator()(int x, int y) const { - if ((x < 0) and (y < 0)) - return 1; - else - return x + y; - } -}; - -det_manip init_dm() { - func f; - std::vector initial_x{-2, 2}, initial_y{-5, 3}; - det_manip dm(f, initial_x, initial_y); - std::cerr << "matrix = " << dm.matrix() << std::endl; - std::cerr << "det = " << dm.determinant() << std::endl; - return dm; -} - -det_manip remove_second_row_col(det_manip dm) { - dm.remove(1, 1); - std::cerr << "matrix = " << dm.matrix() << std::endl; - std::cerr << "det = " << dm.determinant() << std::endl; - return dm; -} - -det_manip insert_second_row_col(det_manip dm) { - dm.insert(1, 1, 6, 4); - std::cerr << "matrix = " << dm.matrix() << std::endl; - std::cerr << "det = " << dm.determinant() << std::endl; - return dm; -} - -det_manip remove_first_row_col(det_manip dm) { - dm.remove(0, 0); - std::cerr << "matrix = " << dm.matrix() << std::endl; - std::cerr << "det = " << dm.determinant() << std::endl; - return dm; -} - -// ----- TESTS ------------------ - -TEST(det_manip, det_manip_zero_mat) { - - // using std::abs; - - auto dm1 = init_dm(); - auto dm2 = remove_second_row_col(dm1); - auto dm3 = insert_second_row_col(dm2); - auto dm4 = remove_first_row_col(dm2); - nda::matrix true_dm1 = {{1, 1}, {-3, 5}}; - nda::matrix true_dm2 = {{1}}; - nda::matrix true_dm3 = {{1, 2}, {1, 10}}; - nda::matrix true_dm4 = {}; - - EXPECT_ARRAY_NEAR(dm1.matrix(), true_dm1); - EXPECT_EQ(dm1.determinant(), 8); - - EXPECT_ARRAY_NEAR(dm2.matrix(), true_dm2); - EXPECT_EQ(dm2.determinant(), 1); - - EXPECT_ARRAY_NEAR(dm3.matrix(), true_dm3); - EXPECT_EQ(dm3.determinant(), 8); - - EXPECT_ARRAY_NEAR(dm4.matrix(), true_dm4); - EXPECT_EQ(dm4.determinant(), 1); -} - -// ------------------------ - -MAKE_MAIN; diff --git a/test/c++/det_manip/det_manip4.cpp b/test/c++/det_manip/det_manip4.cpp deleted file mode 100644 index 67416a510a..0000000000 --- a/test/c++/det_manip/det_manip4.cpp +++ /dev/null @@ -1,109 +0,0 @@ -// Copyright (c) 2016-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) -// Copyright (c) 2016-2018 Centre national de la recherche scientifique (CNRS) -// Copyright (c) 2018-2021 Simons Foundation -// Copyright (c) 2016 Igor Krivenko -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Igor Krivenko, Olivier Parcollet, Nils Wentzell - -#include - -#include -#include -#include -#include -namespace arrays = nda; -using _vector = nda::vector; -using _matrix = nda::matrix; -using triqs::det_manip::det_manip; - -struct func { - double operator()(int x, int y) const { - if ((x < 0) and (y < 0)) - return 0; - else - return 2 / double(x + y); - } -}; - -det_manip init_dm(int size) { - func f; - std::vector initial_x(size), initial_y(size); - std::iota(initial_x.begin(), initial_x.end(), 1); - std::iota(initial_y.begin(), initial_y.end(), 1); - det_manip dm(f, initial_x, initial_y); - std::cerr << "matrix = " << dm.matrix() << std::endl; - std::cerr << "det = " << dm.determinant() << std::endl; - return dm; -} - -// ----- TESTS ------------------ - -TEST(det_manip, det_manip_refill_0) { - auto dm = init_dm(3); - EXPECT_CLOSE(1.0 / 5400, dm.determinant()); - - EXPECT_NEAR(5400, dm.try_refill(_vector{}, _vector{}), 1e-9); - dm.complete_operation(); - - EXPECT_ARRAY_NEAR(_matrix{}, dm.matrix()); - EXPECT_CLOSE(1, dm.determinant()); - EXPECT_ARRAY_NEAR(_matrix{}, dm.inverse_matrix()); -} - -TEST(det_manip, det_manip_refill_2) { - auto dm = init_dm(3); - EXPECT_CLOSE(1.0 / 5400, dm.determinant()); - - EXPECT_NEAR(300.0, dm.try_refill(_vector{1, 2}, _vector{1, 2}), 1e-10); - dm.complete_operation(); - - EXPECT_ARRAY_NEAR(_matrix{{1.0, 2.0 / 3}, {2.0 / 3, 0.5}}, dm.matrix()); - EXPECT_CLOSE(1.0 / 18, dm.determinant()); - EXPECT_ARRAY_NEAR(_matrix{{9, -12}, {-12, 18}}, dm.inverse_matrix()); -} - -TEST(det_manip, det_manip_refill_4) { - auto dm = init_dm(3); - EXPECT_CLOSE(1.0 / 5400, dm.determinant()); - - EXPECT_NEAR(1.0 / 4900, dm.try_refill(_vector{1, 2, 3, 4}, _vector{1, 2, 3, 4}), 1e-10); - dm.complete_operation(); - - EXPECT_ARRAY_NEAR(_matrix{{1.0, 2.0 / 3, 1.0 / 2, 2.0 / 5}, - {2.0 / 3, 1.0 / 2, 2.0 / 5, 1.0 / 3}, - {1.0 / 2, 2.0 / 5, 1.0 / 3, 2.0 / 7}, - {2.0 / 5, 1.0 / 3, 2.0 / 7, 1.0 / 4}}, - dm.matrix()); - EXPECT_CLOSE(1.0 / 26460000, dm.determinant()); - EXPECT_ARRAY_NEAR(_matrix{{100, -600, 1050, -560}, {-600, 4050, -7560, 4200}, {1050, -7560, 14700, -8400}, {-560, 4200, -8400, 4900}}, - dm.inverse_matrix(), 1e-7); -} - -TEST(det_manip, det_manip_refill_empty) { - auto dm = init_dm(0); - EXPECT_CLOSE(1.0, dm.determinant()); - - EXPECT_NEAR(1.0 / 18, dm.try_refill(_vector{1, 2}, _vector{1, 2}), 1e-10); - dm.complete_operation(); - - EXPECT_ARRAY_NEAR(_matrix{{1.0, 2.0 / 3}, {2.0 / 3, 0.5}}, dm.matrix()); - EXPECT_CLOSE(1.0 / 18, dm.determinant()); - EXPECT_ARRAY_NEAR(_matrix{{9, -12}, {-12, 18}}, dm.inverse_matrix()); -} - -// ------------------------ - -MAKE_MAIN; diff --git a/test/c++/det_manip/det_manip_basics.cpp b/test/c++/det_manip/det_manip_basics.cpp new file mode 100644 index 0000000000..07ac6bfbeb --- /dev/null +++ b/test/c++/det_manip/det_manip_basics.cpp @@ -0,0 +1,135 @@ +// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) +// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS) +// Copyright (c) 2018-2023 Simons Foundation +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You may obtain a copy of the License at +// https://www.gnu.org/licenses/gpl-3.0.txt +// +// Authors: Henri Menke, Olivier Parcollet, Nils Wentzell + +#include +#include +#include + +#include + +struct builder { + double operator()(int i, int j) const { return (i == j ? 1 + j : 0); } +}; + +TEST(TRIQSDetManip, DetManipConstructWithMatrixBuilderAndCapacity) { + using namespace triqs::det_manip; + + // create an empty det_manip object + det_manip dm{builder{}, 10}; + EXPECT_EQ(dm.size(), 0); + EXPECT_EQ(dm.determinant(), 1); + EXPECT_EQ(dm.capacity(), 10); + + // try to reserve memory for 5 elements --> should do nothing + dm.reserve(5); + EXPECT_EQ(dm.capacity(), 10); + + // reserve memory for 50 elements + dm.reserve(50); + EXPECT_EQ(dm.capacity(), 50); + + // clear the data + dm.clear(); + EXPECT_EQ(dm.size(), 0); + EXPECT_EQ(dm.determinant(), 1); + EXPECT_EQ(dm.capacity(), 50); + + // check the data + EXPECT_NO_THROW(dm.regenerate_and_check()); +} + +TEST(TRIQSDetManip, DetManipConstructWithMatrixBuilderAndRanges) { + using namespace triqs::det_manip; + std::vector x_args{0, 1, 2, 3, 4}; + auto y_args = x_args; + auto exp_mat = nda::matrix(5, 5); + nda::for_each(exp_mat.shape(), [&](auto i, auto j) { exp_mat(i, j) = builder{}(i, j); }); + auto exp_inv_mat = nda::inverse(exp_mat); + auto exp_det = nda::determinant(exp_mat); + auto check = [&](auto const &dm, auto exp_size, auto exp_cap) { + EXPECT_EQ(dm.size(), exp_size); + EXPECT_EQ(dm.capacity(), exp_cap); + EXPECT_DOUBLE_EQ(exp_det, dm.determinant()); + EXPECT_ARRAY_NEAR(exp_mat, dm.matrix()); + EXPECT_ARRAY_NEAR(exp_inv_mat, dm.inverse_matrix()); + EXPECT_ARRAY_NEAR(exp_inv_mat, dm.inverse_matrix_internal_order()); + for (int i = 0; i < exp_size; ++i) { + EXPECT_EQ(dm.get_x(i), x_args[i]); + EXPECT_EQ(dm.get_y(i), y_args[i]); + EXPECT_EQ(dm.get_x_internal_order()[i], x_args[i]); + EXPECT_EQ(dm.get_y_internal_order()[i], y_args[i]); + } + }; + + // create a det_manip object with a matrix builder and argument ranges + det_manip dm{builder{}, x_args, y_args}; + check(dm, 5, 10); + + // try to reserve memory for 5 elements --> should do nothing + dm.reserve(5); + check(dm, 5, 10); + + // reserve memory for 20 elements + dm.reserve(20); + check(dm, 5, 20); + + // check the data + EXPECT_NO_THROW(dm.regenerate_and_check()); + + // print a warning + dm.set_precision_warning(-1); + EXPECT_NO_THROW(dm.regenerate_and_check()); + + // clear the data + dm.clear(); + EXPECT_EQ(dm.size(), 0); + EXPECT_EQ(dm.determinant(), 1); + EXPECT_EQ(dm.capacity(), 20); + EXPECT_NO_THROW(dm.regenerate_and_check()); +} + +TEST(TRIQSDetManip, DetManipHDF5) { + using namespace triqs::det_manip; + std::vector x_args{0, 1, 2, 3, 4}; + auto y_args = x_args; + + // create a det_manip object with a matrix builder and argument ranges + det_manip dm{builder{}, x_args, y_args}; + + // write to HDF5 + h5::file file("det_manip.h5", 'w'); + h5::write(file, "det", dm); + + // read from HDF5 + det_manip dm_r{builder{}, 1}; + h5::read(file, "det", dm_r); + + // check the read object + EXPECT_EQ(dm_r.size(), dm.size()); + EXPECT_EQ(dm_r.capacity(), dm.capacity()); + EXPECT_EQ(dm_r.determinant(), dm.determinant()); + EXPECT_ARRAY_EQ(dm_r.matrix(), dm.matrix()); + EXPECT_ARRAY_EQ(dm_r.inverse_matrix(), dm.inverse_matrix()); + EXPECT_ARRAY_EQ(dm_r.inverse_matrix_internal_order(), dm.inverse_matrix_internal_order()); + EXPECT_EQ(dm_r.get_x_internal_order(), dm.get_x_internal_order()); + EXPECT_EQ(dm_r.get_y_internal_order(), dm.get_y_internal_order()); + EXPECT_NO_THROW(dm_r.regenerate_and_check()); +} + +MAKE_MAIN diff --git a/test/c++/det_manip/det_manip_constructors.cpp b/test/c++/det_manip/det_manip_constructors.cpp deleted file mode 100644 index aeec7e3d14..0000000000 --- a/test/c++/det_manip/det_manip_constructors.cpp +++ /dev/null @@ -1,47 +0,0 @@ -// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) -// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS) -// Copyright (c) 2018 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Laura Messio, Olivier Parcollet, Nils Wentzell - -#include - -struct fun { - - typedef double result_type; - typedef double argument_type; - - //gives the coefficients of the matrix (function F of the documentation) - double operator()(double x, double y) const { return (x - y + 1) * (x - y); } -}; - -int main() { - try { - fun f; - int init_size = 10; - std::vector initial_x{1, 2, 2.5}, initial_y{3, 4, 9}; - - //creation of an empty class det_manip - triqs::det_manip::det_manip D1(f, init_size); - - //creation of a class det_manip with a 3 by 3 matrix - triqs::det_manip::det_manip D2(f, initial_x, initial_y); - - //the initial matrix: - std::cout << std::endl << "After construction: D.matrix()=" << D1.matrix() << std::endl; - std::cout << std::endl << "After construction: D.matrix()=" << D2.matrix() << std::endl; - } catch (std::exception const &e) { std::cout << "error " << e.what() << std::endl; } -} diff --git a/test/c++/det_manip/det_manip_insert_remove.cpp b/test/c++/det_manip/det_manip_insert_remove.cpp new file mode 100644 index 0000000000..3bb6d74537 --- /dev/null +++ b/test/c++/det_manip/det_manip_insert_remove.cpp @@ -0,0 +1,85 @@ +// Copyright (c) 2016-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) +// Copyright (c) 2016-2018 Centre national de la recherche scientifique (CNRS) +// Copyright (c) 2018-2021 Simons Foundation +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You may obtain a copy of the License at +// https://www.gnu.org/licenses/gpl-3.0.txt +// +// Authors: Olivier Parcollet, Nils Wentzell + +#include "./det_manip_test_utils.hpp" + +#include +#include + +#include + +#include + +// Construct a det_manip object. +auto init_dm() { + std::vector initial_x{-2, 2}, initial_y{-5, 3}; + return triqs::det_manip::det_manip{builder3{}, initial_x, initial_y}; +} + +// Remove the second row and column. +auto remove_second_row_col(auto dm) { + dm.try_remove(1, 1); + dm.complete_operation(); + return dm; +} + +// Insert a row and column at position 1. +auto insert_second_row_col(auto dm) { + dm.try_insert(1, 1, 6, 4); + dm.complete_operation(); + return dm; +} + +// Remove the first row and column. +auto remove_first_row_col(auto dm) { + dm.try_remove(0, 0); + dm.complete_operation(); + return dm; +} + +TEST(TRIQSDetManip, InsertAndRemoveSpecificRowsAndColumns) { + std::vector x_args{-2, 2}, y_args{-5, 3}; + + // det_manip objects + auto dm1 = triqs::det_manip::det_manip{builder3{}, x_args, y_args}; + auto dm2 = remove_second_row_col(dm1); + auto dm3 = insert_second_row_col(dm2); + auto dm4 = remove_first_row_col(dm2); + + // expected matrices + nda::matrix F1 = {{1, 1}, {-3, 5}}; + nda::matrix F2 = {{1}}; + nda::matrix F3 = {{1, 2}, {1, 10}}; + nda::matrix F4 = {}; + + // check results + EXPECT_ARRAY_NEAR(dm1.matrix(), F1); + EXPECT_EQ(dm1.determinant(), 8); + + EXPECT_ARRAY_NEAR(dm2.matrix(), F2); + EXPECT_EQ(dm2.determinant(), 1); + + EXPECT_ARRAY_NEAR(dm3.matrix(), F3); + EXPECT_EQ(dm3.determinant(), 8); + + EXPECT_ARRAY_NEAR(dm4.matrix(), F4); + EXPECT_EQ(dm4.determinant(), 1); +} + +MAKE_MAIN; diff --git a/test/c++/det_manip/det_manip_test_utils.hpp b/test/c++/det_manip/det_manip_test_utils.hpp index b857b9d7cf..91be6197cc 100644 --- a/test/c++/det_manip/det_manip_test_utils.hpp +++ b/test/c++/det_manip/det_manip_test_utils.hpp @@ -50,4 +50,14 @@ struct builder2 { else return 2 / double(x + y); } -}; \ No newline at end of file +}; + +// Matrix builder #3. +struct builder3 { + double operator()(int x, int y) const { + if ((x < 0) and (y < 0)) + return 1; + else + return x + y; + } +}; diff --git a/test/c++/det_manip/det_manipc1.cpp b/test/c++/det_manip/det_manipc1.cpp deleted file mode 100644 index 73993cbdb2..0000000000 --- a/test/c++/det_manip/det_manipc1.cpp +++ /dev/null @@ -1,156 +0,0 @@ -// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) -// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS) -// Copyright (c) 2018-2023 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Henri Menke, Olivier Parcollet, Nils Wentzell - -#include -#include -#include -#include "./old_test_tool.hpp" - -#include -#include - -struct fun { - - typedef std::complex result_type; - typedef double argument_type; - - /*double f(double eps, double tau, double beta) const { - bool s = (tau>0); - tau = (s ? tau : beta + tau); - return (s ? 1 : -1 ) * exp( - eps* (tau)) / (1 + exp(- beta *eps)); - }*/ - - std::complex operator()(double x, double y) const { - const double pi = acos(-1); - const double beta = 10.0; - const double epsi = 0.1; - double tau = x - y; - bool s = (tau > 0); - tau = (s ? tau : beta + tau); - double r = epsi + tau / beta * (1 - 2 * epsi); - return -2 * (pi / beta) / std::sin(pi * r); - - //return - 0.5*(f(0.12, x-y, beta) + f(0.43,x-y,beta) ); //+ f(-0.7,x-y,beta)); - } -}; - -template void assert_close(T1 const &A, T2 const &B, double precision) { - std::cerr << " assert_close " << A << " " << B << " " << std::abs(A - B) << " " << precision << std::endl; - if (std::abs(A - B) > precision) TRIQS_RUNTIME_ERROR << "assert_close error : " << A << "\n" << B; -} -const double PRECISION = 1.e-6; - -struct test { - - fun f; - triqs::det_manip::det_manip D; - std::complex det_old, detratio; - - test() : f(), D(f, 100) {} - - void check() { -#ifndef PRINT_ALL - std::cerr << "det = " << D.determinant() << " == " << std::complex(determinant(D.matrix())) << std::endl; -#else - std::cerr << "det = " << D.determinant() << " == " << std::complex(determinant(D.matrix())) << std::endl - << D.inverse_matrix() << D.matrix() << nda::matrix>(inverse(D.matrix())) << std::endl; - std::cerr << "det_old = " << det_old << "detratio = " << detratio << " determin " << D.determinant() << std::endl; -#endif - assert_close(D.determinant(), long(1) / determinant(D.inverse_matrix()), PRECISION); - nda::assert_all_close(inverse(D.matrix()), D.inverse_matrix(), PRECISION, true); - assert_close(det_old * detratio, D.determinant(), PRECISION); - } - - void run() { - triqs::mc_tools::random_generator RNG("mt19937", 23432); - for (size_t i = 0; i < 100; ++i) { - std::cerr << " ------------------------------------------------" << std::endl; - std::cerr << " i = " << i << " size = " << D.size() << std::endl; - // choose a move - size_t s = D.size(); - det_old = D.determinant(); - detratio = 1; - double x, y; - bool do_something = true; - - switch (RNG((i > 10 ? 4 : 1))) { - case 0: - std::cerr << " insert1" << std::endl; - x = RNG(10.0), y = RNG(10.0); - std::cerr << " x,y = " << x << " " << y << std::endl; - detratio = D.try_insert(RNG(s), RNG(s), x, y); - break; - case 1: - std::cerr << " remove1" << std::endl; - if (s > 0) detratio = D.try_remove(RNG(s), RNG(s)); - break; - case 2: { - long k = 2 + RNG(6); - std::cerr << " Insert" << k << std::endl; - std::vector x_vec(k), y_vec(k); - std::vector i_vec(k), j_vec(k); - for (auto m : itertools::range(k)) { - x_vec.at(m) = RNG(10.0); - y_vec.at(m) = RNG(10.0); - i_vec.at(m) = RNG(s + m); - j_vec.at(m) = RNG(s + m); - } - std::sort(i_vec.begin(), i_vec.end()); - std::sort(j_vec.begin(), j_vec.end()); - if (std::unique(i_vec.begin(), i_vec.end()) == i_vec.end() && std::unique(j_vec.begin(), j_vec.end()) == j_vec.end()) { - detratio = D.try_insert_k(i_vec, j_vec, x_vec, y_vec); - } else - do_something = false; - break; - } - case 3: { - long k = 2 + RNG(6); - std::cerr << " Remove" << k << std::endl; - if (D.size() >= k) { - std::vector i_vec(k), j_vec(k); - for (auto m : itertools::range(k)) { - i_vec.at(m) = RNG(s); - j_vec.at(m) = RNG(s); - } - std::sort(i_vec.begin(), i_vec.end()); - std::sort(j_vec.begin(), j_vec.end()); - if (std::unique(i_vec.begin(), i_vec.end()) == i_vec.end() && std::unique(j_vec.begin(), j_vec.end()) == j_vec.end()) { - detratio = D.try_remove_k(i_vec, j_vec); - } else - do_something = false; - } - break; - } - default: TRIQS_RUNTIME_ERROR << " TEST INTERNAL ERROR"; - }; - - if (do_something) { - if (std::abs(detratio * det_old) > PRECISION) { - D.complete_operation(); - if (D.size() > 0) check(); - } else { - std::cerr << " reject since new det is = " << std::abs(detratio * det_old) << std::endl; - D.reject_last_try(); - } - } - } - } -}; - -int main(int, char **) { test().run(); } diff --git a/test/c++/det_manip/det_manipc2.cpp b/test/c++/det_manip/det_manipc2.cpp deleted file mode 100644 index 5f65c43bb1..0000000000 --- a/test/c++/det_manip/det_manipc2.cpp +++ /dev/null @@ -1,158 +0,0 @@ -// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) -// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS) -// Copyright (c) 2018-2023 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Henri Menke, Olivier Parcollet, Nils Wentzell - -#include -#include -#include - -#include "./old_test_tool.hpp" - -#include -#include - -struct fun { - - typedef std::complex result_type; - typedef double argument_type; - - /*double f(double eps, double tau, double beta) const { - bool s = (tau>0); - tau = (s ? tau : beta + tau); - return (s ? 1 : -1 ) * exp( - eps* (tau)) / (1 + exp(- beta *eps)); - }*/ - - std::complex operator()(double x, double y) const { - const double pi = acos(-1); - const double beta = 10.0; - const double epsi = 0.1; - double tau = x - y; - bool s = (tau > 0); - tau = (s ? tau : beta + tau); - double r = epsi + tau / beta * (1 - 2 * epsi); - std::complex I(0, 1); - return -2 * I * (pi / beta) / std::sin(pi * r); - - //return - 0.5*(f(0.12, x-y, beta) + f(0.43,x-y,beta) ); //+ f(-0.7,x-y,beta)); - } -}; - -template void assert_close(T1 const &A, T2 const &B, double precision) { - std::cerr << " assert_close " << A << " " << B << " " << std::abs(A - B) << " " << precision << std::endl; - if (std::abs(A - B) > precision) TRIQS_RUNTIME_ERROR << "assert_close error : " << A << "\n" << B; -} -const double PRECISION = 1.e-6; - -struct test { - - fun f; - triqs::det_manip::det_manip D; - std::complex det_old, detratio; - - test() : f(), D(f, 100) {} - - void check() { -#ifndef PRINT_ALL - std::cerr << "det = " << D.determinant() << " == " << std::complex(determinant(D.matrix())) << std::endl; -#else - std::cerr << "det = " << D.determinant() << " == " << std::complex(determinant(D.matrix())) << std::endl - << D.inverse_matrix() << D.matrix() << nda::matrix>(inverse(D.matrix())) << std::endl; - std::cerr << "det_old = " << det_old << "detratio = " << detratio << " determin " << D.determinant() << std::endl; -#endif - assert_close(D.determinant(), long(1) / determinant(D.inverse_matrix()), PRECISION); - nda::assert_all_close(inverse(D.matrix()), D.inverse_matrix(), PRECISION, true); - assert_close(det_old * detratio, D.determinant(), PRECISION); - } - - void run() { - triqs::mc_tools::random_generator RNG("mt19937", 23432); - for (size_t i = 0; i < 100; ++i) { - std::cerr << " ------------------------------------------------" << std::endl; - std::cerr << " i = " << i << " size = " << D.size() << std::endl; - // choose a move - size_t s = D.size(); - det_old = D.determinant(); - detratio = 1; - double x, y; - bool do_something = true; - - switch (RNG((i > 10 ? 4 : 1))) { - case 0: - std::cerr << " insert1" << std::endl; - x = RNG(10.0), y = RNG(10.0); - std::cerr << " x,y = " << x << " " << y << std::endl; - detratio = D.try_insert(RNG(s), RNG(s), x, y); - break; - case 1: - std::cerr << " remove1" << std::endl; - if (s > 0) detratio = D.try_remove(RNG(s), RNG(s)); - break; - case 2: { - long k = 2 + RNG(6); - std::cerr << " Insert" << k << std::endl; - std::vector x_vec(k), y_vec(k); - std::vector i_vec(k), j_vec(k); - for (auto m : itertools::range(k)) { - x_vec.at(m) = RNG(10.0); - y_vec.at(m) = RNG(10.0); - i_vec.at(m) = RNG(s + m); - j_vec.at(m) = RNG(s + m); - } - std::sort(i_vec.begin(), i_vec.end()); - std::sort(j_vec.begin(), j_vec.end()); - if (std::unique(i_vec.begin(), i_vec.end()) == i_vec.end() && std::unique(j_vec.begin(), j_vec.end()) == j_vec.end()) { - detratio = D.try_insert_k(i_vec, j_vec, x_vec, y_vec); - } else - do_something = false; - break; - } - case 3: { - long k = 2 + RNG(6); - std::cerr << " Remove" << k << std::endl; - if (D.size() >= k) { - std::vector i_vec(k), j_vec(k); - for (auto m : itertools::range(k)) { - i_vec.at(m) = RNG(s); - j_vec.at(m) = RNG(s); - } - std::sort(i_vec.begin(), i_vec.end()); - std::sort(j_vec.begin(), j_vec.end()); - if (std::unique(i_vec.begin(), i_vec.end()) == i_vec.end() && std::unique(j_vec.begin(), j_vec.end()) == j_vec.end()) { - detratio = D.try_remove_k(i_vec, j_vec); - } else - do_something = false; - } - break; - } - default: TRIQS_RUNTIME_ERROR << " TEST INTERNAL ERROR"; - }; - - if (do_something) { - if (std::abs(detratio * det_old) > PRECISION) { - D.complete_operation(); - if (D.size() > 0) check(); - } else { - std::cerr << " reject since new det is = " << std::abs(detratio * det_old) << std::endl; - D.reject_last_try(); - } - } - } - } -}; - -int main(int, char **) { test().run(); } diff --git a/test/c++/det_manip/detmanip.cpp b/test/c++/det_manip/detmanip.cpp deleted file mode 100644 index c2697db3dc..0000000000 --- a/test/c++/det_manip/detmanip.cpp +++ /dev/null @@ -1,98 +0,0 @@ -// Copyright (c) 2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA) -// Copyright (c) 2018 Centre national de la recherche scientifique (CNRS) -// Copyright (c) 2018-2023 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Olivier Parcollet, Nils Wentzell - -#include -#include -#include -#include -//#include -#include -#include "./old_test_tool.hpp" - -struct fun { - - using result_type = double; - using argument_type = double; - - double operator()(double x, double y) const { - const double pi = acos(-1); - const double beta = 10.0; - const double epsi = 0.1; - double tau = x - y; - bool s = (tau > 0); - tau = (s ? tau : beta + tau); - double r = epsi + tau / beta * (1 - 2 * epsi); - return -2 * (pi / beta) / std::sin(pi * r); - } -}; - -using d_t = triqs::det_manip::det_manip; -const double precision = 1.e-6; -const int seed = 23432; - -TEST(DetManip, ChangeRowCol) { - - for (int N = 1; N < 9; N++) { - std::cerr << "================================\n N = " << N << std::endl; - - for (int i0 = 0; i0 < N; i0++) - for (int j0 = 0; j0 < N; j0++) { - - std::cerr << "------------------------------\n i0 = " << i0 << "j0 = " << j0 << std::endl; - - std::mt19937 gen(seed); - std::uniform_real_distribution<> dis(0.0, 10.0); - - std::vector X(N), Y(N); - for (int n = 0; n < N; ++n) { - X[n] = dis(gen); - Y[n] = dis(gen); - } - - auto d = d_t{fun{}, X, Y}; - - // result - auto x = dis(gen), y = dis(gen); - X[i0] = x; - Y[j0] = y; - auto d2 = d_t{fun{}, X, Y}; - auto det1 = d.determinant(); - if (std::abs(det1) < 1.e-5) { - std::cerr << " Det Too small: not completed"; - continue; - } - - // ops - auto detratio = d.try_change_col_row(i0, j0, x, y); - - d.complete_operation(); - - auto det = d.determinant(); - auto det2 = d2.determinant(); - auto det_check = 1 / determinant(d.inverse_matrix()); - - if (std::abs(detratio - det / det1) > precision) TRIQS_RUNTIME_ERROR << "detratio incorrect : " << detratio << " " << det / det1; - if (std::abs(det - det2) > precision) TRIQS_RUNTIME_ERROR << "Det != d2 : " << det << " " << det2; - if (std::abs(det - det_check) > precision) TRIQS_RUNTIME_ERROR << "Det != det_check : " << det << " " << det_check; - nda::assert_all_close(nda::matrix{inverse(d.matrix())}, d.inverse_matrix(), precision, true); - } - } -} - -MAKE_MAIN; diff --git a/test/c++/det_manip/old_test_tool.hpp b/test/c++/det_manip/old_test_tool.hpp deleted file mode 100644 index 496d16f73e..0000000000 --- a/test/c++/det_manip/old_test_tool.hpp +++ /dev/null @@ -1,56 +0,0 @@ -// Copyright (c) 2020 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Olivier Parcollet - -#include - -namespace nda { - - template inline double assert_abs(T z) { return std::abs(z); } - - template - void assert_all_close(ArrayType1 const &A, ArrayType2 const &B, double precision, bool relative = false) { - typedef typename ArrayType1::value_type F; - auto Abs = map(std::function(assert_abs)); - auto r = max_element(Abs(A - B)); - auto r2 = max_element(Abs(A) + Abs(B)); - if (r > (relative ? precision * r2 : precision)) - TRIQS_RUNTIME_ERROR << "assert_all_close error : \n\n" - << ".. A = " << A << "\n\n" - << ".. B= " << B << "\n\n" - << ".. Residue is r = " << r; - } - - template void assert_equal(T1 const &x1, T2 const &x2) { - if (!(x1 == x2)) TRIQS_RUNTIME_ERROR << "Test assert_equal : failed"; - } - - void assert_close(double x1, double x2) { - if (std::abs(x1 - x2) > 1.e-13) TRIQS_RUNTIME_ERROR << "Test assert_equal : failed"; - } - void assert_close(std::complex x1, std::complex x2) { - if (std::abs(x1 - x2) > 1.e-13) TRIQS_RUNTIME_ERROR << "Test assert_equal : failed"; - } - - // Put here a macro, to print the expression, the line, etc... - void assert_is_true(bool b) { - if (!b) TRIQS_RUNTIME_ERROR << "Test assert_is_true : failed"; - } - - void assert_is_false(bool b) { - if (b) TRIQS_RUNTIME_ERROR << "Test assert_is_false : failed"; - } -} // namespace nda diff --git a/test/c++/det_manip/swap_line_col.cpp b/test/c++/det_manip/swap_line_col.cpp deleted file mode 100644 index 8e14f3c0b4..0000000000 --- a/test/c++/det_manip/swap_line_col.cpp +++ /dev/null @@ -1,119 +0,0 @@ -// Copyright (c) 2022 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Olivier Parcollet - -#include -#include -#include -#include -#include -#include "./old_test_tool.hpp" -#include - -struct fun { - - using result_type = double; - using argument_type = double; - - double operator()(double x, double y) const { - const double pi = acos(-1); - const double beta = 10.0; - const double epsi = 0.1; - double tau = x - y; - bool s = (tau > 0); - tau = (s ? tau : beta + tau); - double r = epsi + tau / beta * (1 - 2 * epsi); - return -2 * (pi / beta) / std::sin(pi * r); - } -}; - -using d_t = triqs::det_manip::det_manip; -// const double precision = 1.e-6; -const int seed = 23432; - -TEST(DetManip, SwapRow) { - - for (int N = 1; N < 9; N++) { - std::cerr << "================================\n N = " << N << std::endl; - - for (int i0 = 0; i0 < N; i0++) - for (int j0 = 0; j0 < N; j0++) { - - std::cerr << "------------------------------\n i0 = " << i0 << "j0 = " << j0 << std::endl; - - std::mt19937 gen(seed); - std::uniform_real_distribution<> dis(0.0, 10.0); - - std::vector X(N), Y(N); - for (int n = 0; n < N; ++n) { - X[n] = dis(gen); - Y[n] = dis(gen); - } - - auto d = d_t{fun{}, X, Y}; - - // result - std::swap(X[i0], X[j0]); - auto d2 = d_t{fun{}, X, Y}; - - // ops - d.swap_row(i0, j0); - - EXPECT_ARRAY_NEAR(d.matrix(), d2.matrix(), 1.e-13); - EXPECT_NEAR(d.determinant(), d2.determinant(), 1.e-13); - - //nda::assert_all_close(d.matrix(), d2.matrix(), precision, true); - } - } -} - -TEST(DetManip, SwapCol) { - - for (int N = 1; N < 9; N++) { - std::cerr << "================================\n N = " << N << std::endl; - - for (int i0 = 0; i0 < N; i0++) - for (int j0 = 0; j0 < N; j0++) { - - std::cerr << "------------------------------\n i0 = " << i0 << "j0 = " << j0 << std::endl; - - std::mt19937 gen(seed); - std::uniform_real_distribution<> dis(0.0, 10.0); - - std::vector X(N), Y(N); - for (int n = 0; n < N; ++n) { - X[n] = dis(gen); - Y[n] = dis(gen); - } - - auto d = d_t{fun{}, X, Y}; - - // result - std::swap(Y[i0], Y[j0]); - auto d2 = d_t{fun{}, X, Y}; - - // ops - d.swap_col(i0, j0); - - EXPECT_ARRAY_NEAR(d.matrix(), d2.matrix(), 1.e-13); - EXPECT_NEAR(d.determinant(), d2.determinant(), 1.e-13); - - //nda::assert_all_close(d.matrix(), d2.matrix(), precision, true); - } - } -} - -MAKE_MAIN;