Skip to content

Commit 2097fa6

Browse files
PeizeLinPeizeLin
committed
1. update LRI::cal_loop3()
1 parent 2a57ade commit 2097fa6

File tree

4 files changed

+137
-10
lines changed

4 files changed

+137
-10
lines changed

‎include/RI/physics/RPA.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -64,7 +64,7 @@ void RPA<TA,Tcell,Ndim,Tdata>::cal_chi0s(
6464

6565
set_Gs_a1(Gs_tau_positive);
6666
set_Gs_a2(Gs_tau_negative);
67-
this->lri.cal({
67+
this->lri.cal_loop3({
6868
Label::ab_ab::a1b1_a2b2,
6969
Label::ab_ab::a1b2_a2b1},
7070
chi0s_vec);
@@ -73,7 +73,7 @@ void RPA<TA,Tcell,Ndim,Tdata>::cal_chi0s(
7373
set_Gs_a2(Gs_tau_positive); // tmp
7474
//set_Gs_a1(conj(Gs_tau_negative));
7575
//set_Gs_a2(conj(Gs_tau_positive));
76-
this->lri.cal({
76+
this->lri.cal_loop3({
7777
Label::ab_ab::a1b1_a2b2,
7878
Label::ab_ab::a1b2_a2b1},
7979
chi0s_vec);

‎include/RI/ri/LRI-cal_loop3.hpp‎

Lines changed: 126 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -29,22 +29,22 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
2929
if(Ds_result.empty())
3030
Ds_result.resize(1);
3131

32-
const bool flag_D_b_transpose = [&labels]() -> bool
32+
const bool flag_D_a_transpose = [&labels]() -> bool
3333
{
3434
for(const Label::ab_ab &label : labels)
3535
switch(label)
3636
{
37-
case Label::ab_ab::a0b0_a2b1: case Label::ab_ab::a0b0_a2b2:
37+
case Label::ab_ab::a0b0_a1b2: case Label::ab_ab::a0b0_a2b2: case Label::ab_ab::a1b2_a2b1:
3838
return true;
3939
}
4040
return false;
4141
}();
42-
const bool flag_D_a_transpose = [&labels]() -> bool
42+
const bool flag_D_b_transpose = [&labels]() -> bool
4343
{
4444
for(const Label::ab_ab &label : labels)
4545
switch(label)
4646
{
47-
case Label::ab_ab::a0b0_a1b2: case Label::ab_ab::a0b0_a2b2:
47+
case Label::ab_ab::a0b0_a2b1: case Label::ab_ab::a0b0_a2b2:
4848
return true;
4949
}
5050
return false;
@@ -326,6 +326,128 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
326326
break;
327327
} // end case a0b0_a2b2
328328

329+
case Label::ab_ab::a1b1_a2b2:
330+
{
331+
const std::vector<TA > list_Aa01 = LRI_Cal_Aux::filter_list_map(
332+
list_Aa01_Da,
333+
this->Ds_ab[Label::ab::a0b0] );
334+
const std::vector<TAC> list_Aa2 = LRI_Cal_Aux::filter_list_map(
335+
list_Aa2_Da,
336+
this->Ds_ab[Label::ab::a2b2] );
337+
const std::vector<TAC> list_Ab01 = LRI_Cal_Aux::filter_list_set(
338+
list_Ab01_Db,
339+
this->index_Ds_ab[Label::ab::a0b0][0]);
340+
const std::vector<TAC> list_Ab2 = LRI_Cal_Aux::filter_list_set(
341+
list_Ab2_Db,
342+
this->index_Ds_ab[Label::ab::a2b2][0]);
343+
344+
for(const TA &Aa01 : list_Aa01)
345+
{
346+
std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
347+
348+
#pragma omp for schedule(dynamic) nowait
349+
for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
350+
{
351+
const TAC &Ab2 = list_Ab2[ib2];
352+
// D_mul = D_a * D_a2b2
353+
Tensor<Tdata> D_mul;
354+
for(const TAC &Aa2 : list_Aa2)
355+
{
356+
const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
357+
if(D_a.empty()) continue;
358+
const Tensor<Tdata> D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
359+
if(D_a2b2.empty()) continue;
360+
361+
// b2a0a1 = a2b2 * a0a1a2
362+
Tensor<Tdata> D_tmp1 = Tensor_Multiply::x1y0y1_ax1_y0y1a(D_a2b2, D_a);
363+
LRI_Cal_Aux::add_Ds(std::move(D_tmp1), D_mul);
364+
}
365+
if(D_mul.empty()) continue;
366+
367+
// D_result = D_mul * D_a1b1 * D_b
368+
for(const TAC &Ab01 : list_Ab01)
369+
{
370+
const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
371+
if(D_b.empty()) continue;
372+
const Tensor<Tdata> D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
373+
if(D_a1b1.empty()) continue;
374+
375+
// b1b2a0 = a1b1 * b2a0a1
376+
const Tensor<Tdata> D_tmp2 = Tensor_Multiply::x1y0y1_ax1_y0y1a(D_a1b1, D_mul);
377+
// a0b0 = b1b2a0 * b0b1b2
378+
Tensor<Tdata> D_tmp3 = Tensor_Multiply::x2y0_abx2_y0ab(D_tmp2, D_b);
379+
LRI_Cal_Aux::add_Ds(std::move(D_tmp3), Ds_result_fixed[Ab01]);
380+
}
381+
} // end for Ab01
382+
383+
if(!Ds_result_fixed.empty())
384+
LRI_Cal_Aux::add_Ds( std::move(Ds_result_fixed),
385+
Ds_result_thread[0][Aa01]);
386+
LRI_Cal_Aux::add_Ds_omp_try(std::move(Ds_result_thread), Ds_result, lock_Ds_result_add, fac_add_Ds);
387+
} // end for Aa01
388+
break;
389+
} // end case a0b0_a2b2
390+
391+
case Label::ab_ab::a1b2_a2b1:
392+
{
393+
const std::vector<TA > list_Aa01 = LRI_Cal_Aux::filter_list_map(
394+
list_Aa01_Da,
395+
this->Ds_ab[Label::ab::a1b2] );
396+
const std::vector<TAC> &list_Aa2 = LRI_Cal_Aux::filter_list_map(
397+
list_Aa2_Da,
398+
this->Ds_ab[Label::ab::a2b1] );
399+
const std::vector<TAC> list_Ab01 = LRI_Cal_Aux::filter_list_set(
400+
list_Ab01_Db,
401+
this->index_Ds_ab[Label::ab::a2b1][0]);
402+
const std::vector<TAC> list_Ab2 = LRI_Cal_Aux::filter_list_set(
403+
list_Ab2_Db,
404+
this->index_Ds_ab[Label::ab::a1b2][0]);
405+
406+
for(const TA &Aa01 : list_Aa01)
407+
{
408+
#pragma omp for schedule(dynamic) nowait
409+
for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
410+
{
411+
const TAC &Ab01 = list_Ab01[ib01];
412+
// D_mul1 = D_b * D_a1b2
413+
Tensor<Tdata> D_mul1;
414+
for(const TAC &Ab2 : list_Ab2)
415+
{
416+
const Tensor<Tdata> D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
417+
if(D_b.empty()) continue;
418+
const Tensor<Tdata> D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
419+
if(D_a1b2.empty()) continue;
420+
421+
// b0b1a1 = b0b1b2 * a1b2
422+
Tensor<Tdata> D_tmp1 = Tensor_Multiply::x0x1y0_x0x1a_y0a(D_b, D_a1b2);
423+
LRI_Cal_Aux::add_Ds(std::move(D_tmp1), D_mul1);
424+
}
425+
if(D_mul1.empty()) continue;
426+
427+
// D_mul2 = D_a2b1 * D_a
428+
Tensor<Tdata> D_mul2;
429+
for(const TAC &Aa2 : list_Aa2)
430+
{
431+
const Tensor<Tdata> &D_a_transpose = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
432+
if(D_a_transpose.empty()) continue;
433+
const Tensor<Tdata> D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
434+
if(D_a2b1.empty()) continue;
435+
// b1a1a0 = a2b1 * a1a0a2
436+
Tensor<Tdata> D_tmp2 = Tensor_Multiply::x1y0y1_ax1_y0y1a(D_a2b1, D_a_transpose);
437+
LRI_Cal_Aux::add_Ds(std::move(D_tmp2), D_mul2);
438+
}
439+
if(D_mul2.empty()) continue;
440+
441+
// D_result = D_mul2 * D_mul1
442+
// a0b0 = b1a1a0 * b0b1a1
443+
Ds_result_thread[0][Aa01][Ab01] = Tensor_Multiply::x2y0_abx2_y0ab(D_mul2, D_mul1);
444+
} // end for Aa01
445+
446+
LRI_Cal_Aux::add_Ds_omp_try(std::move(Ds_result_thread), Ds_result, lock_Ds_result_add, fac_add_Ds);
447+
} // end for Ab01
448+
break;
449+
} // end case a1b2_a2b1
450+
329451
default:
330452
throw std::invalid_argument(std::string(__FILE__)+std::to_string(__LINE__));
331453
} // end switch(label)

‎unittests/Test_All.hpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
#include "unittests/ri/Cell_Nearest-test.hpp"
2121
#include "unittests/physics/Exx-test.hpp"
2222
#include "unittests/physics/RPA-test.hpp"
23+
#include "unittests/physics/GW-test.hpp"
2324

2425
namespace Test_All
2526
{
@@ -84,5 +85,10 @@ namespace Test_All
8485
RPA_Test::main<double>(argc, argv);
8586
RPA_Test::main<std::complex<float>>(argc, argv);
8687
RPA_Test::main<std::complex<double>>(argc, argv);
88+
89+
GW_Test::main<float>(argc, argv);
90+
GW_Test::main<double>(argc, argv);
91+
GW_Test::main<std::complex<float>>(argc, argv);
92+
GW_Test::main<std::complex<double>>(argc, argv);
8793
}
8894
}

‎unittests/ri/LRI_Loop3-test.hpp‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -289,7 +289,7 @@ namespace LRI_Loop3_Test
289289
* Ds_ab[RI::Label::ab::b][Ab01][{Ab2,{0}}](ib0,ib1,ib2);
290290
std::cout<<"a1b0_a2b2\t"<<(Ds_result[0][Aa01][{Ab01,{0}}] - D_test).norm(2)<<std::endl;
291291
}
292-
292+
*/
293293
{
294294
std::vector<T_Ds> Ds_result;
295295
lri.cal_loop3({RI::Label::ab_ab::a1b1_a2b2}, Ds_result);
@@ -302,7 +302,7 @@ namespace LRI_Loop3_Test
302302
* Ds_ab[RI::Label::ab::b][Ab01][{Ab2,{0}}](ib0,ib1,ib2);
303303
std::cout<<"a1b1_a2b2\t"<<(Ds_result[0][Aa01][{Ab01,{0}}] - D_test).norm(2)<<std::endl;
304304
}
305-
305+
/*
306306
{
307307
std::vector<T_Ds> Ds_result;
308308
lri.cal_loop3({RI::Label::ab_ab::a1b2_a2b0}, Ds_result);
@@ -315,7 +315,7 @@ namespace LRI_Loop3_Test
315315
* Ds_ab[RI::Label::ab::b][Ab01][{Ab2,{0}}](ib0,ib1,ib2);
316316
std::cout<<"a1b2_a2b0\t"<<(Ds_result[0][Aa01][{Ab01,{0}}] - D_test).norm(2)<<std::endl;
317317
}
318-
318+
*/
319319
{
320320
std::vector<T_Ds> Ds_result;
321321
lri.cal_loop3({RI::Label::ab_ab::a1b2_a2b1}, Ds_result);
@@ -328,7 +328,6 @@ namespace LRI_Loop3_Test
328328
* Ds_ab[RI::Label::ab::b][Ab01][{Ab2,{0}}](ib0,ib1,ib2);
329329
std::cout<<"a1b2_a2b1\t"<<(Ds_result[0][Aa01][{Ab01,{0}}] - D_test).norm(2)<<std::endl;
330330
}
331-
*/
332331

333332
MPI_Finalize();
334333
}

0 commit comments

Comments
 (0)