Merge branch 'loop3-symmetry-new' into filter_exx

Author: PeizeLin

include/RI/global/Tensor.h

@@ -45,6 +45,7 @@ class Tensor
 	inline T& operator() (const std::size_t i0, const std::size_t i1, const std::size_t i2, const std::size_t i3) const;
 
 	Tensor transpose() const;
+	Tensor dagger() 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

@@ -178,6 +178,17 @@ Tensor<T> Tensor<T>::transpose() const
 	return t;
 }
 
+template<typename T>
+Tensor<T> Tensor<T>::dagger() const
+{
+	assert(this->shape.size() == 2);
+	Tensor<T> t({ this->shape[1], this->shape[0] });
+	for (std::size_t i0 = 0; i0 < this->shape[0]; ++i0)
+		for (std::size_t i1 = 0; i1 < this->shape[1]; ++i1)
+			t(i1, i0) = std::conj((*this)(i0, i1));
+	return t;
+}
+
 template<typename T>
 Global_Func::To_Real_t<T> Tensor<T>::norm(const double p) const
 {

include/RI/symmetry/Symmetry_Filter.h

@@ -0,0 +1,68 @@
+#include <array>
+#include <map>
+#include <set>
+#include <tuple>
+#define NO_SEC_RETURN_TRUE if(this->irreducible_sector_.empty()) return true;
+#include "../global/Array_Operator.h"
+namespace RI
+{
+	using namespace Array_Operator;
+	template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
+	class Symmetry_Filter
+	{
+		using TC = std::array<Tcell, Ndim>;
+		using TAC = std::pair<TA, TC>;
+
+		using TIJ = std::pair<TA, TA>;
+		using TIJR = std::pair<TIJ, TC>;
+		using Tsec = std::map<TIJ, std::set<TC>>;
+	public:
+		Symmetry_Filter(const TC& period_in, const Tsec& irsec)
+			:period(period_in), irreducible_sector_(irsec) {}
+		bool in_irreducible_sector(const TA& Aa, const TAC& Ab) const
+		{
+			NO_SEC_RETURN_TRUE;
+			const TIJ& ap = { Aa, Ab.first };
+			if (irreducible_sector_.find(ap) != irreducible_sector_.end())
+				if (irreducible_sector_.at(ap).find(Ab.second % this->period) != irreducible_sector_.at(ap).end())
+					return true;
+			return false;
+		}
+		bool in_irreducible_sector(const TAC& Aa, const TAC& Ab) const
+		{
+			NO_SEC_RETURN_TRUE;
+			const TC dR = (Ab.second - Aa.second) % this->period;
+			const std::pair<TA, TA> ap = { Aa.first, Ab.first };
+			if (irreducible_sector_.find(ap) != irreducible_sector_.end())
+				if (irreducible_sector_.at(ap).find(dR) != irreducible_sector_.at(ap).end())
+					return true;
+			return false;
+		}
+		bool is_I_in_irreducible_sector(const TA& Aa) const
+		{
+			NO_SEC_RETURN_TRUE;
+			for (const auto& apRs : irreducible_sector_)
+				if (apRs.first.first == Aa)return true;
+			return false;
+		}
+		bool is_J_in_irreducible_sector(const TA& Ab) const
+		{
+			NO_SEC_RETURN_TRUE;
+			for (const auto& apRs : irreducible_sector_)
+				if (apRs.first.second == Ab)return true;
+			return false;
+		}
+		TIJR get_IJR(const TA& I, const TAC& J) const
+		{
+			return { {I,J.first}, J.second % this->period };
+		}
+		TIJR get_IJR(const TAC& I, const TAC& J) const
+		{
+			return { {I.first,J.first}, (J.second - I.second) % this->period };
+		}
+	private:
+		const Tsec& irreducible_sector_;
+		const TC& period;
+	};
+
+}

include/RI/symmetry/Symmetry_Rotation.h

@@ -0,0 +1,55 @@
+#include <RI/global/Tensor.h>
+namespace RI
+{
+	namespace Sym
+	{
+		template <typename T>
+		inline void T1_HR(T* TA, const T* A, const Tensor<T>& T1, const int& n2)
+		{
+			// C' = T1^\dagger * C
+			const int& n1 = T1.shape[0];
+			Blas_Interface::gemm('C', 'N', n1, n2, n1,
+				T(1), T1.ptr(), n1, A, n2, T(0), TA, n2);
+			// zgemm_(&notrans, &dagger, &n12, &nabf, &nabf,
+			// 	&alpha, A, &n12, T1.ptr(), &nabf, &beta, TA, &n12);
+		}
+
+		template <typename T>
+		inline void T1_HR_T2(T* TAT, const T* A, const Tensor<T>& T1, const Tensor<T>& T2)
+		{
+			// H' = T1^\dagger * H * T2
+			const int& n2 = T2.shape[0], & n1 = T1.shape[0];
+			const RI::Shape_Vector& shape = { static_cast<size_t>(n1),static_cast<size_t>(n2) };
+			RI::Tensor<T> AT2(shape);
+			Blas_Interface::gemm('N', 'N', n1, n2, n2,
+				T(1), A, n2, T2.ptr(), n2, T(0), AT2.ptr(), n2);
+			Blas_Interface::gemm('C', 'N', n1, n2, n1,
+				T(1), T1.ptr(), n1, AT2.ptr(), n2, T(0), TAT, n2);
+			// col-major version
+			// zgemm_(&notrans, &notrans, &n2, &n1, &n2,
+			// 	&alpha, T2.ptr(), &n2, A, &n2, &beta, AT2.ptr(), &n2);
+			// zgemm_(&notrans, &dagger, &n2, &n1, &n1,
+			// 	&alpha, AT2.ptr(), &n2, T1.ptr(), &n1, &beta, TAT, &n2);
+		}
+
+		template<typename T>
+		inline void T1_DR_T2(T* TAT, const T* A, const Tensor<T>& T1, const Tensor<T>& T2)
+		{
+			// D' = T1^T * D * T2^*  = T1^T * [T2^\dagger * D^T]^T
+			const int& n2 = T2.shape[0], & n1 = T1.shape[0];
+			const RI::Shape_Vector& shape = { static_cast<size_t>(n1),static_cast<size_t>(n2) };
+			RI::Tensor<T> AT2(shape);
+			BlasConnector::gemm('C', 'T', n2, n1, n2,
+				T(1), T2.ptr(), n2, A, n2, T(0), AT2.ptr(), n1);
+			BlasConnector::gemm('T', 'T', n1, n2, n1,
+				T(1), T1.ptr(), n1, AT2.ptr(), n1, T(0), TAT, n2);
+			// col-major version
+			// zgemm_(&transpose, &dagger, &nw1, &nw2, &nw2,
+			// 	&alpha, A, &nw2, T2.ptr(), &nw2, &beta, AT2.ptr(), &nw1);
+			// zgemm_(&transpose, &transpose, &nw2, &nw1, &nw1,
+			// 	&alpha, AT2.ptr(), &nw1, T1.ptr(), &nw1, &beta, TAT, &nw2);
+		}
+
+	}
+
+}