rotation funcs

Author: maki49

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_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);
+		}
+
+	}
+
+}