• R/O
  • HTTP
  • SSH
  • HTTPS

python: 提交

libtetrabz python package


Commit MetaInfo

修訂704a16eb12039769d615179b86b7cb759ff7d246 (tree)
時間2022-03-18 14:13:50
作者Mitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Backup. It orks fine.

Change Summary

差異

--- a/src/c/libtetrabz_dblstep.c
+++ b/src/c/libtetrabz_dblstep.c
@@ -33,17 +33,17 @@ void libtetrabz_dblstep2(
3333 /*
3434 Tetrahedron method for theta( - de)
3535 */
36- int ii, jj, ib, indx[4];
36+ int i4, j4, ib, indx[4];
3737 double v, thr, e[4], tsmall[4][4];
3838
3939 thr = 1.0e-8;
4040
4141 for (ib = 0; ib < nb; ib++) {
4242
43- for (ii = 0; ii < 4; ii++)
44- w1[ib][ii] = 0.0;
43+ for (i4 = 0; i4 < 4; i4++)
44+ w1[ib][i4] = 0.0;
4545
46- for (ii = 0; ii < 4; ii++) e[ii] = - ei[ii] + ej[ib][ii];
46+ for (i4 = 0; i4 < 4; i4++) e[i4] = - ei[i4] + ej[ib][i4];
4747 eig_sort(4, e, indx);
4848
4949 if (fabs(e[0]) < thr && fabs(e[3]) < thr) {
@@ -51,52 +51,52 @@ void libtetrabz_dblstep2(
5151 Theta(0) = 0.5
5252 */
5353 v = 0.5;
54- for (ii = 0; ii < 4; ii++)
55- w1[ib][ii] += v * 0.25;
54+ for (i4 = 0; i4 < 4; i4++)
55+ w1[ib][i4] += v * 0.25;
5656 }
5757 else if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
5858 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
59- for (ii = 0; ii < 4; ii++)
60- for (jj = 0; jj < 4; jj++)
61- w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
59+ for (i4 = 0; i4 < 4; i4++)
60+ for (j4 = 0; j4 < 4; j4++)
61+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
6262 }
6363 else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
6464
6565 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
66- for (ii = 0; ii < 4; ii++)
67- for (jj = 0; jj < 4; jj++)
68- w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
66+ for (i4 = 0; i4 < 4; i4++)
67+ for (j4 = 0; j4 < 4; j4++)
68+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
6969
7070 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
71- for (ii = 0; ii < 4; ii++)
72- for (jj = 0; jj < 4; jj++)
73- w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
71+ for (i4 = 0; i4 < 4; i4++)
72+ for (j4 = 0; j4 < 4; j4++)
73+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
7474
7575 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
76- for (ii = 0; ii < 4; ii++)
77- for (jj = 0; jj < 4; jj++)
78- w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
76+ for (i4 = 0; i4 < 4; i4++)
77+ for (j4 = 0; j4 < 4; j4++)
78+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
7979 }
8080 else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
8181
8282 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
83- for (ii = 0; ii < 4; ii++)
84- for (jj = 0; jj < 4; jj++)
85- w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
83+ for (i4 = 0; i4 < 4; i4++)
84+ for (j4 = 0; j4 < 4; j4++)
85+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
8686
8787 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
88- for (ii = 0; ii < 4; ii++)
89- for (jj = 0; jj < 4; jj++)
90- w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
88+ for (i4 = 0; i4 < 4; i4++)
89+ for (j4 = 0; j4 < 4; j4++)
90+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
9191
9292 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
93- for (ii = 0; ii < 4; ii++)
94- for (jj = 0; jj < 4; jj++)
95- w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
93+ for (i4 = 0; i4 < 4; i4++)
94+ for (j4 = 0; j4 < 4; j4++)
95+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
9696 }
9797 else if (e[3] <= 0.0) {
98- for (ii = 0; ii < 4; ii++)
99- w1[ib][ii] += 0.25;
98+ for (i4 = 0; i4 < 4; i4++)
99+ w1[ib][i4] += 0.25;
100100 }
101101 }
102102 }
@@ -114,7 +114,7 @@ void dblstep(
114114 /*
115115 Main SUBROUTINE for Theta(- E1) * Theta(E1 - E2)
116116 */
117- int it, ik, ib, ii, jj, jb, **ikv, indx[4];
117+ int it, ik, ib, jb, i20, i4, j4, **ikv, indx[4];
118118 double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], *** w1, ** w2, v, tsmall[4][4], thr;
119119
120120 ikv = (int**)malloc(6 * nk * sizeof(int*));
@@ -127,9 +127,15 @@ void dblstep(
127127 ej1 = (double**)malloc(4 * sizeof(double*));
128128 ei1[0] = (double*)malloc(4 * nb * sizeof(double));
129129 ej1[0] = (double*)malloc(4 * nb * sizeof(double));
130- for (ii = 0; ii < 4; ii++) {
131- ei1[ii] = ei1[0] + ii * nb;
132- ej1[ii] = ej1[0] + ii * nb;
130+ for (i4 = 0; i4 < 4; i4++) {
131+ ei1[i4] = ei1[0] + i4 * nb;
132+ ej1[i4] = ej1[0] + i4 * nb;
133+ }
134+
135+ ej2 = (double**)malloc(nb * sizeof(double*));
136+ ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
137+ for (ib = 0; ib < nb; ib++) {
138+ ej2[ib] = ej2[0] + ib * 4;
133139 }
134140
135141 w1 = (double***)malloc(nb * sizeof(double**));
@@ -137,17 +143,11 @@ void dblstep(
137143 w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
138144 for (ib = 0; ib < nb; ib++) {
139145 w1[ib] = w1[0] + ib * 4;
140- for (ii = 0; ii < 4; ii++) {
141- w1[ib][ii] = w1[0][0] + ib * 4 * nb + ii * nb;
146+ for (i4 = 0; i4 < 4; i4++) {
147+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
142148 }
143149 }
144150
145- ej2 = (double**)malloc(nb * sizeof(double*));
146- ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
147- for (ib = 0; ib < nb; ib++) {
148- ej2[ib] = ej2[0] + ib * 4;
149- }
150-
151151 w2 = (double**)malloc(nb * sizeof(double*));
152152 w2[0] = (double*)malloc(nb * 4 * sizeof(double));
153153 for (ib = 0; ib < nb; ib++) {
@@ -165,27 +165,27 @@ void dblstep(
165165
166166 for(it = 0; it < 6*nk; it++) {
167167
168- for (ii = 0; ii < 4; ii++)
168+ for (i4 = 0; i4 < 4; i4++)
169169 for (ib = 0; ib < nb; ib++) {
170- ei1[ii][ib] = 0.0;
171- ej1[ii][ib] = 0.0;
170+ ei1[i4][ib] = 0.0;
171+ ej1[i4][ib] = 0.0;
172172 }
173- for (jj = 0; jj < 20; jj++) {
174- for (ii = 0; ii < 4; ii++) {
173+ for (i20 = 0; i20 < 20; i20++) {
174+ for (i4 = 0; i4 < 4; i4++) {
175175 for (ib = 0; ib < nb; ib++) {
176- ei1[ii][ib] += eig1[ikv[it][jj]][ib] * wlsm[jj][ii];
177- ej1[ii][ib] += eig2[ikv[it][jj]][ib] * wlsm[jj][ii];
176+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
177+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
178178 }
179179 }
180180 }
181-
181+
182182 for (ib = 0; ib < nb; ib++) {
183183
184- for (ii = 0; ii < 4; ii++)
184+ for (i4 = 0; i4 < 4; i4++)
185185 for (jb = 0; jb < nb; jb++)
186- w1[ib][ii][jb] = 0.0;
186+ w1[ib][i4][jb] = 0.0;
187187
188- for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
188+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
189189 eig_sort(4, e, indx);
190190
191191 if (e[0] <= 0.0 && 0.0 < e[1]) {
@@ -193,21 +193,21 @@ void dblstep(
193193 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
194194
195195 if (v > thr) {
196- for (jj = 0; jj < 4; jj++) {
197- ei2[jj] = 0.0;
198- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
196+ for (j4 = 0; j4 < 4; j4++) {
197+ ei2[j4] = 0.0;
198+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
199199 }
200- for (ii = 0; ii < 4; ii++)
201- for (jj = 0; jj < 4; jj++) {
202- ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
200+ for (i4 = 0; i4 < 4; i4++)
201+ for (j4 = 0; j4 < 4; j4++) {
202+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
203203 for (jb = 0; jb < nb; jb++)
204- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
204+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
205205 }
206206 libtetrabz_dblstep2(nb, ei2, ej2, w2);
207- for (ii = 0; ii < 4; ii++)
207+ for (i4 = 0; i4 < 4; i4++)
208208 for (jb = 0; jb < nb; jb++)
209- for (jj = 0; jj < 4; jj++)
210- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
209+ for (j4 = 0; j4 < 4; j4++)
210+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
211211 }
212212 }
213213 else if (e[1] <= 0.0 && 0.0 < e[2]) {
@@ -215,61 +215,61 @@ void dblstep(
215215 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
216216
217217 if (v > thr) {
218- for (jj = 0; jj < 4; jj++) {
219- ei2[jj] = 0.0;
220- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
218+ for (j4 = 0; j4 < 4; j4++) {
219+ ei2[j4] = 0.0;
220+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
221221 }
222- for (ii = 0; ii < 4; ii++)
223- for (jj = 0; jj < 4; jj++) {
224- ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
222+ for (i4 = 0; i4 < 4; i4++)
223+ for (j4 = 0; j4 < 4; j4++) {
224+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
225225 for (jb = 0; jb < nb; jb++)
226- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
226+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
227227 }
228228 libtetrabz_dblstep2(nb, ei2, ej2, w2);
229- for (ii = 0; ii < 4; ii++)
229+ for (i4 = 0; i4 < 4; i4++)
230230 for (jb = 0; jb < nb; jb++)
231- for (jj = 0; jj < 4; jj++)
232- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
231+ for (j4 = 0; j4 < 4; j4++)
232+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
233233 }
234234
235235 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
236236
237237 if (v > thr) {
238- for (jj = 0; jj < 4; jj++) {
239- ei2[jj] = 0.0;
240- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
238+ for (j4 = 0; j4 < 4; j4++) {
239+ ei2[j4] = 0.0;
240+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
241241 }
242- for (ii = 0; ii < 4; ii++)
243- for (jj = 0; jj < 4; jj++) {
244- ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
242+ for (i4 = 0; i4 < 4; i4++)
243+ for (j4 = 0; j4 < 4; j4++) {
244+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
245245 for (jb = 0; jb < nb; jb++)
246- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
246+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
247247 }
248248 libtetrabz_dblstep2(nb, ei2, ej2, w2);
249- for (ii = 0; ii < 4; ii++)
249+ for (i4 = 0; i4 < 4; i4++)
250250 for (jb = 0; jb < nb; jb++)
251- for (jj = 0; jj < 4; jj++)
252- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
251+ for (j4 = 0; j4 < 4; j4++)
252+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
253253 }
254254
255255 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
256256
257257 if (v > thr) {
258- for (jj = 0; jj < 4; jj++) {
259- ei2[jj] = 0.0;
260- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
258+ for (j4 = 0; j4 < 4; j4++) {
259+ ei2[j4] = 0.0;
260+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
261261 }
262- for (ii = 0; ii < 4; ii++)
263- for (jj = 0; jj < 4; jj++) {
264- ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
262+ for (i4 = 0; i4 < 4; i4++)
263+ for (j4 = 0; j4 < 4; j4++) {
264+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
265265 for (jb = 0; jb < nb; jb++)
266- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
266+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
267267 }
268268 libtetrabz_dblstep2(nb, ei2, ej2, w2);
269- for (ii = 0; ii < 4; ii++)
269+ for (i4 = 0; i4 < 4; i4++)
270270 for (jb = 0; jb < nb; jb++)
271- for (jj = 0; jj < 4; jj++)
272- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
271+ for (j4 = 0; j4 < 4; j4++)
272+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
273273 }
274274 }
275275 else if (e[2] <= 0.0 && 0.0 < e[3]) {
@@ -277,83 +277,83 @@ void dblstep(
277277 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
278278
279279 if (v > thr) {
280- for (jj = 0; jj < 4; jj++) {
281- ei2[jj] = 0.0;
282- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
280+ for (j4 = 0; j4 < 4; j4++) {
281+ ei2[j4] = 0.0;
282+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
283283 }
284- for (ii = 0; ii < 4; ii++)
285- for (jj = 0; jj < 4; jj++) {
286- ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
284+ for (i4 = 0; i4 < 4; i4++)
285+ for (j4 = 0; j4 < 4; j4++) {
286+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
287287 for (jb = 0; jb < nb; jb++)
288- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
288+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
289289 }
290290 libtetrabz_dblstep2(nb, ei2, ej2, w2);
291- for (ii = 0; ii < 4; ii++)
291+ for (i4 = 0; i4 < 4; i4++)
292292 for (jb = 0; jb < nb; jb++)
293- for (jj = 0; jj < 4; jj++)
294- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
293+ for (j4 = 0; j4 < 4; j4++)
294+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
295295 }
296296
297297 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
298298
299299 if (v > thr) {
300- for (jj = 0; jj < 4; jj++) {
301- ei2[jj] = 0.0;
302- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
300+ for (j4 = 0; j4 < 4; j4++) {
301+ ei2[j4] = 0.0;
302+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
303303 }
304- for (ii = 0; ii < 4; ii++)
305- for (jj = 0; jj < 4; jj++) {
306- ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
304+ for (i4 = 0; i4 < 4; i4++)
305+ for (j4 = 0; j4 < 4; j4++) {
306+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
307307 for (jb = 0; jb < nb; jb++)
308- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
308+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
309309 }
310310 libtetrabz_dblstep2(nb, ei2, ej2, w2);
311- for (ii = 0; ii < 4; ii++)
311+ for (i4 = 0; i4 < 4; i4++)
312312 for (jb = 0; jb < nb; jb++)
313- for (jj = 0; jj < 4; jj++)
314- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
313+ for (j4 = 0; j4 < 4; j4++)
314+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
315315 }
316316
317317 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
318318
319319 if (v > thr) {
320- for (jj = 0; jj < 4; jj++) {
321- ei2[jj] = 0.0;
322- for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
320+ for (j4 = 0; j4 < 4; j4++) {
321+ ei2[j4] = 0.0;
322+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
323323 }
324- for (ii = 0; ii < 4; ii++)
325- for (jj = 0; jj < 4; jj++) {
326- ei2[jj] += ei1[indx[ii]][ib] * tsmall[ii][jj];
324+ for (i4 = 0; i4 < 4; i4++)
325+ for (j4 = 0; j4 < 4; j4++) {
326+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
327327 for (jb = 0; jb < nb; jb++)
328- ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
328+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
329329 }
330330 libtetrabz_dblstep2(nb, ei2, ej2, w2);
331- for (ii = 0; ii < 4; ii++)
331+ for (i4 = 0; i4 < 4; i4++)
332332 for (jb = 0; jb < nb; jb++)
333- for (jj = 0; jj < 4; jj++)
334- w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
333+ for (j4 = 0; j4 < 4; j4++)
334+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
335335 }
336336 }
337337 else if (e[3] <= 0.0) {
338- for (ii = 0; ii < 4; ii++) {
339- ei2[ii] = ei1[ii][ib];
338+ for (i4 = 0; i4 < 4; i4++) {
339+ ei2[i4] = ei1[i4][ib];
340340 for (jb = 0; jb < nb; jb++)
341- ej2[jb][ii] = ej1[ii][jb];
341+ ej2[jb][i4] = ej1[i4][jb];
342342 }
343343 libtetrabz_dblstep2(nb, ei2, ej2, w2);
344- for (ii = 0; ii < 4; ii++)
344+ for (i4 = 0; i4 < 4; i4++)
345345 for (jb = 0; jb < nb; jb++)
346- w1[ib][ii][jb] += w2[jb][ii];
346+ w1[ib][i4][jb] += w2[jb][i4];
347347 }
348348 else {
349349 continue;
350350 }
351351 }
352- for (ii = 0; ii < 20; ii++)
353- for (jj = 0; jj < 4; jj++)
354- for (ib = 0; ib < nb; ib++)
352+ for (i20 = 0; i20 < 20; i20++)
353+ for (ib = 0; ib < nb; ib++)
354+ for (i4 = 0; i4 < 4; i4++)
355355 for (jb = 0; jb < nb; jb++)
356- wght[ikv[it][ii]][ib][jb] += w1[ib][jj][jb] * wlsm[ii][jj];
356+ wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
357357 }
358358 for (ik = 0; ik < nk; ik++)
359359 for (ib = 0; ib < nb; ib++)
--- a/src/c/libtetrabz_fermigr.c
+++ b/src/c/libtetrabz_fermigr.c
@@ -1,25 +1,25 @@
1-//
2-// Copyright (c) 2014 Mitsuaki Kawamura
3-//
4-// Permission is hereby granted, free of charge, to any person obtaining a
5-// copy of this software and associated documentation files (the
6-// "Software"), to deal in the Software without restriction, including
7-// without limitation the rights to use, copy, modify, merge, publish,
8-// distribute, sublicense, and/or sell copies of the Software, and to
9-// permit persons to whom the Software is furnished to do so, subject to
10-// the following conditions:
11-//
12-// The above copyright notice and this permission notice shall be included
13-// in all copies or substantial portions of the Software.
14-//
15-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16-// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17-// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18-// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19-// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20-// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21-// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22-//
1+/*
2+ Copyright (c) 2014 Mitsuaki Kawamura
3+
4+ Permission is hereby granted, free of charge, to any person obtaining a
5+ copy of this software and associated documentation files (the
6+ "Software"), to deal in the Software without restriction, including
7+ without limitation the rights to use, copy, modify, merge, publish,
8+ distribute, sublicense, and/or sell copies of the Software, and to
9+ permit persons to whom the Software is furnished to do so, subject to
10+ the following conditions:
11+
12+ The above copyright notice and this permission notice shall be included
13+ in all copies or substantial portions of the Software.
14+
15+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+*/
2323 #include "libtetrabz_common.h"
2424 #include <stdlib.h>
2525
@@ -29,7 +29,7 @@ void libtetrabz_fermigr3(
2929 double de[4],
3030 double **w1
3131 ) {
32- int i4, j3, ie, indx[3];
32+ int i4, j3, ie, indx[4];
3333 double e[4], tsmall[4][3], v;
3434
3535 for (i4 = 0; i4 < 4; i4++) e[i4] = de[i4];
@@ -44,26 +44,26 @@ void libtetrabz_fermigr3(
4444 libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
4545 for (i4 = 0; i4 < 4; i4++)
4646 for (j3 = 0; j3 < 3; j3++)
47- w1[i4][ie] += v * tsmall[i4][j3] / 3.0;
47+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
4848 }
4949 else if (e[1] < e0[ie] && e0[ie] <= e[2]) {
5050
5151 libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
5252 for (i4 = 0; i4 < 4; i4++)
5353 for (j3 = 0; j3 < 3; j3++)
54- w1[i4][ie] += v * tsmall[i4][j3] / 3.0;
54+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
5555
5656 libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
5757 for (i4 = 0; i4 < 4; i4++)
5858 for (j3 = 0; j3 < 3; j3++)
59- w1[i4][ie] += v * tsmall[i4][j3] / 3.0;
59+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
6060 }
6161 else if (e[2] < e0[ie] && e0[ie] < e[3]) {
6262
6363 libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
6464 for (i4 = 0; i4 < 4; i4++)
6565 for (j3 = 0; j3 < 3; j3++)
66- w1[i4][ie] += v * tsmall[i4][j3] / 3.0;
66+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
6767 }
6868 }
6969 }
@@ -98,7 +98,7 @@ void libtetrabz_fermigr2(
9898
9999 if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
100100
101- libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
101+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
102102
103103 if (v > thr) {
104104 for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
@@ -114,7 +114,7 @@ void libtetrabz_fermigr2(
114114 }
115115 else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
116116
117- libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
117+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
118118
119119 if (v > thr) {
120120 for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
@@ -128,7 +128,7 @@ void libtetrabz_fermigr2(
128128 w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
129129 }
130130
131- libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
131+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
132132
133133 if (v > thr) {
134134 for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
@@ -142,7 +142,7 @@ void libtetrabz_fermigr2(
142142 w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
143143 }
144144
145- libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
145+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
146146
147147 if (v > thr) {
148148 for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
@@ -158,8 +158,8 @@ void libtetrabz_fermigr2(
158158 }
159159 else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
160160
161- libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
162-
161+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
162+
163163 if (v > thr) {
164164 for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
165165 for (i4 = 0; i4 < 4; i4++)
@@ -172,7 +172,7 @@ void libtetrabz_fermigr2(
172172 w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
173173 }
174174
175- libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
175+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
176176
177177 if (v > thr) {
178178 for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
@@ -186,7 +186,7 @@ void libtetrabz_fermigr2(
186186 w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
187187 }
188188
189- libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
189+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
190190
191191 if (v > thr) {
192192 for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
@@ -201,9 +201,8 @@ void libtetrabz_fermigr2(
201201 }
202202 }
203203 else if (e[3] <= 0.0) {
204- for (i4 = 0; i4 < 4; i4++) {
204+ for (i4 = 0; i4 < 4; i4++)
205205 de[i4] = ej1[ib][i4] - ei1[i4];
206- }
207206 libtetrabz_fermigr3(ne, e0, de, w2);
208207 for (i4 = 0; i4 < 4; i4++)
209208 for (ie = 0; ie < ne; ie++)
@@ -287,7 +286,7 @@ void fermigr(
287286
288287 thr = 1.0e-10;
289288
290- for(it = 0; it < 6*nk; it++) {
289+ for (it = 0; it < 6*nk; it++) {
291290
292291 for (i4 = 0; i4 < 4; i4++)
293292 for (ib = 0; ib < nb; ib++) {
@@ -302,13 +301,13 @@ void fermigr(
302301 }
303302 }
304303 }
305-
304+
306305 for (ib = 0; ib < nb; ib++) {
307306
308- for (jb = 0; jb < nb; jb++)
309- for (i4 = 0; i4 < 4; i4++)
307+ for (i4 = 0; i4 < 4; i4++)
308+ for (jb = 0; jb < nb; jb++)
310309 for (ie = 0; ie < ne; ie++)
311- w1[ib][jb][ie][i4] = 0.0;
310+ w1[ib][i4][jb][ie] = 0.0;
312311
313312 for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
314313 eig_sort(4, e, indx);
@@ -324,7 +323,7 @@ void fermigr(
324323 }
325324 for (i4 = 0; i4 < 4; i4++)
326325 for (j4 = 0; j4 < 4; j4++) {
327- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
326+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
328327 for (jb = 0; jb < nb; jb++)
329328 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
330329 }
@@ -347,7 +346,7 @@ void fermigr(
347346 }
348347 for (i4 = 0; i4 < 4; i4++)
349348 for (j4 = 0; j4 < 4; j4++) {
350- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
349+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
351350 for (jb = 0; jb < nb; jb++)
352351 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
353352 }
@@ -368,7 +367,7 @@ void fermigr(
368367 }
369368 for (i4 = 0; i4 < 4; i4++)
370369 for (j4 = 0; j4 < 4; j4++) {
371- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
370+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
372371 for (jb = 0; jb < nb; jb++)
373372 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
374373 }
@@ -389,7 +388,7 @@ void fermigr(
389388 }
390389 for (i4 = 0; i4 < 4; i4++)
391390 for (j4 = 0; j4 < 4; j4++) {
392- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
391+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
393392 for (jb = 0; jb < nb; jb++)
394393 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
395394 }
@@ -412,7 +411,7 @@ void fermigr(
412411 }
413412 for (i4 = 0; i4 < 4; i4++)
414413 for (j4 = 0; j4 < 4; j4++) {
415- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
414+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
416415 for (jb = 0; jb < nb; jb++)
417416 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
418417 }
@@ -433,7 +432,7 @@ void fermigr(
433432 }
434433 for (i4 = 0; i4 < 4; i4++)
435434 for (j4 = 0; j4 < 4; j4++) {
436- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
435+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
437436 for (jb = 0; jb < nb; jb++)
438437 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
439438 }
@@ -454,7 +453,7 @@ void fermigr(
454453 }
455454 for (i4 = 0; i4 < 4; i4++)
456455 for (j4 = 0; j4 < 4; j4++) {
457- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
456+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
458457 for (jb = 0; jb < nb; jb++)
459458 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
460459 }
@@ -467,12 +466,11 @@ void fermigr(
467466 }
468467 }
469468 else if (e[3] <= 0.0) {
470-
471- for (i4 = 0; i4 < 4; i4++) {
472- ei2[i4] = ej1[i4][ib];
473- for (jb = 0; jb < nb; jb++)
474- ej2[jb][i4] = ej1[i4][jb];
475- }
469+ for (i4 = 0; i4 < 4; i4++) {
470+ ei2[i4] = ei1[i4][ib];
471+ for (jb = 0; jb < nb; jb++)
472+ ej2[jb][i4] = ej1[i4][jb];
473+ }
476474 libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
477475 for (i4 = 0; i4 < 4; i4++)
478476 for (jb = 0; jb < nb; jb++)
@@ -483,16 +481,32 @@ void fermigr(
483481 continue;
484482 }
485483 }
486- for (i4 = 0; i4 < 20; i4++)
484+ for (i20 = 0; i20 < 20; i20++)
487485 for (ib = 0; ib < nb; ib++)
488- for (j4 = 0; j4 < 4; j4++)
486+ for (i4 = 0; i4 < 4; i4++)
489487 for (jb = 0; jb < nb; jb++)
490488 for (ie = 0; ie < ne; ie++)
491- wght[ikv[it][i20]][ib][jb][ie] += wlsm[i20][j4] * w1[ib][j4][jb][ie];
489+ wght[ikv[it][i20]][ib][jb][ie] += wlsm[i20][i4] * w1[ib][i4][jb][ie];
492490 }
493491 for (ik = 0; ik < nk; ik++)
494492 for (ib = 0; ib < nb; ib++)
495493 for (jb = 0; jb < nb; jb++)
496494 for (ie = 0; ie < ne; ie++)
497495 wght[ik][ib][jb][ie] /= (6.0 * (double) nk);
496+
497+ free(ikv[0]);
498+ free(ikv);
499+ free(ei1[0]);
500+ free(ei1);
501+ free(ej1[0]);
502+ free(ej1);
503+ free(ej2[0]);
504+ free(ej2);
505+ free(w1[0][0][0]);
506+ free(w1[0][0]);
507+ free(w1[0]);
508+ free(w1);
509+ free(w2[0][0]);
510+ free(w2[0]);
511+ free(w2);
498512 }
--- a/src/c/libtetrabz_polcmplx.c
+++ b/src/c/libtetrabz_polcmplx.c
@@ -257,7 +257,7 @@ void libtetrabz_polcmplx3(
257257 Tetrahedron method for delta(om - ep + e)
258258 */
259259 int i4, ir, indx[4], ie;
260- double e[4], x[4], thr, w2[4][2];
260+ double e[4], x[4], thr, w2[4][2], denom;
261261
262262 for (i4 = 0; i4 < 3; i4++) e[i4] = de[i4];
263263 eig_sort(4, e, indx);
@@ -268,10 +268,11 @@ void libtetrabz_polcmplx3(
268268 The former is more stable.
269269 The latter is more accurate ?
270270 */
271- for (i4 = 0; i4 < 4; i4++)
272- for (ir = 0; ir < 2; ir++) {
273- w1[i4][ie][ir] = 0.25 / (de[i4] + e0[ie][ir]);
274- }
271+ for (i4 = 0; i4 < 4; i4++) {
272+ denom = (de[i4] + e0[ie][0]) * (de[i4] + e0[ie][0]) + e0[ie][1] * e0[ie][1];
273+ w1[i4][ie][0] = 0.25 * (de[i4] + e0[ie][0]) / denom;
274+ w1[i4][ie][1] = 0.25 * (-e0[ie][1]) / denom;
275+ }
275276 continue;
276277
277278 for (i4 = 0; i4 < 4; i4++)
@@ -572,9 +573,8 @@ void libtetrabz_polcmplx2(
572573 }
573574 }
574575 else if (e[3] <= 0.0) {
575- for (i4 = 0; i4 < 4; i4++) {
576+ for (i4 = 0; i4 < 4; i4++)
576577 de[i4] = ej1[ib][i4] - ei1[i4];
577- }
578578 libtetrabz_polcmplx3(ne, e0, de, w2);
579579 for (i4 = 0; i4 < 4; i4++)
580580 for (ie = 0; ie < ne; ie++)
@@ -621,7 +621,7 @@ void polcmplx(
621621
622622 ej2 = (double**)malloc(nb * sizeof(double*));
623623 ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
624- for (ib = 0; ib < nb; ib++)
624+ for (ib = 0; ib < nb; ib++)
625625 ej2[ib] = ej2[0] + ib * 4;
626626
627627 w1 = (double*****)malloc(nb * sizeof(double****));
@@ -685,11 +685,11 @@ void polcmplx(
685685
686686 for (ib = 0; ib < nb; ib++) {
687687
688- for (jb = 0; jb < nb; jb++)
689- for (i4 = 0; i4 < 4; i4++)
688+ for (i4 = 0; i4 < 4; i4++)
689+ for (jb = 0; jb < nb; jb++)
690690 for (ie = 0; ie < ne; ie++)
691691 for (ir = 0; ir < 2; ir++)
692- w1[ib][jb][ie][i4][ir] = 0.0;
692+ w1[ib][i4][jb][ie][ir] = 0.0;
693693
694694 for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
695695 eig_sort(4, e, indx);
@@ -705,7 +705,7 @@ void polcmplx(
705705 }
706706 for (i4 = 0; i4 < 4; i4++)
707707 for (j4 = 0; j4 < 4; j4++) {
708- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
708+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
709709 for (jb = 0; jb < nb; jb++)
710710 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
711711 }
@@ -729,7 +729,7 @@ void polcmplx(
729729 }
730730 for (i4 = 0; i4 < 4; i4++)
731731 for (j4 = 0; j4 < 4; j4++) {
732- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
732+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
733733 for (jb = 0; jb < nb; jb++)
734734 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
735735 }
@@ -751,7 +751,7 @@ void polcmplx(
751751 }
752752 for (i4 = 0; i4 < 4; i4++)
753753 for (j4 = 0; j4 < 4; j4++) {
754- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
754+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
755755 for (jb = 0; jb < nb; jb++)
756756 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
757757 }
@@ -773,7 +773,7 @@ void polcmplx(
773773 }
774774 for (i4 = 0; i4 < 4; i4++)
775775 for (j4 = 0; j4 < 4; j4++) {
776- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
776+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
777777 for (jb = 0; jb < nb; jb++)
778778 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
779779 }
@@ -797,7 +797,7 @@ void polcmplx(
797797 }
798798 for (i4 = 0; i4 < 4; i4++)
799799 for (j4 = 0; j4 < 4; j4++) {
800- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
800+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
801801 for (jb = 0; jb < nb; jb++)
802802 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
803803 }
@@ -819,7 +819,7 @@ void polcmplx(
819819 }
820820 for (i4 = 0; i4 < 4; i4++)
821821 for (j4 = 0; j4 < 4; j4++) {
822- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
822+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
823823 for (jb = 0; jb < nb; jb++)
824824 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
825825 }
@@ -841,7 +841,7 @@ void polcmplx(
841841 }
842842 for (i4 = 0; i4 < 4; i4++)
843843 for (j4 = 0; j4 < 4; j4++) {
844- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
844+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
845845 for (jb = 0; jb < nb; jb++)
846846 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
847847 }
@@ -853,11 +853,10 @@ void polcmplx(
853853 for (ir = 0; ir < 2; ir++)
854854 w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
855855 }
856- }
856+ }
857857 else if (e[3] <= 0.0) {
858-
859858 for (i4 = 0; i4 < 4; i4++) {
860- ei2[i4] = ej1[i4][ib];
859+ ei2[i4] = ei1[i4][ib];
861860 for (jb = 0; jb < nb; jb++)
862861 ej2[jb][i4] = ej1[i4][jb];
863862 }
@@ -872,13 +871,13 @@ void polcmplx(
872871 continue;
873872 }
874873 }
875- for (i4 = 0; i4 < 20; i4++)
874+ for (i20 = 0; i20 < 20; i20++)
876875 for (ib = 0; ib < nb; ib++)
877- for (j4 = 0; j4 < 4; j4++)
876+ for (i4 = 0; i4 < 4; i4++)
878877 for (jb = 0; jb < nb; jb++)
879878 for (ie = 0; ie < ne; ie++)
880879 for (ir = 0; ir < 2; ir++)
881- wght[ikv[it][i20]][ib][jb][ie][ir] += wlsm[i20][j4] * w1[ib][j4][jb][ie][ir];
880+ wght[ikv[it][i20]][ib][jb][ie][ir] += wlsm[i20][i4] * w1[ib][i4][jb][ie][ir];
882881 }
883882 for (ik = 0; ik < nk; ik++)
884883 for (ib = 0; ib < nb; ib++)
--- a/src/c/libtetrabz_polstat.c
+++ b/src/c/libtetrabz_polstat.c
@@ -1,25 +1,25 @@
1-//
2-// Copyright (c) 2014 Mitsuaki Kawamura
3-//
4-// Permission is hereby granted, free of charge, to any person obtaining a
5-// copy of this software and associated documentation files (the
6-// "Software"), to deal in the Software without restriction, including
7-// without limitation the rights to use, copy, modify, merge, publish,
8-// distribute, sublicense, and/or sell copies of the Software, and to
9-// permit persons to whom the Software is furnished to do so, subject to
10-// the following conditions:
11-//
12-// The above copyright notice and this permission notice shall be included
13-// in all copies or substantial portions of the Software.
14-//
15-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16-// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17-// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18-// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19-// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20-// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21-// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22-//
1+/*
2+ Copyright (c) 2014 Mitsuaki Kawamura
3+
4+ Permission is hereby granted, free of charge, to any person obtaining a
5+ copy of this software and associated documentation files (the
6+ "Software"), to deal in the Software without restriction, including
7+ without limitation the rights to use, copy, modify, merge, publish,
8+ distribute, sublicense, and/or sell copies of the Software, and to
9+ permit persons to whom the Software is furnished to do so, subject to
10+ the following conditions:
11+
12+ The above copyright notice and this permission notice shall be included
13+ in all copies or substantial portions of the Software.
14+
15+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19+ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20+ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21+ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22+*/
2323 #include "libtetrabz_common.h"
2424 #include <math.h>
2525 #include <stdio.h>
@@ -149,8 +149,8 @@ double libtetrabz_polstat_1211(
149149 }
150150
151151 void libtetrabz_polstat3(
152- double *de,
153- double *w1
152+ double de[4],
153+ double w1[4]
154154 )
155155 {
156156 /*
@@ -168,7 +168,7 @@ void libtetrabz_polstat3(
168168 for(i4 =0; i4 <4; i4++){
169169 if (e[i4] < thr2) {
170170 if (i4 == 3) {
171- printf(" Nesting # ");
171+ printf(" Nesting # \n");
172172 }
173173 ln[i4] = 0.0;
174174 e[i4] = 0.0;
@@ -200,7 +200,7 @@ void libtetrabz_polstat3(
200200 if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
201201 printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
202202 printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
203- printf("weighting 4=3=2");
203+ printf("weighting 4=3=2\n");
204204 }
205205 }
206206 }
@@ -216,7 +216,7 @@ void libtetrabz_polstat3(
216216 if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
217217 printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
218218 printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
219- printf("weighting 4=3 2=1");
219+ printf("weighting 4=3 2=1\n");
220220 }
221221 }
222222 else {
@@ -231,7 +231,7 @@ void libtetrabz_polstat3(
231231 if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
232232 printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
233233 printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
234- printf("weighting 4=3");
234+ printf("weighting 4=3\n");
235235 }
236236 }
237237 }
@@ -248,7 +248,7 @@ void libtetrabz_polstat3(
248248 if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
249249 printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
250250 printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
251- printf("weighting 3=2=1");
251+ printf("weighting 3=2=1\n");
252252 }
253253 }
254254 else {
@@ -263,7 +263,7 @@ void libtetrabz_polstat3(
263263 if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
264264 printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
265265 printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
266- printf("weighting 3=2");
266+ printf("weighting 3=2\n");
267267 }
268268 }
269269 }
@@ -279,7 +279,7 @@ void libtetrabz_polstat3(
279279 if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
280280 printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
281281 printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
282- printf("weighting 2=1");
282+ printf("weighting 2=1\n");
283283 }
284284 }
285285 else {
@@ -294,7 +294,7 @@ void libtetrabz_polstat3(
294294 if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
295295 printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
296296 printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
297- printf("weighting");
297+ printf("weighting\n");
298298 }
299299 }
300300 }
@@ -318,6 +318,9 @@ void libtetrabz_polstat2(
318318
319319 for (ib = 0; ib < nb; ib++) {
320320
321+ for (i4 = 0; i4 < 4; i4++)
322+ w1[ib][i4] = 0.0;
323+
321324 for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
322325 eig_sort(4, e, indx);
323326
@@ -326,7 +329,7 @@ void libtetrabz_polstat2(
326329 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
327330
328331 if (v > thr) {
329- for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
332+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
330333 for (i4 = 0; i4 < 4; i4++)
331334 for (j4 = 0; j4 < 4; j4++)
332335 de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
@@ -341,7 +344,7 @@ void libtetrabz_polstat2(
341344 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
342345
343346 if (v > thr) {
344- for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
347+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
345348 for (i4 = 0; i4 < 4; i4++)
346349 for (j4 = 0; j4 < 4; j4++)
347350 de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
@@ -354,7 +357,7 @@ void libtetrabz_polstat2(
354357 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
355358
356359 if (v > thr) {
357- for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
360+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
358361 for (i4 = 0; i4 < 4; i4++)
359362 for (j4 = 0; j4 < 4; j4++)
360363 de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
@@ -367,7 +370,7 @@ void libtetrabz_polstat2(
367370 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
368371
369372 if (v > thr) {
370- for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
373+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
371374 for (i4 = 0; i4 < 4; i4++)
372375 for (j4 = 0; j4 < 4; j4++)
373376 de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
@@ -380,9 +383,9 @@ void libtetrabz_polstat2(
380383 else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
381384
382385 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
383-
386+
384387 if (v > thr) {
385- for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
388+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
386389 for (i4 = 0; i4 < 4; i4++)
387390 for (j4 = 0; j4 < 4; j4++)
388391 de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
@@ -395,7 +398,7 @@ void libtetrabz_polstat2(
395398 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
396399
397400 if (v > thr) {
398- for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
401+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
399402 for (i4 = 0; i4 < 4; i4++)
400403 for (j4 = 0; j4 < 4; j4++)
401404 de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
@@ -408,7 +411,7 @@ void libtetrabz_polstat2(
408411 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
409412
410413 if (v > thr) {
411- for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
414+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
412415 for (i4 = 0; i4 < 4; i4++)
413416 for (j4 = 0; j4 < 4; j4++)
414417 de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
@@ -422,7 +425,8 @@ void libtetrabz_polstat2(
422425 for (i4 = 0; i4 < 4; i4++)
423426 de[i4] = ej1[ib][i4] - ei1[i4];
424427 libtetrabz_polstat3(de, w2);
425- for (i4 = 0; i4 < 4; i4++) w1[ib][i4] += w2[i4];
428+ for (i4 = 0; i4 < 4; i4++)
429+ w1[ib][i4] += w2[i4];
426430 }
427431 }
428432 }
@@ -447,8 +451,8 @@ void polstat(
447451
448452 ikv = (int**)malloc(6 * nk * sizeof(int*));
449453 ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
450- for (ik = 0; ik < 6 * nk; ik++) {
451- ikv[ik] = ikv[0] + ik * 20;
454+ for (it = 0; it < 6 * nk; it++) {
455+ ikv[it] = ikv[0] + it * 20;
452456 }
453457
454458 ei1 = (double**)malloc(4 * sizeof(double*));
@@ -462,9 +466,8 @@ void polstat(
462466
463467 ej2 = (double**)malloc(nb * sizeof(double*));
464468 ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
465- for (ib = 0; ib < nb; ib++) {
469+ for (ib = 0; ib < nb; ib++)
466470 ej2[ib] = ej2[0] + ib * 4;
467- }
468471
469472 w1 = (double***)malloc(nb * sizeof(double**));
470473 w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
@@ -504,7 +507,7 @@ void polstat(
504507 }
505508 }
506509 }
507-
510+
508511 for (ib = 0; ib < nb; ib++) {
509512
510513 for (i4 = 0; i4 < 4; i4++)
@@ -525,7 +528,7 @@ void polstat(
525528 }
526529 for (i4 = 0; i4 < 4; i4++)
527530 for (j4 = 0; j4 < 4; j4++) {
528- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
531+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
529532 for (jb = 0; jb < nb; jb++)
530533 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
531534 }
@@ -547,7 +550,7 @@ void polstat(
547550 }
548551 for (i4 = 0; i4 < 4; i4++)
549552 for (j4 = 0; j4 < 4; j4++) {
550- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
553+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
551554 for (jb = 0; jb < nb; jb++)
552555 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
553556 }
@@ -567,7 +570,7 @@ void polstat(
567570 }
568571 for (i4 = 0; i4 < 4; i4++)
569572 for (j4 = 0; j4 < 4; j4++) {
570- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
573+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
571574 for (jb = 0; jb < nb; jb++)
572575 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
573576 }
@@ -587,7 +590,7 @@ void polstat(
587590 }
588591 for (i4 = 0; i4 < 4; i4++)
589592 for (j4 = 0; j4 < 4; j4++) {
590- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
593+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
591594 for (jb = 0; jb < nb; jb++)
592595 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
593596 }
@@ -609,7 +612,7 @@ void polstat(
609612 }
610613 for (i4 = 0; i4 < 4; i4++)
611614 for (j4 = 0; j4 < 4; j4++) {
612- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
615+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
613616 for (jb = 0; jb < nb; jb++)
614617 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
615618 }
@@ -629,7 +632,7 @@ void polstat(
629632 }
630633 for (i4 = 0; i4 < 4; i4++)
631634 for (j4 = 0; j4 < 4; j4++) {
632- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
635+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
633636 for (jb = 0; jb < nb; jb++)
634637 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
635638 }
@@ -649,7 +652,7 @@ void polstat(
649652 }
650653 for (i4 = 0; i4 < 4; i4++)
651654 for (j4 = 0; j4 < 4; j4++) {
652- ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
655+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
653656 for (jb = 0; jb < nb; jb++)
654657 ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
655658 }
@@ -659,10 +662,10 @@ void polstat(
659662 for (j4 = 0; j4 < 4; j4++)
660663 w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
661664 }
662- }
665+ }
663666 else if (e[3] <= 0.0) {
664667 for (i4 = 0; i4 < 4; i4++) {
665- ei2[i4] = ej1[i4][ib];
668+ ei2[i4] = ei1[i4][ib];
666669 for (jb = 0; jb < nb; jb++)
667670 ej2[jb][i4] = ej1[i4][jb];
668671 }
@@ -698,5 +701,5 @@ void polstat(
698701 free(w1[0]);
699702 free(w1);
700703 free(w2[0]);
701- free(w2);
704+ free(w2);
702705 }
--- a/src/c/test.c
+++ b/src/c/test.c
@@ -391,14 +391,14 @@ void test_polstat(
391391 for (ib = 0; ib < nb; ib++)
392392 for (jb = 0; jb < nb; jb++)
393393 val[ib][jb] += wght[ik][ib][jb] * mat[ik];
394- printf("# libtetrabz_polstat");
395- printf(" %15.5e %15.5e", pi * (68.0 + 45.0 * log(3.0)) / 96.0, val[0][0] * vbz);
396- printf(" %15.5e %15.5e", pi * 8.0 / 5.0, val[0][1] * vbz);
397- printf(" %15.5e %15.5e", pi * (228.0 + 22.0 * sqrt(2.0) - 96.0 * log(2.0)
394+ printf("# libtetrabz_polstat\n");
395+ printf(" %15.5e %15.5e\n", pi * (68.0 + 45.0 * log(3.0)) / 96.0, val[0][0] * vbz);
396+ printf(" %15.5e %15.5e\n", pi * 8.0 / 5.0, val[0][1] * vbz);
397+ printf(" %15.5e %15.5e\n", pi * (228.0 + 22.0 * sqrt(2.0) - 96.0 * log(2.0)
398398 + 192.0 * log(4.0 + sqrt(2.0))
399399 - 3.0 * log(1.0 + 2.0 * sqrt(2.0))) / 1536.0,
400400 val[1][0] * vbz);
401- printf(" %15.5e %15.5e", pi * sqrt(8.0) / 5.0, val[1][1] * vbz);
401+ printf(" %15.5e %15.5e\n", pi * sqrt(8.0) / 5.0, val[1][1] * vbz);
402402 printf("\n");
403403
404404 free(wght[0][0]);
@@ -564,8 +564,10 @@ void test_polcmplx(
564564
565565 for (ie = 0; ie < ne; ie++) {
566566 for (ir = 0; ir < 2; ir++) {
567- val0[0][1][ie][ir] = (8.0 * pi) / (5.0 * (1.0 + 2.0 * e0[ie][ir]));
568- val0[1][1][ie][ir] = (sqrt(8.0) * pi) / (5.0 * (1.0 + 4.0 * e0[ie][ir]));
567+ val0[0][1][ie][0] = (8.0 * pi * (1.0 + 2.0 * e0[ie][0])) / (5.0 * (pow(1.0 + 2.0 * e0[ie][0], 2) + pow(2.0 * e0[ie][1], 2)));
568+ val0[0][1][ie][1] = (- 8.0 * pi * 2.0 * e0[ie][1]) / (5.0 * (pow(1.0 + 2.0 * e0[ie][0], 2) + pow(2.0 * e0[ie][1], 2)));
569+ val0[1][1][ie][0] = (sqrt(8.0) * pi * (1.0 + 4.0 * e0[ie][0])) / (5.0 * (pow(1.0 + 4.0 * e0[ie][0],2) + pow(4.0 * e0[ie][1], 2)));
570+ val0[1][1][ie][1] = (- sqrt(8.0) * pi * 4.0 * e0[ie][1]) / (5.0 * (pow(1.0 + 4.0 * e0[ie][0], 2) + pow(4.0 * e0[ie][1], 2)));
569571 }
570572 }
571573
@@ -669,7 +671,7 @@ int main()
669671 }
670672
671673 printf("# Ideal Result\n");
672- /*
674+
673675 test_occ(ng, nk, nb, bvec, vbz, eig1, mat);
674676
675677 test_fermieng(ng, nk, nb, bvec, vbz, eig1, mat);
@@ -677,15 +679,14 @@ int main()
677679 test_dos(ng, nk, nb, bvec, vbz, eig1, mat);
678680
679681 test_intdos(ng, nk, nb, bvec, vbz, eig1, mat);
680- */
682+
681683 test_dblstep(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
682- /*
684+
683685 test_dbldelta(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
684-
686+
685687 test_polstat(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
686-
688+
687689 test_fermigr(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
688690
689691 test_polcmplx(ng, nk, nb, bvec, vbz, eig1, eig2, mat);
690- */
691692 }
Show on old repository browser