Fix real-space symmetry filtering in GPU contractions

include/RI/physics/symmetry/Filter_Atom_Symmetry.h

@@ -45,6 +45,14 @@ 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__));
 		}
@@ -62,6 +70,11 @@ 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__));
@@ -100,6 +113,14 @@ 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/ri/LRI-cal_loop3_GPU.hpp

@@ -1554,7 +1554,8 @@ 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, Ab01, Ab2))	continue;
+									// The output block is (Aa01, Ab01); Ab2 is only an intermediate atom.
+									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	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);
@@ -1693,7 +1694,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, Ab01, Ab2))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	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);
@@ -1832,7 +1833,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, Ab01, Ab2))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	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);
@@ -1971,7 +1972,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, Ab01, Ab2))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	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);
@@ -2106,10 +2107,11 @@ 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,10 +2272,11 @@ 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;
@@ -2434,10 +2437,11 @@ 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;
@@ -2598,10 +2602,11 @@ 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;

include/RI/ri/LRI_Cal_Aux.h

@@ -95,7 +95,7 @@ namespace LRI_Cal_Aux
 		if (D_result.empty())
 		{
 			if(Tdata(1.0)==fac)
-				D_result = D_add;
+				D_result = D_add.copy();
 			else if(Tdata(-1.0)==fac)
 				D_result = -D_add;
 			else

include/RI/ri/LRI_k-cal_cvc_mo.hpp

@@ -467,7 +467,7 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvc_mo_k_hartree_onthefly(
 		{
 			for (int k1_index=0; k1_index<k1_indices.size(); ++k1_index)
 			{
-				const TA& mu = list_J[mu_index];
+				const TA& mu = list_I[mu_index];
 				const int k1 = k1_indices[k1_index];
 				// 2.1 calculate C^mu (i,a^*)[k1,k1] on-the-fly
 				const Tensor<Tdata> C_mu_ia = Cs_ao_mo_to_Cs_mo(Cs_ao_mo, map_psi, k1, k1,

include/RI/ri/LRI_k-cal_hartree.hpp

@@ -121,7 +121,7 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvcd_k_hartree(
 #pragma omp for schedule(dynamic) collapse(2)
 		for (int v_index=0; v_index<list_I.size(); v_index++)
 		{
-			for (int u_index=0; u_index<list_I.size(); u_index++)
+			for (int u_index=0; u_index<list_J.size(); u_index++)
 			{
 				const TA v = list_I[v_index];
 				const TA u = list_J[u_index];

unittests/gpu_symmetry/CMakeLists.txt

@@ -0,0 +1,36 @@
+cmake_minimum_required(VERSION 3.18)
+project(libri_gpu_symmetry_tests LANGUAGES CXX)
+
+option(ENABLE_GPU "Run the DDLA/CUDA contraction checks" OFF)
+set(LIBRI_ROOT "${CMAKE_CURRENT_SOURCE_DIR}/../.." CACHE PATH "LibRI source tree")
+set(LIBCOMM_ROOT "" CACHE PATH "LibComm installation")
+set(CEREAL_ROOT "" CACHE PATH "cereal installation")
+set(DDLA_ROOT "" CACHE PATH "DDLA installation")
+
+find_package(MPI REQUIRED COMPONENTS CXX)
+find_package(OpenMP REQUIRED COMPONENTS CXX)
+find_package(BLAS REQUIRED)
+find_package(LAPACK REQUIRED)
+
+add_executable(test_lri_gpu_symmetry test_lri_gpu_symmetry.cpp)
+target_compile_features(test_lri_gpu_symmetry PRIVATE cxx_std_14)
+target_include_directories(test_lri_gpu_symmetry PRIVATE
+  "${LIBRI_ROOT}/include" "${LIBCOMM_ROOT}/include" "${CEREAL_ROOT}/include")
+target_link_libraries(test_lri_gpu_symmetry PRIVATE
+  MPI::MPI_CXX OpenMP::OpenMP_CXX ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
+
+if(ENABLE_GPU)
+  find_package(CUDAToolkit REQUIRED)
+  find_library(DDLA_LIBRARY NAMES ddla PATHS "${DDLA_ROOT}/lib" REQUIRED NO_DEFAULT_PATH)
+  target_compile_features(test_lri_gpu_symmetry PRIVATE cxx_std_17)
+  target_compile_definitions(test_lri_gpu_symmetry PRIVATE __GPU_RI __DDLA_RI DDLA_USE_CUDA)
+  target_include_directories(test_lri_gpu_symmetry PRIVATE "${DDLA_ROOT}/include")
+  target_link_libraries(test_lri_gpu_symmetry PRIVATE
+    ${DDLA_LIBRARY} CUDA::cudart CUDA::cublas CUDA::curand CUDA::cusolver)
+endif()
+
+enable_testing()
+add_test(NAME lri_symmetry COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 1
+  ${MPIEXEC_PREFLAGS} $<TARGET_FILE:test_lri_gpu_symmetry> ${MPIEXEC_POSTFLAGS})
+set_tests_properties(lri_symmetry PROPERTIES ENVIRONMENT
+  "OMP_NUM_THREADS=1;OPENBLAS_NUM_THREADS=1;MKL_NUM_THREADS=1")

unittests/gpu_symmetry/README.md

@@ -0,0 +1,35 @@
+# GPU symmetry regression
+
+Compare all 18 contractions, individually and together, against full CPU
+contractions with an independently applied output-sector mask. The fixture uses
+two unequal atoms and three periodic cells, real and complex tensors, both CPU
+modes, and optionally DDLA/CUDA. It checks unrestricted, restricted and empty
+output, as well as seeded signed, zero and fractional accumulation.
+
+```sh
+cmake -S unittests/gpu_symmetry -B build-symmetry \
+  -DLIBCOMM_ROOT=/path/to/LibComm -DCEREAL_ROOT=/path/to/cereal \
+  -DENABLE_GPU=ON -DDDLA_ROOT=/path/to/DDLA
+cmake --build build-symmetry --parallel 2
+ctest --test-dir build-symmetry --output-on-failure
+```
+
+The GPU test requires one GPU, one MPI rank, and DDLA built for CUDA. Omit
+`ENABLE_GPU` and `DDLA_ROOT` to check a CPU-only build. The library headers can
+be selected independently with `LIBRI_ROOT`, allowing the same test to reproduce
+failures in an earlier revision.
+
+For OpenMP coverage or an isolated failing contraction:
+
+```sh
+OMP_NUM_THREADS=4 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1 \
+  mpirun -np 1 build-symmetry/test_lri_gpu_symmetry
+OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1 \
+  mpirun -np 1 build-symmetry/test_lri_gpu_symmetry 1 gpu real
+```
+
+Arguments select the contraction index (`0`–`17`, `18` for the combined sum,
+or `-1` for all), backend (`cpu`, `gpu`, `all`) and scalar (`real`, `complex`,
+`all`). Index 1 is `a1b2_a2b1`: filtering `(Aa01,Aa2)` while summing over `Aa2`
+incorrectly removes contributions to the accepted output pair `(Aa01,Ab01)`.
+The off-diagonal sector with a nonzero translation detects this error.