Separate GPU symmetry filters and align callbacks with loop order

include/RI/physics/Exx.hpp

@@ -9,7 +9,7 @@
 #include "../ri/Cell_Nearest.h"
 #include "../ri/Label.h"
 #include "../global/Map_Operator.h"
-#include "./symmetry/Filter_Atom_Symmetry.h"
+#include "./symmetry/Filter_Atom_Symmetry_Wrapper.h"
 
 #include <cassert>
 
@@ -42,8 +42,8 @@ void Exx<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
 	if(flag_symmetry)
-		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, irreducible_sector);
+		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>>(
+			this->period, irreducible_sector, this->lri.cal_mode);
 	else
 		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 }

include/RI/physics/GW.hpp

@@ -6,7 +6,7 @@
 
 #include "GW.h"
 #include "../ri/Label.h"
-#include "./symmetry/Filter_Atom_Symmetry.h"
+#include "./symmetry/Filter_Atom_Symmetry_Wrapper.h"
 
 namespace RI
 {
@@ -37,8 +37,8 @@ void GW<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
 	if(flag_symmetry)
-		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, irreducible_sector);
+		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>>(
+			this->period, irreducible_sector, this->lri.cal_mode);
 	else
 		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 }

include/RI/physics/RPA.hpp

@@ -7,7 +7,7 @@
 
 #include "RPA.h"
 #include "../ri/Label.h"
-#include "./symmetry/Filter_Atom_Symmetry.h"
+#include "./symmetry/Filter_Atom_Symmetry_Wrapper.h"
 
 #include <cassert>
 
@@ -38,8 +38,8 @@ void RPA<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
 	if(flag_symmetry)
-		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, irreducible_sector);
+		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>>(
+			this->period, irreducible_sector, this->lri.cal_mode);
 	else
 		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 }

include/RI/physics/symmetry/Filter_Atom_Symmetry.h

@@ -45,14 +45,6 @@ class Filter_Atom_Symmetry: public Filter_Atom<TA, std::pair<TA, TC>>
 				case Label::ab_ab::a1b2_a2b0:
 				case Label::ab_ab::a1b2_a2b1:
 					return !this->symmetry.is_Aa_in_irreducible_sector(A1);
-			// GPU visits Aa01 before the output Aa2 is known in these cases.
-				case Label::ab_ab::a0b0_a1b1:
-				case Label::ab_ab::a0b1_a1b0:
-				case Label::ab_ab::a0b0_a1b2:
-				case Label::ab_ab::a0b1_a1b2:
-				case Label::ab_ab::a0b2_a1b0:
-				case Label::ab_ab::a0b2_a1b1:
-					return false;
 			default:
 				throw std::invalid_argument("label "+Label_Tools::get_name(label)+" error in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
 		}
@@ -70,11 +62,6 @@ class Filter_Atom_Symmetry: public Filter_Atom<TA, std::pair<TA, TC>>
 				case Label::ab_ab::a0b1_a1b2:
 				case Label::ab_ab::a0b2_a1b0:
 				case Label::ab_ab::a0b2_a1b1:
-			// a01b2_a2b01: Ab01 is the output right atom on GPU.
-				case Label::ab_ab::a0b2_a2b0:
-				case Label::ab_ab::a0b2_a2b1:
-				case Label::ab_ab::a1b2_a2b0:
-				case Label::ab_ab::a1b2_a2b1:
 					return !this->symmetry.is_Ab_in_irreducible_sector(A1.first);
 			default:
 				throw std::invalid_argument("label "+Label_Tools::get_name(label)+" error in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
@@ -113,14 +100,6 @@ class Filter_Atom_Symmetry: public Filter_Atom<TA, std::pair<TA, TC>>
 				case Label::ab_ab::a1b0_a2b2:
 				case Label::ab_ab::a1b1_a2b2:
 					return !this->symmetry.in_irreducible_sector(A1, A3);
-			// GPU output is (Aa2, Ab2) or (Aa2, Ab01); Aa01 is intermediate.
-				case Label::ab_ab::a0b0_a1b1:
-				case Label::ab_ab::a0b1_a1b0:
-				case Label::ab_ab::a0b0_a1b2:
-				case Label::ab_ab::a0b1_a1b2:
-				case Label::ab_ab::a0b2_a1b0:
-				case Label::ab_ab::a0b2_a1b1:
-					return !this->symmetry.in_irreducible_sector(A2, A3);
 			// a01b2_a2b01:                                        Aa01,          Ab01,          Aa2
 			default:
 					return false;

include/RI/physics/symmetry/Filter_Atom_Symmetry_GPU.h

@@ -0,0 +1,79 @@
+// ===================
+//  Author: maki49
+//  date: 2026.11.10
+// ===================
+
+#pragma once
+
+#include "../../ri/Filter_Atom.h"
+#include "../../ri/Label_Tools.h"
+#include "Symmetry_Filter.h"
+
+namespace RI
+{
+
+// Arguments follow the GPU loop nesting, from outermost to innermost.
+// Only output atoms may be screened; intermediate atoms are summed over.
+template<typename TA, typename TC, typename Tdata>
+class Filter_Atom_Symmetry_GPU: public Filter_Atom<TA, std::pair<TA, TC>>
+{
+  public:
+    using TAC = std::pair<TA, TC>;
+    using Group = Label::Aab_Aab;
+
+    Filter_Atom_Symmetry_GPU(
+        const TC &period,
+        const std::map<std::pair<TA,TA>, std::set<TC>> &irsec)
+        : symmetry(period, irsec) {}
+
+    bool filter_for1(const Label::ab_ab &label, const TA &A1) const override
+    {
+        switch (Label_Tools::to_Aab_Aab(label))
+        {
+            case Group::a01b01_a01b01: // Aa01 is intermediate.
+            case Group::a01b01_a01b2:
+                return false;
+            default: // Aa01 is the left output atom.
+                return !symmetry.is_Aa_in_irreducible_sector(A1);
+        }
+    }
+
+    bool filter_for1(const Label::ab_ab &, const TAC &A1) const override
+    {
+        // The first loop with a cell-bearing atom visits output Ab01.
+        return !symmetry.is_Ab_in_irreducible_sector(A1.first);
+    }
+
+    bool filter_for2(const Label::ab_ab &label, const TA &A1, const TAC &A2) const override
+    {
+        // Only the final stage of a01b2_a2b01 has both outputs at level 2.
+        if (Label_Tools::to_Aab_Aab(label) == Group::a01b2_a2b01)
+            return !symmetry.in_irreducible_sector(A1, A2);
+        return false;
+    }
+
+    bool filter_for32(const Label::ab_ab &label,
+                      const TA &A1, const TAC &A2, const TAC &A3) const override
+    {
+        switch (Label_Tools::to_Aab_Aab(label))
+        {
+            case Group::a01b01_a01b01: // Aa01, Aa2, Ab2
+                return !symmetry.in_irreducible_sector(A2, A3);
+            case Group::a01b01_a01b2: // Aa01, Ab01, Aa2
+                return !symmetry.in_irreducible_sector(A3, A2);
+            case Group::a01b01_a2b01: // Aa01, Ab01, Ab2
+            case Group::a01b2_a2b01: // Aa01, Aa2, Ab01
+                return !symmetry.in_irreducible_sector(A1, A3);
+            case Group::a01b01_a2b2: // Aa01, Ab01, Ab2
+                return !symmetry.in_irreducible_sector(A1, A2);
+            default:
+                throw std::invalid_argument("Unknown GPU symmetry label");
+        }
+    }
+
+    // filter_for31 builds intermediates and deliberately never screens pairs.
+    // All other callbacks retain the base class's no-op behavior.
+    Symmetry_Filter<TA,TC,Tdata> symmetry;
+};
+
+}

include/RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h

@@ -0,0 +1,101 @@
+// ===================
+//  Author: Peize Lin
+//  date: 2026.11.10
+// ===================
+
+#pragma once
+
+#ifndef __GPU_RI
+
+#include "Filter_Atom_Symmetry.h"
+
+namespace RI
+{
+
+template<typename TA, typename TC, typename Tdata>
+class Filter_Atom_Symmetry_Wrapper: public Filter_Atom_Symmetry<TA,TC,Tdata>
+{
+  public:
+	Filter_Atom_Symmetry_Wrapper(
+		const TC &period,
+		const std::map<std::pair<TA,TA>, std::set<TC>> &irsec,
+		const LRI_Cal_Mode &cal_mode_in)
+		: Filter_Atom_Symmetry<TA,TC,Tdata>(period, irsec){}
+};
+
+}
+
+#else
+
+#include "../../ri/Filter_Atom.h"
+#include "Filter_Atom_Symmetry.h"
+#include "Filter_Atom_Symmetry_GPU.h"
+#include "../../ri/Label_Tools.h"
+
+#include <stdexcept>
+
+namespace RI
+{
+
+template<typename TA, typename TC, typename Tdata>
+class Filter_Atom_Symmetry_Wrapper: public Filter_Atom<TA, std::pair<TA, TC>>
+{
+  public:
+	using TAC = std::pair<TA, TC>;
+
+	Filter_Atom_Symmetry_Wrapper(
+		const TC &period,
+		const std::map<std::pair<TA,TA>, std::set<TC>> &irsec,
+		const LRI_Cal_Mode &cal_mode_in)
+		:cal_mode(cal_mode_in),
+		 filter_CPU(period, irsec),
+		 filter_GPU(period, irsec){}
+
+	virtual bool filter_for1(const Label::ab_ab &label, const TA &A1) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? this->filter_GPU.filter_for1(label, A1)
+			: this->filter_CPU.filter_for1(label, A1);
+	}
+	virtual bool filter_for1(const Label::ab_ab &label, const TAC &A1) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? this->filter_GPU.filter_for1(label, A1)
+			: this->filter_CPU.filter_for1(label, A1);
+	}
+
+	virtual bool filter_for2(const Label::ab_ab &label, const TA &A1, const TAC &A2) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? this->filter_GPU.filter_for2(label, A1, A2)
+			: this->filter_CPU.filter_for2(label, A1, A2);
+	}
+
+	virtual bool filter_for32(const Label::ab_ab &label, const TA &A1, const TAC &A2, const TAC &A3) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? this->filter_GPU.filter_for32(label, A1, A2, A3)
+			: this->filter_CPU.filter_for32(label, A1, A2, A3);
+	}
+	virtual bool filter_for32(const Label::ab_ab &label, const TAC &A1, const TA &A2, const TAC &A3) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? false
+			: this->filter_CPU.filter_for32(label, A1, A2, A3);
+	}
+	virtual bool filter_for32(const Label::ab_ab &label, const TAC &A1, const TAC &A2, const TAC &A3) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? false
+			: this->filter_CPU.filter_for32(label, A1, A2, A3);
+	}
+
+  public:	// private
+	const LRI_Cal_Mode &cal_mode;
+	Filter_Atom_Symmetry<TA,TC,Tdata> filter_CPU;
+	Filter_Atom_Symmetry_GPU<TA,TC,Tdata> filter_GPU;
+};
+
+}
+
+#endif

include/RI/ri/LRI-cal_loop3_GPU.hpp

@@ -1552,8 +1552,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									// The output block is (Aa01, Ab01); Ab2 is only an intermediate atom.
-									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = Global_Func::find(Ds_b_transpose, Ab01.first, TAC{Ab2.first, (Ab2.second-Ab01.second)%this->period});
 									if(D_b.empty())	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
@@ -1692,7 +1691,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
@@ -1831,7 +1830,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = Global_Func::find(Ds_b_transpose, Ab01.first, TAC{Ab2.first, (Ab2.second-Ab01.second)%this->period});
 									if(D_b.empty())	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
@@ -1970,7 +1969,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
@@ -2105,11 +2104,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								const TAC &Aa2 = list_Aa2[ia2];
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
+								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									// Screen the output pair; Aa2 is summed over, not an output atom.
-									if(this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 									if(this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
@@ -2270,11 +2268,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								const TAC &Aa2 = list_Aa2[ia2];
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
+								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									// Screen the output pair; Aa2 is summed over, not an output atom.
-									if(this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 									if(this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
@@ -2435,11 +2432,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								const TAC &Aa2 = list_Aa2[ia2];
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
+								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									// Screen the output pair; Aa2 is summed over, not an output atom.
-									if(this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 									if(this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
@@ -2600,11 +2596,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								const TAC &Aa2 = list_Aa2[ia2];
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
+								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									// Screen the output pair; Aa2 is summed over, not an output atom.
-									if(this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 									if(this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;