Implement CVC algorithm and enhance LR class

include/RI/distribute/Distribute_Equally.h

@@ -56,6 +56,27 @@ namespace Distribute_Equally
 		const std::array<Tcell,Ndim> &period,
 		const std::size_t num_index,
 		const bool flag_task_repeatable);
+
+	template<typename Tindex>
+	extern void distribute_atom_and_k_pair(
+		const MPI_Comm &mpi_comm,
+		const std::size_t nat,
+		const std::size_t nk,
+		std::vector<Tindex> &list_I,
+		std::vector<Tindex> &list_J,
+		std::vector<Tindex> &list_k1_index,
+		std::vector<Tindex> &list_k2_index,
+		const bool flag_task_repeatable);
+
+	template<typename Tindex>
+	extern void distribute_atom_pair_and_k(
+		const MPI_Comm &mpi_comm,
+		const std::size_t nat,
+		const std::size_t nk,
+		std::vector<Tindex> &list_I,
+		std::vector<Tindex> &list_J,
+		std::vector<Tindex> &list_k_index,
+		const bool flag_task_repeatable);
 }
 
 }

include/RI/distribute/Distribute_Equally.hpp

@@ -126,6 +126,136 @@ namespace Distribute_Equally
 				period);
 		return atoms_split_list;
 	}
+
+	// 均分{atomI,atomJ,k1,k2}
+	template<typename Tindex>
+	void distribute_atom_and_k_pair(
+		const MPI_Comm &mpi_comm,
+		const std::size_t nat,
+		const std::size_t nk,
+		std::vector<Tindex> &list_I,
+		std::vector<Tindex> &list_J,
+		std::vector<Tindex> &list_k1_index,
+		std::vector<Tindex> &list_k2_index,
+		const bool flag_task_repeatable)
+	{
+		// task_sizes的顺序必须从小到大，否则在split中会出现rank_size<group_size，所以先判断nat和nk的大小
+		std::size_t ntaskA, ntaskB;
+		std::vector<Tindex> *A1_ptr, *A2_ptr, *B1_ptr, *B2_ptr;
+		if (nk >= nat)
+		{
+			ntaskA = nat;
+			ntaskB = nk;
+			A1_ptr = &list_I;
+			A2_ptr = &list_J;
+			B1_ptr = &list_k1_index;
+			B2_ptr = &list_k2_index;
+		}
+		else
+		{
+			ntaskA = nk;
+			ntaskB = nat;
+			A1_ptr = &list_k1_index;
+			A2_ptr = &list_k2_index;
+			B1_ptr = &list_I;
+			B2_ptr = &list_J;
+		}
+		const std::vector<std::size_t> task_sizes{ntaskA, ntaskA, ntaskB, ntaskB};
+		const std::vector<std::tuple<MPI_Wrapper::mpi_comm, std::size_t, std::size_t>>
+			comm_color_sizes = Split_Processes::split_all(mpi_comm, task_sizes);
+
+		if(!flag_task_repeatable)
+			if(RI::MPI_Wrapper::mpi_get_rank(std::get<0>(comm_color_sizes.back())()))
+				return;
+
+		std::vector<Tindex> indicesA, indicesB;
+		for(Tindex i=0; i<ntaskA; ++i)
+			indicesA.push_back(i);
+		for(Tindex i=0; i<ntaskB; ++i)
+			indicesB.push_back(i);
+
+		*A1_ptr = Divide_Atoms::divide_atoms(
+			std::get<1>(comm_color_sizes[1]),
+			std::get<2>(comm_color_sizes[1]),
+			indicesA);
+		*A2_ptr = Divide_Atoms::divide_atoms(
+			std::get<1>(comm_color_sizes[2]),
+			std::get<2>(comm_color_sizes[2]),
+			indicesA);
+		*B1_ptr = Divide_Atoms::divide_atoms(
+			std::get<1>(comm_color_sizes[3]),
+			std::get<2>(comm_color_sizes[3]),
+			indicesB);
+		*B2_ptr = Divide_Atoms::divide_atoms(
+			std::get<1>(comm_color_sizes[4]),
+			std::get<2>(comm_color_sizes[4]),
+			indicesB);
+	}
+
+	// 均分{atomI,atomJ,k}
+	template<typename Tindex>
+	void distribute_atom_pair_and_k(
+		const MPI_Comm &mpi_comm,
+		const std::size_t nat,
+		const std::size_t nk,
+		std::vector<Tindex> &list_I,
+		std::vector<Tindex> &list_J,
+		std::vector<Tindex> &list_k_index,
+		const bool flag_task_repeatable)
+	{
+		std::vector<std::size_t> task_sizes;
+		std::vector<Tindex> indices_atom, indices_k;
+		for(Tindex i=0; i<nat; ++i)
+			indices_atom.push_back(i);
+		for(Tindex i=0; i<nk; ++i)
+			indices_k.push_back(i);
+	// task_sizes的顺序必须从小到大，否则在split中会出现rank_size<group_size，所以先判断nat和nk的大小
+		if (nk >= nat)
+		{
+			task_sizes = {nat, nat, nk};
+		}
+		else
+		{
+			task_sizes = {nk, nat, nat};
+		}
+		const std::vector<std::tuple<MPI_Wrapper::mpi_comm, std::size_t, std::size_t>>
+			comm_color_sizes = Split_Processes::split_all(mpi_comm, task_sizes);
+
+		if(!flag_task_repeatable)
+			if(RI::MPI_Wrapper::mpi_get_rank(std::get<0>(comm_color_sizes.back())()))
+				return;
+
+		if (nk >= nat)
+		{
+			list_I = Divide_Atoms::divide_atoms(
+				std::get<1>(comm_color_sizes[1]),
+				std::get<2>(comm_color_sizes[1]),
+				indices_atom);
+			list_J = Divide_Atoms::divide_atoms(
+				std::get<1>(comm_color_sizes[2]),
+				std::get<2>(comm_color_sizes[2]),
+				indices_atom);
+			list_k_index = Divide_Atoms::divide_atoms(
+				std::get<1>(comm_color_sizes[3]),
+				std::get<2>(comm_color_sizes[3]),
+				indices_k);
+		}
+		else
+		{
+			list_k_index = Divide_Atoms::divide_atoms(
+				std::get<1>(comm_color_sizes[1]),
+				std::get<2>(comm_color_sizes[1]),
+				indices_k);
+			list_I = Divide_Atoms::divide_atoms(
+				std::get<1>(comm_color_sizes[2]),
+				std::get<2>(comm_color_sizes[2]),
+				indices_atom);
+			list_J = Divide_Atoms::divide_atoms(
+				std::get<1>(comm_color_sizes[3]),
+				std::get<2>(comm_color_sizes[3]),
+				indices_atom);
+		}
+	}
 }
 
 }

include/RI/global/Array_Operator.h

@@ -13,7 +13,12 @@ namespace RI
 namespace Array_Operator
 {
 	template<typename T, std::size_t N>
-	extern std::array<T,N> operator%(const std::array<T,N> &v1, const std::array<T,N> &v2);
+	typename std::enable_if<std::is_integral<T>::value, std::array<T,N>>::type
+	operator%(const std::array<T,N> &v1, const std::array<T,N> &v2);
+
+    template<typename T, std::size_t N>
+    typename std::enable_if<std::is_floating_point<T>::value, std::array<T,N>>::type
+    operator%(const std::array<T,N> &v1, const std::array<T,N> &v2);
 
 	template<typename T, std::size_t N>
 	extern std::array<T,N> operator+(const std::array<T,N> &v1, const std::array<T,N> &v2);

include/RI/global/Array_Operator.hpp

@@ -13,7 +13,8 @@ namespace RI
 namespace Array_Operator
 {
 	template<typename T, std::size_t N>
-	std::array<T,N> operator%(const std::array<T,N> &v1, const std::array<T,N> &v2)
+	typename std::enable_if<std::is_integral<T>::value, std::array<T,N>>::type
+	operator%(const std::array<T,N> &v1, const std::array<T,N> &v2)
 	{
 		auto mod = [](const T i, const T n){ return (i%n+3*n/2)%n-n/2; };			// [-n/2,n/2]
 	//	auto mod = [](const T i, const T n){ return (i%n+n)%n; };					// [0,n)
@@ -23,6 +24,24 @@ namespace Array_Operator
 			v[i] = mod(v1[i], v2[i]);
 		return v;
 	}
+	template<typename T, std::size_t N>
+	typename std::enable_if<std::is_floating_point<T>::value, std::array<T,N>>::type
+	operator%(const std::array<T,N> &v1, const std::array<T,N> &v2)
+	{
+		constexpr double epsilon = 1e-6;
+		auto mod_f = [&](T x, T period){      // [0, period)
+			T r = std::fmod(x, period);
+			if (std::abs(r) < epsilon) r = 0.0;
+			if (std::abs(r - period) < epsilon) r = 0.0;
+			if (r < 0) r += period;
+
+			return r;
+		};
+		std::array<T,N> v;
+		for(std::size_t i=0; i<N; ++i)
+			v[i] = mod_f(v1[i], v2[i]);
+		return v;
+	}
 
 	template<typename T, std::size_t N>
 	std::array<T,N> operator+(const std::array<T,N> &v1, const std::array<T,N> &v2)

include/RI/global/Global_Func-1.h

@@ -83,6 +83,43 @@ namespace Global_Func
 			return find( ptr->second, keys... );
 	}
 
+	// These functions are designed to replace Global_Func::find,
+	// the find result can be arithmetic type, Tensor, vector, array, or map.
+
+	// Helper: compute return type of find_map with N keys (Map::mapped_type applied N times)
+	template<typename Map, typename... Keys>
+	struct find_map_result;
+
+	template<typename Map, typename Key>
+	struct find_map_result<Map, Key> {
+		using type = typename Map::mapped_type;
+	};
+
+	template<typename Map, typename Key0, typename... Keys>
+	struct find_map_result<Map, Key0, Keys...> {
+		using type = typename find_map_result<typename Map::mapped_type, Keys...>::type;
+	};
+
+	// Base case: 1 key — returns const& to Map::mapped_type (ZERO sentinel if key not found)
+	template<class Map, class Key>
+	static inline const typename Map::mapped_type& find_map(const Map& map, const Key& key)
+	{
+		const auto& it = map.find(key);
+		if (it != map.end()) return it->second;
+		return ZERO<typename Map::mapped_type>;
+	}
+
+	// Recursive case: 2+ keys — peels off key0, recurses on remaining keys
+	template<class Map, class Key0, class... Keys>
+	static inline const typename find_map_result<Map, Key0, Keys...>::type&
+	find_map(const Map& map, const Key0& key0, const Keys&... keys)
+	{
+		const auto it = map.find(key0);
+		if (it != map.end())
+			return find_map(it->second, keys...);
+		return ZERO<typename find_map_result<Map, Key0, Keys...>::type>;
+	}
+
 	// in_set(3, {2,3,5,7})
 	// Peize Lin add 2022.05.26
 	template<typename T>
@@ -121,6 +158,21 @@ namespace Global_Func
 	{
 		return std::vector<T>(v.begin(), v.end());
 	}
+
+	/// @brief sorted unique union of two containers
+	/// @tparam C1,C2  containers supporting .begin()/.end() and iterator-pair construction
+	///               (e.g., std::vector, std::deque, std::list)
+	template<typename C1, typename C2>
+	static auto set_union(const C1& c1, const C2& c2)
+		-> C1
+	{
+		static_assert(
+			std::is_same<typename C1::value_type, typename C2::value_type>::value,
+			"set_union: both containers must have the same value_type");
+		std::set<typename C1::value_type> s(c1.begin(), c1.end());
+		s.insert(c2.begin(), c2.end());
+		return C1(s.begin(), s.end());
+	}
 }
 
 }

include/RI/global/Global_Func-2.h

@@ -49,6 +49,21 @@ namespace Global_Func
 		typename std::enable_if<!Global_Func::is_complex<Tout>::value,int>::type =0>
 	Tout convert(const Tin &t)
 	{ return t.real(); }
+
+	template<
+		typename T,
+		typename std::enable_if<!Global_Func::is_complex<T>::value,int>::type =0>
+	inline T get_conj(const T& x)
+	{
+		return x;
+	}
+	template<
+		typename T,
+		typename std::enable_if< Global_Func::is_complex<T>::value,int>::type =0>
+	inline T get_conj(const T& x)
+	{
+		return std::conj(x);
+	}
 }
 
 }

include/RI/global/Shape_Vector.h

@@ -29,7 +29,12 @@ class Shape_Vector
 		for(auto ptr_in=v_in.begin(); ptr_in<v_in.end(); )
 			*(ptr_this++) = *(ptr_in++);
 	}
-
+	Shape_Vector(const std::vector<std::size_t>& v_in)
+	:size_(v_in.size())
+	{
+		assert(v_in.size() <= sizeof(v) / sizeof(*v));
+		for (std::size_t i = 0;i < size_;++i) this->v[i] = v_in[i];
+	}
 	const std::size_t* begin() const noexcept { return this->v; }
 	const std::size_t* end() const noexcept { return this->v+size_; }
 	std::size_t size() const noexcept { return size_; }

include/RI/global/Tensor.h

@@ -45,6 +45,7 @@ class Tensor
 
 	Tensor transpose() const;
 	Tensor dagger() const;
+	Tensor conjugate() const;
 
 	// ||d||_p = (|d_1|^p+|d_2|^p+...)^{1/p}
 	// if(p==std::numeric_limits<double>::max())    ||d||_max = max_i |d_i|

include/RI/global/Tensor.hpp

@@ -13,6 +13,7 @@
 #include <vector>
 #include <functional>
 #include <cassert>
+#include <iostream>
 #include <limits>
 
 namespace RI
@@ -95,10 +96,18 @@ template<typename T1, typename T2>
 bool same_shape (const Tensor<T1> &t1, const Tensor<T2> &t2)
 {
 	if(t1.shape.size() != t2.shape.size())
+	{
+		std::cerr << "same_shape: ndim mismatch ("
+		          << t1.shape.size() << " vs " << t2.shape.size() << ")" << std::endl;
 		return false;
+	}
 	for(std::size_t ishape=0; ishape<t1.shape.size(); ++ishape)
 		if(t1.shape[ishape] != t2.shape[ishape])
+		{
+			std::cerr << "same_shape: dim[" << ishape << "] mismatch ("
+			          << t1.shape[ishape] << " vs " << t2.shape[ishape] << ")" << std::endl;
 			return false;
+		}
 	return true;
 }
 
@@ -182,6 +191,15 @@ Tensor<T> Tensor<T>::dagger() const
 	return t;
 }
 
+template<typename T>
+Tensor<T> Tensor<T>::conjugate() const
+{
+	Tensor<T> t(this->shape);
+	for (std::size_t i = 0; i < this->data->size(); ++i)
+        (*t.data)[i] = Global_Func::get_conj((*this->data)[i]);
+	return t;
+}
+
 template<typename T>
 Global_Func::To_Real_t<T> Tensor<T>::norm(const double p) const
 {

include/RI/global/Tensor_Multiply-32.hpp

@@ -79,7 +79,6 @@ namespace Tensor_Multiply
 			Tdata(0.0), Txy.ptr());
 		return Txy;
 	}
-
 }
 
 }