Skip to content

Commit 375686d

Browse files
Dengda98Dengda98
authored
REFAC: use macro to manage static loop of ChnlGrid and IntegGrid (#169)
1 parent 5e645f5 commit 375686d

File tree

13 files changed

+323
-405
lines changed

13 files changed

+323
-405
lines changed

‎pygrt/C_extension/include/grt/common/const.h‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -132,6 +132,16 @@ typedef cplx_t cplxIntegGrid[GRT_SRC_M_NUM][GRT_INTEG_NUM];
132132
typedef real_t realIntegGrid[GRT_SRC_M_NUM][GRT_INTEG_NUM];
133133
typedef size_t sizeIntegGrid[GRT_SRC_M_NUM][GRT_INTEG_NUM];
134134

135+
// 为以上静态多维数组的循环创建宏，谨慎使用，留意作用域
136+
#define GRT_LOOP_ChnlGrid(im, c) \
137+
for(int im = 0; im < GRT_SRC_M_NUM; ++im) \
138+
for(int c = 0; c < GRT_CHANNEL_NUM; ++c)
139+
140+
#define GRT_LOOP_IntegGrid(im, v) \
141+
for(int im = 0; im < GRT_SRC_M_NUM; ++im) \
142+
for(int v = 0; v < GRT_INTEG_NUM; ++v)
143+
144+
135145
/** 不同震源类型在大小为 GRT_SRC_M_NUM 的数组中的索引 */
136146
enum {
137147
GRT_SRC_M_EX_INDEX = 0,

‎pygrt/C_extension/src/common/dwm.c‎

Lines changed: 12 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -67,23 +67,21 @@ real_t grt_discrete_integ(
6767
for(size_t ir=0; ir<nr; ++ir){
6868
if(iendkrs[ir]) continue; // 该震中距下的波数k积分已收敛
6969

70-
memset(SUM, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_INTEG_NUM);
70+
memset(SUM, 0, sizeof(cplxIntegGrid));
7171

7272
// 计算被积函数一项 F(k,w)Jm(kr)k
7373
grt_int_Pk(k, rs[ir], QWV, false, SUM);
7474

7575
iendk0 = true;
76-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
77-
int modr = GRT_SRC_M_ORDERS[i];
7876

79-
for(int v=0; v<GRT_INTEG_NUM; ++v){
80-
sum_J[ir][i][v] += SUM[i][v];
77+
GRT_LOOP_IntegGrid(im, v){
78+
int modr = GRT_SRC_M_ORDERS[im];
79+
sum_J[ir][im][v] += SUM[im][v];
8180

82-
// 是否提前判断达到收敛
83-
if(keps <= 0.0 || (modr==0 && v!=0 && v!=2)) continue;
84-
85-
iendk0 = iendk0 && (fabs(SUM[i][v])/ fabs(sum_J[ir][i][v]) <= keps);
86-
}
81+
// 是否提前判断达到收敛
82+
if(keps <= 0.0 || (modr==0 && v!=0 && v!=2)) continue;
83+
84+
iendk0 = iendk0 && (fabs(SUM[im][v])/ fabs(sum_J[ir][im][v]) <= keps);
8785
}
8886

8987
if(keps > 0.0){
@@ -101,21 +99,17 @@ real_t grt_discrete_integ(
10199
grt_int_Pk(k, rs[ir], QWV_uiz, false, SUM);
102100

103101
// keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛
104-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
105-
for(int v=0; v<GRT_INTEG_NUM; ++v){
106-
sum_uiz_J[ir][i][v] += SUM[i][v];
107-
}
102+
GRT_LOOP_IntegGrid(im, v){
103+
sum_uiz_J[ir][im][v] += SUM[im][v];
108104
}
109105

110106
// ------------------------------- ui_r -----------------------------------
111107
// 计算被积函数一项 F(k,w)Jm(kr)k
112108
grt_int_Pk(k, rs[ir], QWV, true, SUM);
113109

114110
// keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛
115-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
116-
for(int v=0; v<GRT_INTEG_NUM; ++v){
117-
sum_uir_J[ir][i][v] += SUM[i][v];
118-
}
111+
GRT_LOOP_IntegGrid(im, v){
112+
sum_uir_J[ir][im][v] += SUM[im][v];
119113
}
120114
} // END if calc_upar
121115

‎pygrt/C_extension/src/common/fim.c‎

Lines changed: 30 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
#include <stdio.h>
1414
#include <complex.h>
1515
#include <stdlib.h>
16+
#include <string.h>
1617

1718
#include "grt/common/fim.h"
1819
#include "grt/common/integral.h"
@@ -71,27 +72,21 @@ real_t grt_linear_filon_integ(
7172
for(size_t ir=0; ir<nr; ++ir){
7273
if(iendkrs[ir]) continue; // 该震中距下的波数k积分已收敛
7374

74-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
75-
for(int v=0; v<GRT_INTEG_NUM; ++v){
76-
SUM[i][v] = 0.0;
77-
}
78-
}
75+
memset(SUM, 0, sizeof(cplxIntegGrid));
7976

8077
// F(k, w)*Jm(kr)k 的近似公式, sqrt(k) * F(k,w) * cos
8178
grt_int_Pk_filon(k, rs[ir], true, QWV, false, SUM);
8279

8380
iendk0 = true;
84-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
85-
int modr = GRT_SRC_M_ORDERS[i];
8681

87-
for(int v=0; v<GRT_INTEG_NUM; ++v){
88-
sum_J[ir][i][v] += SUM[i][v];
89-
90-
// 是否提前判断达到收敛
91-
if(keps <= 0.0 || (modr==0 && v!=0 && v!=2)) continue;
82+
GRT_LOOP_IntegGrid(im, v){
83+
int modr = GRT_SRC_M_ORDERS[im];
84+
sum_J[ir][im][v] += SUM[im][v];
9285

93-
iendk0 = iendk0 && (fabs(SUM[i][v])/ fabs(sum_J[ir][i][v]) <= keps);
94-
}
86+
// 是否提前判断达到收敛
87+
if(keps <= 0.0 || (modr==0 && v!=0 && v!=2)) continue;
88+
89+
iendk0 = iendk0 && (fabs(SUM[im][v])/ fabs(sum_J[ir][im][v]) <= keps);
9590
}
9691

9792
if(keps > 0.0){
@@ -109,21 +104,17 @@ real_t grt_linear_filon_integ(
109104
grt_int_Pk_filon(k, rs[ir], true, QWV_uiz, false, SUM);
110105

111106
// keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛
112-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
113-
for(int v=0; v<GRT_INTEG_NUM; ++v){
114-
sum_uiz_J[ir][i][v] += SUM[i][v];
115-
}
107+
GRT_LOOP_IntegGrid(im, v){
108+
sum_uiz_J[ir][im][v] += SUM[im][v];
116109
}
117110

118111
// ------------------------------- ui_r -----------------------------------
119112
// 计算被积函数一项 F(k,w)Jm(kr)k
120113
grt_int_Pk_filon(k, rs[ir], true, QWV, true, SUM);
121114

122115
// keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uir也收敛
123-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
124-
for(int v=0; v<GRT_INTEG_NUM; ++v){
125-
sum_uir_J[ir][i][v] += SUM[i][v];
126-
}
116+
GRT_LOOP_IntegGrid(im, v){
117+
sum_uir_J[ir][im][v] += SUM[im][v];
127118
}
128119
} // END if calc_upar
129120

@@ -143,14 +134,12 @@ real_t grt_linear_filon_integ(
143134
for(size_t ir=0; ir<nr; ++ir){
144135
real_t tmp = 2.0*(1.0 - cos(dk*rs[ir])) / (rs[ir]*rs[ir]*dk);
145136

146-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
147-
for(int v=0; v<GRT_INTEG_NUM; ++v){
148-
sum_J[ir][i][v] *= tmp;
137+
GRT_LOOP_IntegGrid(im, v){
138+
sum_J[ir][im][v] *= tmp;
149139

150-
if(calc_upar){
151-
sum_uiz_J[ir][i][v] *= tmp;
152-
sum_uir_J[ir][i][v] *= tmp;
153-
}
140+
if(calc_upar){
141+
sum_uiz_J[ir][im][v] *= tmp;
142+
sum_uir_J[ir][im][v] *= tmp;
154143
}
155144
}
156145
}
@@ -189,10 +178,8 @@ real_t grt_linear_filon_integ(
189178
real_t tmpc = tmp * (1.0 - cos(dk*rs[ir]));
190179
real_t tmps = sgn * tmp * sin(dk*rs[ir]);
191180

192-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
193-
for(int v=0; v<GRT_INTEG_NUM; ++v){
194-
sum_J[ir][i][v] += (- tmpc*SUM_Gc[iik][i][v] + tmps*SUM_Gs[iik][i][v] - sgn*SUM_Gs[iik][i][v]/rs[ir]);
195-
}
181+
GRT_LOOP_IntegGrid(im, v){
182+
sum_J[ir][im][v] += (- tmpc*SUM_Gc[iik][im][v] + tmps*SUM_Gs[iik][im][v] - sgn*SUM_Gs[iik][im][v]/rs[ir]);
196183
}
197184

198185
// ---------------- 位移空间导数，SUM_Gc/s数组重复利用 --------------------------
@@ -205,13 +192,10 @@ real_t grt_linear_filon_integ(
205192
// Gs
206193
grt_int_Pk_filon(k0N, rs[ir], false, QWV_uiz, false, SUM_Gs[iik]);
207194

208-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
209-
for(int v=0; v<GRT_INTEG_NUM; ++v){
210-
sum_uiz_J[ir][i][v] += (- tmpc*SUM_Gc[iik][i][v] + tmps*SUM_Gs[iik][i][v] - sgn*SUM_Gs[iik][i][v]/rs[ir]);
211-
}
195+
GRT_LOOP_IntegGrid(im, v){
196+
sum_uiz_J[ir][im][v] += (- tmpc*SUM_Gc[iik][im][v] + tmps*SUM_Gs[iik][im][v] - sgn*SUM_Gs[iik][im][v]/rs[ir]);
212197
}
213198

214-
215199
// ------------------------------- ui_r -----------------------------------
216200
// 计算被积函数一项 F(k,w)Jm(kr)k
217201
// Gc
@@ -220,10 +204,8 @@ real_t grt_linear_filon_integ(
220204
// Gs
221205
grt_int_Pk_filon(k0N, rs[ir], false, QWV, true, SUM_Gs[iik]);
222206

223-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
224-
for(int v=0; v<GRT_INTEG_NUM; ++v){
225-
sum_uir_J[ir][i][v] += (- tmpc*SUM_Gc[iik][i][v] + tmps*SUM_Gs[iik][i][v] - sgn*SUM_Gs[iik][i][v]/rs[ir]);
226-
}
207+
GRT_LOOP_IntegGrid(im, v){
208+
sum_uir_J[ir][im][v] += (- tmpc*SUM_Gc[iik][im][v] + tmps*SUM_Gs[iik][im][v] - sgn*SUM_Gs[iik][im][v]/rs[ir]);
227209
}
228210
} // END if calc_upar
229211

@@ -234,14 +216,13 @@ real_t grt_linear_filon_integ(
234216
// 乘上总系数 sqrt(2.0/(PI*r)) / dk0, 除dks0是在该函数外还会再乘dk0, 并将结果加到原数组中
235217
for(size_t ir=0; ir<nr; ++ir){
236218
real_t tmp = sqrt(2.0/(PI*rs[ir])) / dk0;
237-
for(int i=0; i<GRT_SRC_M_NUM; ++i){
238-
for(int v=0; v<GRT_INTEG_NUM; ++v){
239-
sum_J0[ir][i][v] += sum_J[ir][i][v] * tmp;
240219

241-
if(calc_upar){
242-
sum_uiz_J0[ir][i][v] += sum_uiz_J[ir][i][v] * tmp;
243-
sum_uir_J0[ir][i][v] += sum_uir_J[ir][i][v] * tmp;
244-
}
220+
GRT_LOOP_IntegGrid(im, v){
221+
sum_J0[ir][im][v] += sum_J[ir][im][v] * tmp;
222+
223+
if(calc_upar){
224+
sum_uiz_J0[ir][im][v] += sum_uiz_J[ir][im][v] * tmp;
225+
sum_uir_J0[ir][im][v] += sum_uir_J[ir][im][v] * tmp;
245226
}
246227
}
247228
}

‎pygrt/C_extension/src/common/integral.c‎

Lines changed: 6 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -162,17 +162,14 @@ void grt_int_Pk_sa_filon(const real_t k3[3], real_t r, const cplxChnlGrid QWV3[3
162162
cplxChnlGrid quad_a={0};
163163
cplxChnlGrid quad_b={0};
164164
cplxChnlGrid quad_c={0};
165-
for(int im=0; im<GRT_SRC_M_NUM; ++im){
165+
GRT_LOOP_ChnlGrid(im, c){
166166
int modr = GRT_SRC_M_ORDERS[im];
167-
for(int c=0; c<GRT_CHANNEL_NUM; ++c){
168-
if(modr==0 && GRT_QWV_CODES[c] == 'v') continue;
169-
170-
cplx_t F3[3];
171-
for(int d=0; d<3; ++d) F3[d] = QWV3[d][im][c] * sqrt(k3[d]) * sgn;
167+
if(modr==0 && GRT_QWV_CODES[c] == 'v') continue;
168+
cplx_t F3[3];
169+
for(int d=0; d<3; ++d) F3[d] = QWV3[d][im][c] * sqrt(k3[d]) * sgn;
172170

173-
// 拟合参数
174-
grt_quad_term(k3, F3, &quad_a[im][c], &quad_b[im][c], &quad_c[im][c]);
175-
}
171+
// 拟合参数
172+
grt_quad_term(k3, F3, &quad_a[im][c], &quad_b[im][c], &quad_c[im][c]);
176173
}
177174

178175
real_t kmin = k3[0], kmax = k3[2];

‎pygrt/C_extension/src/common/iostats.c‎

Lines changed: 39 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -21,13 +21,11 @@ void grt_write_stats(FILE *f0, real_t k, const cplxChnlGrid QWV)
2121
{
2222
fwrite(&k, sizeof(real_t), 1, f0);
2323

24-
for(int im=0; im<GRT_SRC_M_NUM; ++im){
24+
GRT_LOOP_ChnlGrid(im, c){
2525
int modr = GRT_SRC_M_ORDERS[im];
26-
for(int c=0; c<GRT_CHANNEL_NUM; ++c){
27-
if(modr == 0 && GRT_QWV_CODES[c] == 'v') continue;
26+
if(modr == 0 && GRT_QWV_CODES[c] == 'v') continue;
2827

29-
fwrite(&QWV[im][c], sizeof(cplx_t), 1, f0);
30-
}
28+
fwrite(&QWV[im][c], sizeof(cplx_t), 1, f0);
3129
}
3230
}
3331

@@ -39,14 +37,12 @@ int grt_extract_stats(FILE *bf0, FILE *af0){
3937
snprintf(K, sizeof(K), GRT_STRING_FMT, "k"); K[0]=GRT_COMMENT_HEAD;
4038
fprintf(af0, "%s", K);
4139

42-
for(int im=0; im<GRT_SRC_M_NUM; ++im){
40+
GRT_LOOP_ChnlGrid(im, c){
4341
int modr = GRT_SRC_M_ORDERS[im];
44-
for(int c=0; c<GRT_CHANNEL_NUM; ++c){
45-
if(modr == 0 && GRT_QWV_CODES[c] == 'v') continue;
42+
if(modr == 0 && GRT_QWV_CODES[c] == 'v') continue;
4643

47-
snprintf(K, sizeof(K), "%s_%c", GRT_SRC_M_NAME_ABBR[im], GRT_QWV_CODES[c]);
48-
fprintf(af0, GRT_STR_CMPLX_FMT, K);
49-
}
44+
snprintf(K, sizeof(K), "%s_%c", GRT_SRC_M_NAME_ABBR[im], GRT_QWV_CODES[c]);
45+
fprintf(af0, GRT_STR_CMPLX_FMT, K);
5046
}
5147

5248
return 0;
@@ -58,14 +54,12 @@ int grt_extract_stats(FILE *bf0, FILE *af0){
5854
if(1 != fread(&k, sizeof(real_t), 1, bf0)) return -1;
5955
fprintf(af0, GRT_REAL_FMT, k);
6056

61-
for(int im=0; im<GRT_SRC_M_NUM; ++im){
57+
GRT_LOOP_ChnlGrid(im, c){
6258
int modr = GRT_SRC_M_ORDERS[im];
63-
for(int c=0; c<GRT_CHANNEL_NUM; ++c){
64-
if(modr == 0 && GRT_QWV_CODES[c] == 'v') continue;
65-
66-
if(1 != fread(&val, sizeof(cplx_t), 1, bf0)) return -1;
67-
fprintf(af0, GRT_CMPLX_FMT, creal(val), cimag(val));
68-
}
59+
if(modr == 0 && GRT_QWV_CODES[c] == 'v') continue;
60+
61+
if(1 != fread(&val, sizeof(cplx_t), 1, bf0)) return -1;
62+
fprintf(af0, GRT_CMPLX_FMT, creal(val), cimag(val));
6963
}
7064

7165
return 0;
@@ -75,21 +69,18 @@ int grt_extract_stats(FILE *bf0, FILE *af0){
7569
void grt_write_stats_ptam(
7670
FILE *f0,
7771
realIntegGrid Kpt[GRT_PTAM_PT_MAX],
78-
cplxIntegGrid Fpt[GRT_PTAM_PT_MAX]
79-
){
80-
72+
cplxIntegGrid Fpt[GRT_PTAM_PT_MAX])
73+
{
8174
for(int i=0; i<GRT_PTAM_PT_MAX; ++i){
82-
for(int im=0; im<GRT_SRC_M_NUM; ++im){
75+
76+
GRT_LOOP_IntegGrid(im, v){
8377
int modr = GRT_SRC_M_ORDERS[im];
84-
for(int v=0; v<GRT_INTEG_NUM; ++v){
85-
if(modr == 0 && v!=0 && v!=2) continue;
86-
87-
fwrite(&Kpt[i][im][v], sizeof(real_t), 1, f0);
88-
fwrite(&Fpt[i][im][v], sizeof(cplx_t), 1, f0);
89-
}
78+
if(modr == 0 && v!=0 && v!=2) continue;
79+
80+
fwrite(&Kpt[i][im][v], sizeof(real_t), 1, f0);
81+
fwrite(&Fpt[i][im][v], sizeof(cplx_t), 1, f0);
9082
}
9183
}
92-
9384
}
9485

9586

@@ -99,23 +90,21 @@ int grt_extract_stats_ptam(FILE *bf0, FILE *af0){
9990
char K[20], K2[20];
10091
int icol=0;
10192

102-
for(int im=0; im<GRT_SRC_M_NUM; ++im){
93+
GRT_LOOP_IntegGrid(im, v){
10394
int modr = GRT_SRC_M_ORDERS[im];
104-
for(int v=0; v<GRT_INTEG_NUM; ++v){
105-
if(modr == 0 && v!=0 && v!=2) continue;
106-
107-
snprintf(K2, sizeof(K2), "sum_%s_%d_k", GRT_SRC_M_NAME_ABBR[im], v);
108-
if(icol==0){
109-
snprintf(K, sizeof(K), GRT_STRING_FMT, K2); K2[0]=GRT_COMMENT_HEAD;
110-
fprintf(af0, "%s", K);
111-
} else {
112-
fprintf(af0, GRT_STRING_FMT, K2);
113-
}
114-
snprintf(K2, sizeof(K2), "sum_%s_%d", GRT_SRC_M_NAME_ABBR[im], v);
115-
fprintf(af0, GRT_STR_CMPLX_FMT, K2);
116-
117-
icol++;
95+
if(modr == 0 && v!=0 && v!=2) continue;
96+
97+
snprintf(K2, sizeof(K2), "sum_%s_%d_k", GRT_SRC_M_NAME_ABBR[im], v);
98+
if(icol==0){
99+
snprintf(K, sizeof(K), GRT_STRING_FMT, K2); K2[0]=GRT_COMMENT_HEAD;
100+
fprintf(af0, "%s", K);
101+
} else {
102+
fprintf(af0, GRT_STRING_FMT, K2);
118103
}
104+
snprintf(K2, sizeof(K2), "sum_%s_%d", GRT_SRC_M_NAME_ABBR[im], v);
105+
fprintf(af0, GRT_STR_CMPLX_FMT, K2);
106+
107+
icol++;
119108
}
120109

121110
return 0;
@@ -125,16 +114,14 @@ int grt_extract_stats_ptam(FILE *bf0, FILE *af0){
125114
real_t k;
126115
cplx_t val;
127116

128-
for(int im=0; im<GRT_SRC_M_NUM; ++im){
117+
GRT_LOOP_IntegGrid(im, v){
129118
int modr = GRT_SRC_M_ORDERS[im];
130-
for(int v=0; v<GRT_INTEG_NUM; ++v){
131-
if(modr == 0 && v!=0 && v!=2) continue;
119+
if(modr == 0 && v!=0 && v!=2) continue;
132120

133-
if(1 != fread(&k, sizeof(real_t), 1, bf0)) return -1;
134-
fprintf(af0, GRT_REAL_FMT, k);
135-
if(1 != fread(&val, sizeof(cplx_t), 1, bf0)) return -1;
136-
fprintf(af0, GRT_CMPLX_FMT, creal(val), cimag(val));
137-
}
121+
if(1 != fread(&k, sizeof(real_t), 1, bf0)) return -1;
122+
fprintf(af0, GRT_REAL_FMT, k);
123+
if(1 != fread(&val, sizeof(cplx_t), 1, bf0)) return -1;
124+
fprintf(af0, GRT_CMPLX_FMT, creal(val), cimag(val));
138125
}
139126

140127
return 0;

0 commit comments

Comments
 (0)