Squashed commit of the following:

commit bd40883
Author: maki49 <[email protected]>
Date:   Tue Sep 30 17:49:30 2025 +0800

    remove timer and title

commit ae7a76d
Author: maki49 <[email protected]>
Date:   Fri Sep 26 21:15:49 2025 +0800

    minimize the index access

commit cb68103
Author: maki49 <[email protected]>
Date:   Fri Sep 26 21:03:58 2025 +0800

    wait after each term to save memory

commit 4793d3c
Author: maki49 <[email protected]>
Date:   Wed Sep 24 13:57:26 2025 +0800

    malloc_trim(0)

commit 150b929
Author: maki49 <[email protected]>
Date:   Wed Sep 24 13:43:32 2025 +0800

    fix: filter zero-tensors

commit 3a773e4
Author: maki49 <[email protected]>
Date:   Tue Sep 23 21:27:58 2025 +0800

    symmetry filter

commit 811e6da
Author: maki49 <[email protected]>
Date:   Tue Sep 23 21:27:22 2025 +0800

    default cvc (only for test)

commit 3f4b6b9
Author: maki49 <[email protected]>
Date:   Tue Sep 23 19:02:13 2025 +0800

    omp parallel method 3

commit d91c80d
Author: maki49 <[email protected]>
Date:   Tue Sep 23 17:31:38 2025 +0800

    omp parallel method 2

commit 0a8a529
Author: maki49 <[email protected]>
Date:   Tue Sep 23 16:08:48 2025 +0800

    fix omp bug

commit e329128
Author: maki49 <[email protected]>
Date:   Sat Sep 20 16:58:40 2025 +0800

    CVC (D-perspective) algorithm (draft)

    fix adding tensor maps

    fix compile problems

commit d4f7320
Author: Minye Zhang <[email protected]>
Date:   Mon May 19 17:25:38 2025 +0200

    remove redundant member variables

Author: maki49

include/RI/global/Shared_Vector.h

@@ -28,6 +28,18 @@ 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];
+	}
+	Shape_Vector(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_; }

include/RI/global/Tensor.h

@@ -47,6 +47,11 @@ class Tensor
 	Tensor transpose() const;
 	Tensor dagger() const;
 
+	/// permute from the input index order to {0, 1, 2, ..., N}
+	Tensor permute_from(const std::vector<std::size_t>& order) const;
+	/// permute from {0, 1, 2, ..., N } to the input new index order 
+	// Tensor permute_to(const std::vector<std::size_t>) const;
+
 	// ||d||_p = (|d_1|^p+|d_2|^p+...)^{1/p}
 	// if(p==std::numeric_limits<double>::max())    ||d||_max = max_i |d_i|
 	Global_Func::To_Real_t<T> norm(const double p) const;

include/RI/global/Tensor.hpp

@@ -324,4 +324,74 @@ std::array<std::array<std::array<std::array<T,N3>,N2>,N1>,N0> to_array(const Ten
 	return a;
 }
 
+// check whether an array is a permutation of {from, from+1, ..., to}
+inline bool is_permutation(const std::vector<size_t>& arr, const std::size_t from, const std::size_t to)
+{
+	std::vector<bool> seen(arr.size(), false);
+	if (arr.size() != to - from + 1) return false;
+	for (std::size_t i = 0; i < arr.size(); ++i) {
+		if (arr[i] < from || arr[i] > to) return false;
+		const std::size_t idx = arr[i] - from;
+		if (seen[idx]) return false;
+		seen[idx] = true;
+	}
+	return true;
+}
+
+template<typename T>
+Tensor<T> Tensor<T>::permute_from(const std::vector<std::size_t>& order) const
+{
+	const std::size_t N = order.size();
+	assert(this->shape.size() == N);
+	assert(is_permutation(order, 0, N-1));
+	Shape_Vector new_shape(this->shape);
+	for (std::size_t d = 0; d < N; ++d)
+		new_shape[d] = this->shape[order[d]];
+
+	Tensor<T> t(new_shape);
+
+	/*	N-permutation
+	for (std::size_t i0 = 0; i0 < new_shape[0], ++i0)
+		for (std::size_t i1 = 0, i1 < new_shape[1], ++i1)
+			...
+				for (std::size_t iN = 0, iN < new_shape[N], ++iN)
+					t(i0, i1, ..., iN) = this->operator()(order(i0), order(i1), ..., order(N))
+	*/
+	std::vector<std::size_t> strides(N);
+	strides[N - 1] = 1;
+	for (int d = N - 2; d >= 0; --d)
+		strides[d] = strides[d + 1] * this->shape[d + 1];
+
+	std::vector<std::size_t> new_strides(N);
+	new_strides[N - 1] = 1;
+	for (int d = N - 2; d >= 0; --d)
+		new_strides[d] = new_strides[d + 1] * new_shape[d + 1];
+
+
+	std::vector<std::size_t> new_idx(N, 0);
+	std::vector<std::size_t> idx(N, 0);
+	for (std::size_t i = 0; i < t.get_shape_all(); ++i)	// new 1D-index
+	{
+		// new ND-index
+		std::size_t tmp = i;
+		for (std::size_t d = 0; d < N; ++d)
+		{
+			new_idx[d] = tmp / new_strides[d];
+			tmp %= new_strides[d];
+		}
+		// old ND-index
+		for (std::size_t d = 0; d < N; ++d)
+			idx[order[d]] = new_idx[d];
+
+		// old 1D-index
+		std::size_t i0 = 0;
+		for (std::size_t d = 0; d < N; ++d)
+			i0 += idx[d] * strides[d];
+
+		(*t.data)[i] = (*this->data)[i0];
+	}
+
+	return t;
+}
+
 }

include/RI/global/Tensor_Multiply-32.hpp

@@ -80,6 +80,28 @@ namespace Tensor_Multiply
 		return Txy;
 	}
 
+	// Txy(x0,x1,..., xM) = Tx(x0, x1,x2, ..., xM, y0, y1, ..., yN) * Vy(y0, y1, ..., yN)
+	template<typename Tdata>
+	Tensor<Tdata> gemv(const Tensor<Tdata>& Tx, const Tensor<Tdata>& Vy)
+	{
+		assert(Tx.shape.size() >= Vy.shape.size());
+		const std::size_t ny = Vy.get_shape_all();
+		assert(Tx.get_shape_all() % ny == 0);
+		const std::size_t dim = Tx.shape.size() - Vy.shape.size();
+		std::vector<std::size_t> shape_vector;
+		if (dim == 0)
+			shape_vector.push_back(1);
+		else
+			for (int d = 0; d < dim;++d)
+				shape_vector.push_back(Tx.shape[d]);
+		Tensor<Tdata> Txy(shape_vector);
+		Blas_Interface::gemv(
+			'N', Txy.get_shape_all(), ny,
+			Tdata(1.0), Tx.ptr(), ny, Vy.ptr(), 1,
+			Tdata(0.0), Txy.ptr(), 1);
+		return Txy;
+	}
+
 }
 
 }

include/RI/physics/Exx.h

@@ -126,10 +126,9 @@ class Exx
 	};
 	Flag_Save_Result flag_save_result;
 
-	MPI_Comm mpi_comm;
-	std::map<TA,Tatom_pos> atoms_pos;
-	std::array<Tatom_pos,Ndim> latvec;
-	std::array<Tcell,Ndim> period;
+	std::string method = "loop3";		// "loop3", "cvc"
+	std::map<TA, std::map<TAC, std::map<TA, std::map<TAC, Tensor<Tdata>>>>> cvc_;
+	std::map<TA, std::map<TAC, std::map<TA, std::map<TAC, Tensor<Tdata>>>>>& get_cvc() const { return this->cvc_; }
 };
 
 }

include/RI/physics/Exx.hpp

@@ -18,22 +18,17 @@ namespace RI
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
 void Exx<TA,Tcell,Ndim,Tdata>::set_parallel(
-	const MPI_Comm &mpi_comm_in,
-	const std::map<TA,Tatom_pos> &atoms_pos_in,
-	const std::array<Tatom_pos,Ndim> &latvec_in,
-	const std::array<Tcell,Ndim> &period_in)
+	const MPI_Comm &mpi_comm,
+	const std::map<TA,Tatom_pos> &atoms_pos,
+	const std::array<Tatom_pos,Ndim> &latvec,
+	const std::array<Tcell,Ndim> &period)
 {
-	this->mpi_comm = mpi_comm_in;
-	this->atoms_pos = atoms_pos_in;
-	this->latvec = latvec_in;
-	this->period = period_in;
-
 	this->lri.set_parallel(
-		this->mpi_comm, this->atoms_pos, this->latvec, this->period,
+		mpi_comm, atoms_pos, latvec, period,
 		{Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2, Label::ab_ab::a0b0_a2b1, Label::ab_ab::a0b0_a2b2});
 	this->flag_finish.stru = true;
 	//if()
-		this->post_2D.set_parallel(this->mpi_comm, this->atoms_pos, this->period);
+		this->post_2D.set_parallel(mpi_comm, atoms_pos, period);
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
@@ -43,7 +38,7 @@ void Exx<TA,Tcell,Ndim,Tdata>::set_symmetry(
 {
 	if(flag_symmetry)
 		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, irreducible_sector);
+			this->lri.period, irreducible_sector);
 	else
 		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 }
@@ -271,12 +266,21 @@ void Exx<TA,Tcell,Ndim,Tdata>::cal_Hs(
 
 	if(!this->flag_finish.Ds_delta)
 		this->Hs.clear();
-	this->lri.cal_loop3(
-		{Label::ab_ab::a0b0_a1b1,
-		 Label::ab_ab::a0b0_a1b2,
-		 Label::ab_ab::a0b0_a2b1,
-		 Label::ab_ab::a0b0_a2b2},
-		this->Hs);
+	if (this->method == "loop3")
+		this->lri.cal_loop3(
+			{ Label::ab_ab::a0b0_a1b1,
+			 Label::ab_ab::a0b0_a1b2,
+			 Label::ab_ab::a0b0_a2b1,
+			 Label::ab_ab::a0b0_a2b2 },
+			this->Hs);
+	else if (this->method == "cvc")
+	{
+		if (this->cvc_.empty())
+			this->cvc_ = this->lri.cal_cvc();
+		this->Hs = this->lri.constract_cvc_ds(this->cvc_);
+	}
+	else
+		assert(false && "Unknown method in Exx::cal_Hs");
 
 	//if()
 		this->energy = this->post_2D.cal_energy(

include/RI/physics/GW.h

@@ -76,11 +76,6 @@ class GW
 		bool Gs=false;
 	};
 	Flag_Finish flag_finish;
-
-	MPI_Comm mpi_comm;
-	std::map<TA,Tatom_pos> atoms_pos;
-	std::array<Tatom_pos,Ndim> latvec;
-	std::array<Tcell,Ndim> period;
 };
 
 }

include/RI/physics/GW.hpp

@@ -13,18 +13,13 @@ namespace RI
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
 void GW<TA,Tcell,Ndim,Tdata>::set_parallel(
-	const MPI_Comm &mpi_comm_in,
-	const std::map<TA,Tatom_pos> &atoms_pos_in,
-	const std::array<Tatom_pos,Ndim> &latvec_in,
-	const std::array<Tcell,Ndim> &period_in)
+	const MPI_Comm &mpi_comm,
+	const std::map<TA,Tatom_pos> &atoms_pos,
+	const std::array<Tatom_pos,Ndim> &latvec,
+	const std::array<Tcell,Ndim> &period)
 {
-	this->mpi_comm = mpi_comm_in;
-	this->atoms_pos = atoms_pos_in;
-	this->latvec = latvec_in;
-	this->period = period_in;
-
 	this->lri.set_parallel(
-		this->mpi_comm, this->atoms_pos, this->latvec, this->period,
+		mpi_comm, atoms_pos, latvec, period,
 		{Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2, Label::ab_ab::a0b0_a2b1, Label::ab_ab::a0b0_a2b2});
 	this->flag_finish.stru = true;
 	//if()
@@ -38,7 +33,7 @@ void GW<TA,Tcell,Ndim,Tdata>::set_symmetry(
 {
 	if(flag_symmetry)
 		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, irreducible_sector);
+			this->lri.period, irreducible_sector);
 	else
 		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 }

include/RI/physics/RPA.hpp

@@ -34,7 +34,7 @@ void RPA<TA,Tcell,Ndim,Tdata>::set_symmetry(
 {
 	if(flag_symmetry)
 		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, irreducible_sector);
+			this->lri.period, irreducible_sector);
 	else
 		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 }