• R/O
  • HTTP
  • SSH
  • HTTPS

python: Commit

libtetrabz python package


Commit MetaInfo

Revisionf163836f2012ce195962046ad5aa99ab04ce2c79 (tree)
Zeit2022-03-16 02:06:38
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Backup

Ändern Zusammenfassung

Diff

--- a/src/c/libtetrabz_common.c
+++ b/src/c/libtetrabz_common.c
@@ -33,27 +33,27 @@ void libtetrabz_initialize(
3333 define shortest diagonal line & define type of tetragonal
3434 */
3535 int ii, i0, i1, i2, i3, it, itype, ivvec0[4], divvec[4][4], ivvec[6][20][3], ikv0[3], nt;
36- double bvec2[3][3], bvec3[4][3], *bnorm;
36+ double bvec2[3][3], bvec3[4][3], bnorm[4];
3737 double wlsm1[4][4] = { {1440.0, 0.0, 30.0, 0.0},
38- {0.0, 1440.0, 0.0, 30.0},
39- {30.0, 0.0, 1440.0, 0.0},
40- {0.0, 30.0, 0.0, 1440.0} },
38+ {0.0, 1440.0, 0.0, 30.0},
39+ {30.0, 0.0, 1440.0, 0.0},
40+ {0.0, 30.0, 0.0, 1440.0} },
4141 wlsm2[4][4] = {{-38.0, 7.0, 17.0, -28.0},
42- {-28.0, -38.0, 7.0, 17.0},
43- {17.0, -28.0, -38.0, 7.0},
44- {7.0, 17.0, -28.0, -38.0}},
42+ {-28.0, -38.0, 7.0, 17.0},
43+ {17.0, -28.0, -38.0, 7.0},
44+ {7.0, 17.0, -28.0, -38.0}},
4545 wlsm3[4][4] = {{-56.0, 9.0, -46.0, 9.0},
46- {9.0, -56.0, 9.0, -46.0},
47- {-46.0, 9.0, -56.0, 9.0},
48- {9.0, -46.0, 9.0, -56.0}},
46+ {9.0, -56.0, 9.0, -46.0},
47+ {-46.0, 9.0, -56.0, 9.0},
48+ {9.0, -46.0, 9.0, -56.0}},
4949 wlsm4[4][4] = {{-38.0, -28.0, 17.0, 7.0},
50- {7.0, -38.0, -28.0, 17.0},
51- {17.0, 7.0, -38.0, -28.0},
52- {-28.0, 17.0, 7.0, -38.0}},
50+ {7.0, -38.0, -28.0, 17.0},
51+ {17.0, 7.0, -38.0, -28.0},
52+ {-28.0, 17.0, 7.0, -38.0}},
5353 wlsm5[4][4] = {{-18.0, -18.0, 12.0, -18.0},
54- {-18.0, -18.0, -18.0, 12.0},
55- {12.0, -18.0, -18.0, -18.0},
56- {-18.0, 12.0, -18.0, -18.0}};
54+ {-18.0, -18.0, -18.0, 12.0},
55+ {12.0, -18.0, -18.0, -18.0},
56+ {-18.0, 12.0, -18.0, -18.0}};
5757
5858 for (i1 = 0; i1 < 3; i1++)
5959 for (i2 = 0; i2 < 3; i2++)
@@ -136,11 +136,11 @@ void libtetrabz_initialize(
136136
137137 for (i1 = 0; i1 < 4; i1++) {
138138 for (i2 = 0; i2 < 4; i2++) {
139- wlsm[i1][i2] = wlsm1[i1][i2] /= 1260.0;
140- wlsm[i1][i2 + 4] = wlsm2[i1][i2] /= 1260.0;
141- wlsm[i1][i2 + 8] = wlsm3[i1][i2] /= 1260.0;
142- wlsm[i1][i2 + 12] = wlsm4[i1][i2] /= 1260.0;
143- wlsm[i1][i2 + 16] = wlsm5[i1][i2] /= 1260.0;
139+ wlsm[i2][i1] = wlsm1[i1][i2] /= 1260.0;
140+ wlsm[i2 + 4][i1] = wlsm2[i1][i2] /= 1260.0;
141+ wlsm[i2 + 8][i1] = wlsm3[i1][i2] /= 1260.0;
142+ wlsm[i2 + 12][i1] = wlsm4[i1][i2] /= 1260.0;
143+ wlsm[i2 + 16][i1] = wlsm5[i1][i2] /= 1260.0;
144144 }
145145 }
146146 /*
@@ -184,18 +184,18 @@ void libtetrabz_tsmall_a1(
184184 *v = a10 * a20 * a30;
185185
186186 tsmall[0][0] = 1.0;
187- tsmall[0][1] = 0.0;
188- tsmall[0][2] = 0.0;
189- tsmall[0][3] = 0.0;
190- tsmall[1][0] = 1.0 - a10;
187+ tsmall[0][1] = 1.0 - a10;
188+ tsmall[0][2] = 1.0 - a20;
189+ tsmall[0][3] = 1.0 - a30;
190+ tsmall[1][0] = 0.0;
191191 tsmall[1][1] = a10;
192192 tsmall[1][2] = 0.0;
193193 tsmall[1][3] = 0.0;
194- tsmall[2][0] = 1.0 - a20;
194+ tsmall[2][0] = 0.0;
195195 tsmall[2][1] = 0.0;
196196 tsmall[2][2] = a20;
197197 tsmall[2][3] = 0.0;
198- tsmall[3][0] = 1.0 - a30;
198+ tsmall[3][0] = 0.0;
199199 tsmall[3][1] = 0.0;
200200 tsmall[3][2] = 0.0;
201201 tsmall[3][3] = a30;
@@ -222,20 +222,20 @@ void libtetrabz_tsmall_b1(
222222 *v = a20 * a30 * a13;
223223
224224 tsmall[0][0] = 1.0;
225- tsmall[0][1] = 0.0;
226- tsmall[0][2] = 0.0;
225+ tsmall[0][1] = 1.0 - a20;
226+ tsmall[0][2] = 1.0 - a30;
227227 tsmall[0][3] = 0.0;
228- tsmall[1][0] = 1.0 - a20;
228+ tsmall[1][0] = 0.0;
229229 tsmall[1][1] = 0.0;
230- tsmall[1][2] = a20;
231- tsmall[1][3] = 0.0;
232- tsmall[2][0] = 1.0 - a30;
233- tsmall[2][1] = 0.0;
230+ tsmall[1][2] = 0.0;
231+ tsmall[1][3] = a13;
232+ tsmall[2][0] = 0.0;
233+ tsmall[2][1] = a20;
234234 tsmall[2][2] = 0.0;
235- tsmall[2][3] = a30;
235+ tsmall[2][3] = 0.0;
236236 tsmall[3][0] = 0.0;
237- tsmall[3][1] = a13;
238- tsmall[3][2] = 0.0;
237+ tsmall[3][1] = 0.0;
238+ tsmall[3][2] = a30;
239239 tsmall[3][3] = 1.0 - a13;
240240 }
241241
@@ -264,14 +264,14 @@ void libtetrabz_tsmall_b2(
264264 tsmall[0][3] = 0.0;
265265 tsmall[1][0] = 0.0;
266266 tsmall[1][1] = 1.0;
267- tsmall[1][2] = 0.0;
268- tsmall[1][3] = 0.0;
267+ tsmall[1][2] = 1.0 - a21;
268+ tsmall[1][3] = 1.0 - a31;
269269 tsmall[2][0] = 0.0;
270- tsmall[2][1] = 1.0 - a21;
270+ tsmall[2][1] = 0.0;
271271 tsmall[2][2] = a21;
272272 tsmall[2][3] = 0.0;
273273 tsmall[3][0] = 0.0;
274- tsmall[3][1] = 1.0 - a31;
274+ tsmall[3][1] = 0.0;
275275 tsmall[3][2] = 0.0;
276276 tsmall[3][3] = a31;
277277 }
@@ -296,19 +296,19 @@ void libtetrabz_tsmall_b3(
296296 *v = a12 * a20 * a31;
297297
298298 tsmall[0][0] = 1.0;
299- tsmall[0][1] = 0.0;
299+ tsmall[0][1] = 1.0 - a20;
300300 tsmall[0][2] = 0.0;
301301 tsmall[0][3] = 0.0;
302- tsmall[1][0] = 1.0 - a20;
302+ tsmall[1][0] = 0.0;
303303 tsmall[1][1] = 0.0;
304- tsmall[1][2] = a20;
305- tsmall[1][3] = 0.0;
304+ tsmall[1][2] = a12;
305+ tsmall[1][3] = 1.0 - a31;
306306 tsmall[2][0] = 0.0;
307- tsmall[2][1] = a12;
307+ tsmall[2][1] = a20;
308308 tsmall[2][2] = 1.0 - a12;
309309 tsmall[2][3] = 0.0;
310310 tsmall[3][0] = 0.0;
311- tsmall[3][1] = 1.0 - a31;
311+ tsmall[3][1] = 0.0;
312312 tsmall[3][2] = 0.0;
313313 tsmall[3][3] = a31;
314314 }
@@ -342,10 +342,10 @@ void libtetrabz_tsmall_c1(
342342 tsmall[2][0] = 0.0;
343343 tsmall[2][1] = 0.0;
344344 tsmall[2][2] = 1.0;
345- tsmall[2][3] = 0.0;
345+ tsmall[2][3] = 1.0 - a32;
346346 tsmall[3][0] = 0.0;
347347 tsmall[3][1] = 0.0;
348- tsmall[3][2] = 1.0 - a32;
348+ tsmall[3][2] = 0.0;
349349 tsmall[3][3] = a32;
350350 }
351351
@@ -374,15 +374,15 @@ void libtetrabz_tsmall_c2(
374374 tsmall[0][3] = 0.0;
375375 tsmall[1][0] = 0.0;
376376 tsmall[1][1] = 1.0;
377- tsmall[1][2] = 0.0;
377+ tsmall[1][2] = 1.0 - a31;
378378 tsmall[1][3] = 0.0;
379379 tsmall[2][0] = 0.0;
380- tsmall[2][1] = 1.0 - a31;
380+ tsmall[2][1] = 0.0;
381381 tsmall[2][2] = 0.0;
382- tsmall[2][3] = a31;
382+ tsmall[2][3] = a23;
383383 tsmall[3][0] = 0.0;
384384 tsmall[3][1] = 0.0;
385- tsmall[3][2] = a23;
385+ tsmall[3][2] = a31;
386386 tsmall[3][3] = 1.0 - a23;
387387 }
388388
@@ -406,20 +406,20 @@ void libtetrabz_tsmall_c3(
406406 *v = a23 * a13 * a30;
407407
408408 tsmall[0][0] = 1.0;
409- tsmall[0][1] = 0.0;
409+ tsmall[0][1] = 1.0 - a30;
410410 tsmall[0][2] = 0.0;
411411 tsmall[0][3] = 0.0;
412- tsmall[1][0] = 1.0 - a30;
412+ tsmall[1][0] = 0.0;
413413 tsmall[1][1] = 0.0;
414- tsmall[1][2] = 0.0;
415- tsmall[1][3] = a30;
414+ tsmall[1][2] = a13;
415+ tsmall[1][3] = 0.0;
416416 tsmall[2][0] = 0.0;
417- tsmall[2][1] = a13;
417+ tsmall[2][1] = 0.0;
418418 tsmall[2][2] = 0.0;
419- tsmall[2][3] = 1.0 - a13;
419+ tsmall[2][3] = a23;
420420 tsmall[3][0] = 0.0;
421- tsmall[3][1] = 0.0;
422- tsmall[3][2] = a23;
421+ tsmall[3][1] = a30;
422+ tsmall[3][2] = 1.0 - a13;
423423 tsmall[3][3] = 1.0 - a23;
424424 }
425425
@@ -444,17 +444,17 @@ void libtetrabz_triangle_a1(
444444 *v = 3.0 * a10 * a20 / (e[3] - e[0]);
445445
446446 tsmall[0][0] = 1.0 - a10;
447- tsmall[0][1] = a10;
448- tsmall[0][2] = 0.0;
449- tsmall[0][3] = 0.0;
450- tsmall[1][0] = 1.0 - a20;
447+ tsmall[0][1] = 1.0 - a20;
448+ tsmall[0][2] = 1.0 - a30;
449+ tsmall[1][0] = a10;
451450 tsmall[1][1] = 0.0;
452- tsmall[1][2] = a20;
453- tsmall[1][3] = 0.0;
454- tsmall[2][0] = 1.0 - a30;
455- tsmall[2][1] = 0.0;
451+ tsmall[1][2] = 0.0;
452+ tsmall[2][0] = 0.0;
453+ tsmall[2][1] = a20;
456454 tsmall[2][2] = 0.0;
457- tsmall[2][3] = a30;
455+ tsmall[3][0] = 0.0;
456+ tsmall[3][1] = 0.0;
457+ tsmall[3][2] = a30;
458458 }
459459
460460 void libtetrabz_triangle_b1(
@@ -476,17 +476,17 @@ void libtetrabz_triangle_b1(
476476 *v = 3.0 * a30 * a13 / (e[2] - e[0]);
477477
478478 tsmall[0][0] = 1.0 - a20;
479- tsmall[0][1] = 0.0;
480- tsmall[0][2] = a20;
481- tsmall[0][3] = 0.0;
482- tsmall[1][0] = 1.0 - a30;
479+ tsmall[0][1] = 1.0 - a30;
480+ tsmall[0][2] = 0.0;
481+ tsmall[1][0] = 0.0;
483482 tsmall[1][1] = 0.0;
484- tsmall[1][2] = 0.0;
485- tsmall[1][3] = a30;
486- tsmall[2][0] = 0.0;
487- tsmall[2][1] = a13;
483+ tsmall[1][2] = a13;
484+ tsmall[2][0] = a20;
485+ tsmall[2][1] = 0.0;
488486 tsmall[2][2] = 0.0;
489- tsmall[2][3] = 1.0 - a13;
487+ tsmall[3][0] = 0.0;
488+ tsmall[3][1] = a30;
489+ tsmall[3][2] = 1.0 - a13;
490490 }
491491
492492
@@ -513,16 +513,16 @@ void libtetrabz_triangle_b2(
513513
514514 tsmall[0][0] = 1.0 - a20;
515515 tsmall[0][1] = 0.0;
516- tsmall[0][2] = a20;
517- tsmall[0][3] = 0.0;
516+ tsmall[0][2] = 0.0;
518517 tsmall[1][0] = 0.0;
519518 tsmall[1][1] = a12;
520- tsmall[1][2] = 1.0 - a12;
521- tsmall[1][3] = 0.0;
522- tsmall[2][0] = 0.0;
523- tsmall[2][1] = 1.0 - a31;
519+ tsmall[1][2] = 1.0 - a31;
520+ tsmall[2][0] = a20;
521+ tsmall[2][1] = 1.0 - a12;
524522 tsmall[2][2] = 0.0;
525- tsmall[2][3] = a31;
523+ tsmall[3][0] = 0.0;
524+ tsmall[3][1] = 0.0;
525+ tsmall[3][2]= a31;
526526 }
527527
528528 void libtetrabz_triangle_c1(
@@ -545,15 +545,15 @@ Cut triangle C1
545545 tsmall[0][0] = a03;
546546 tsmall[0][1] = 0.0;
547547 tsmall[0][2] = 0.0;
548- tsmall[0][3] = 1.0 - a03;
549548 tsmall[1][0] = 0.0;
550549 tsmall[1][1] = a13;
551550 tsmall[1][2] = 0.0;
552- tsmall[1][3] = 1.0 - a13;
553551 tsmall[2][0] = 0.0;
554552 tsmall[2][1] = 0.0;
555553 tsmall[2][2] = a23;
556- tsmall[2][3] = 1.0 - a23;
554+ tsmall[3][0] = 1.0 - a03;
555+ tsmall[3][1] = 1.0 - a13;
556+ tsmall[3][2] = 1.0 - a23;
557557 }
558558
559559 void eig_sort(
--- a/src/c/libtetrabz_dbldelta.c
+++ b/src/c/libtetrabz_dbldelta.c
@@ -37,10 +37,7 @@ void libtetrabz_dbldelta2(
3737
3838 for (ib = 0; ib < nb; ib++) {
3939
40- for (ii = 0; ii < 3; ii++)
41- w[ib][ii] = 0.0;
42-
43- for (ii = 0; ii < 3; ii++) e[ii] = ej[ii][ib];
40+ for (ii = 0; ii < 3; ii++) e[ii] = ej[ib][ii];
4441 eig_sort(3, e, indx);
4542
4643 if (e[2] < 1.0e-10) {
@@ -69,6 +66,10 @@ void libtetrabz_dbldelta2(
6966 w[ib][indx[1]] = v * a12;
7067 w[ib][indx[2]] = v * (2.0 - a02 - a12);
7168 }
69+ else {
70+ for (ii = 0; ii < 3; ii++)
71+ w[ib][ii] = 0.0;
72+ }
7273 }
7374 }
7475
@@ -86,43 +87,43 @@ void dbldelta(
8687 Main SUBROUTINE for Delta(E1) * Delta(E2)
8788 */
8889 int it, ik, ib, ii, jj, jb, ** ikv, indx[4];
89- double wlsm[4][20], ** ei1, ** ej1, ** ej2, e[4], *** w1, **w2, v, tsmall[3][4], thr;
90+ double wlsm[20][4], ** ei1, ** ej1, ** ej2, e[4], *** w1, **w2, v, tsmall[4][3], thr;
9091
9192 thr = 1.0e-10;
9293
93- ikv = (int**)calloc(6 * nk, sizeof(int*));
94- ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
94+ ikv = (int**)malloc(6 * nk * sizeof(int*));
95+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
9596 for (ik = 0; ik < 6 * nk; ik++) {
9697 ikv[ik] = ikv[0] + ik * 20;
9798 }
9899
99- ei1 = (double**)calloc(4, sizeof(double*));
100- ej1 = (double**)calloc(4, sizeof(double*));
101- ei1[0] = (double*)calloc(4 * nb, sizeof(double));
102- ej1[0] = (double*)calloc(4 * nb, sizeof(double));
100+ ei1 = (double**)malloc(4 * sizeof(double*));
101+ ej1 = (double**)malloc(4 * sizeof(double*));
102+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
103+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
103104 for (ii = 0; ii < 4; ii++) {
104105 ei1[ii] = ei1[0] + ii * nb;
105106 ej1[ii] = ej1[0] + ii * nb;
106107 }
107108
108- w1 = (double***)malloc(nb, sizeof(double**));
109- w1[0] = (double**)malloc(nb * nb, sizeof(double*));
110- w1[0][0] = (double*)malloc(nb * nb * 4, sizeof(double));
109+ w1 = (double***)malloc(nb * sizeof(double**));
110+ w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
111+ w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
111112 for (ib = 0; ib < nb; ib++) {
112- w1[ii] = w1[0] + ib * nb;
113- for (jb = 0; jb < nb; jb++) {
114- w1[ib][jb] = w1[0][0] + ib * nb * 4 + jb * 4;
113+ w1[ib] = w1[0] + ib * 4;
114+ for (ii = 0; ii < 4; ii++) {
115+ w1[ib][ii] = w1[0][0] + ib * 4 * nb + ii * nb;
115116 }
116117 }
117118
118- ej2 = (double**)calloc(3, sizeof(double*));
119- ej2[0] = (double*)calloc(3 * nb, sizeof(double));
120- for (ii = 0; ii < 3; ii++) {
121- ej2[ii] = ej2[0] + ii * nb;
119+ ej2 = (double**)malloc(nb * sizeof(double*));
120+ ej2[0] = (double*)malloc(nb * 3 * sizeof(double));
121+ for (ib = 0; ib < nb; ib++) {
122+ ej2[ib] = ej2[0] + ib * 3;
122123 }
123124
124- w2 = (double**)calloc(nb, sizeof(double*));
125- w2[0] = (double*)calloc(nb*3, sizeof(double));
125+ w2 = (double**)malloc(nb * sizeof(double*));
126+ w2[0] = (double*)malloc(nb*3 * sizeof(double));
126127 for (ib = 0; ib < nb; ib++) {
127128 w2[ib] = w2[0] + ib * 3;
128129 }
@@ -136,24 +137,26 @@ void dbldelta(
136137
137138 for (it = 0; it < 6 * nk; it++) {
138139
139- for (ii = 0; ii < 4; ii++) {
140+ for (ii = 0; ii < 4; ii++)
140141 for (ib = 0; ib < nb; ib++) {
141142 ei1[ii][ib] = 0.0;
142143 ej1[ii][ib] = 0.0;
143- for (jj = 0; jj < 20; jj++) {
144- ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
145- ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
144+ }
145+ for (jj = 0; jj < 20; jj++) {
146+ for (ii = 0; ii < 4; ii++) {
147+ for (ib = 0; ib < nb; ib++) {
148+ ei1[ii][ib] += eig1[ikv[it][jj]][ib] * wlsm[jj][ii];
149+ ej1[ii][ib] += eig2[ikv[it][jj]][ib] * wlsm[jj][ii];
146150 }
147151 }
148152 }
149153
150- for (ib = 0; ib < nb; ib++)
151- for (jb = 0; jb < nb; jb++)
152- for (ii = 0; ii < 4; ii++)
153- w1[ib][jb][ii] = 0.0;
154-
155154 for (ib = 0; ib < nb; ib++) {
156155
156+ for (ii = 0; ii < 4; ii++)
157+ for (jb = 0; jb < nb; jb++)
158+ w1[ib][ii][jb] = 0.0;
159+
157160 for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
158161 eig_sort(4, e, indx);
159162
@@ -162,15 +165,17 @@ void dbldelta(
162165 libtetrabz_triangle_a1(e, 0.0, &v, tsmall);
163166
164167 if (v > thr) {
165- for (ii = 0; ii < 3; ii++)
166- for (jj = 0; jj < 4; jj++)
168+ for (jj = 0; jj < 3; jj++)
169+ for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
170+ for (ii = 0; ii < 4; ii++)
171+ for (jj = 0; jj < 3; jj++)
167172 for (jb = 0; jb < nb; jb++)
168- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
173+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
169174 libtetrabz_dbldelta2(nb, ej2, w2);
170- for (ii = 0; ii < 4; ii++)
175+ for (ii = 0; ii < 4; ii++)
176+ for (jb = 0; jb < nb; jb++)
171177 for (jj = 0; jj < 3; jj++)
172- for (jb = 0; jb < nb; jb++)
173- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
178+ w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
174179 }
175180 }
176181 else if (e[1] < 0.0 && 0.0 <= e[2]) {
@@ -178,28 +183,33 @@ void dbldelta(
178183 libtetrabz_triangle_b1(e, 0.0, &v, tsmall);
179184
180185 if (v > thr) {
181- for (ii = 0; ii < 3; ii++)
182- for (jj = 0; jj < 4; jj++)
183- for (jb = 0; jb < nb; jb++)
184- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
185- libtetrabz_dbldelta2(nb, ej2, w2);
186+ for (jj = 0; jj < 3; jj++)
187+ for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
186188 for (ii = 0; ii < 4; ii++)
187189 for (jj = 0; jj < 3; jj++)
188190 for (jb = 0; jb < nb; jb++)
189- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
191+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
192+ libtetrabz_dbldelta2(nb, ej2, w2);
193+ for (ii = 0; ii < 4; ii++)
194+ for (jb = 0; jb < nb; jb++)
195+ for (jj = 0; jj < 3; jj++)
196+ w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
190197 }
198+
191199 libtetrabz_triangle_b2(e, 0.0, &v, tsmall);
192200
193201 if (v > thr) {
194- for (ii = 0; ii < 3; ii++)
195- for (jj = 0; jj < 4; jj++)
196- for (jb = 0; jb < nb; jb++)
197- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
198- libtetrabz_dbldelta2(nb, ej2, w2);
202+ for (jj = 0; jj < 3; jj++)
203+ for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
199204 for (ii = 0; ii < 4; ii++)
200205 for (jj = 0; jj < 3; jj++)
201206 for (jb = 0; jb < nb; jb++)
202- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
207+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
208+ libtetrabz_dbldelta2(nb, ej2, w2);
209+ for (ii = 0; ii < 4; ii++)
210+ for (jb = 0; jb < nb; jb++)
211+ for (jj = 0; jj < 3; jj++)
212+ w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
203213 }
204214 }
205215 else if (e[2] < 0.0 && 0.0 < e[3]) {
@@ -207,15 +217,17 @@ void dbldelta(
207217 libtetrabz_triangle_c1(e, 0.0, &v, tsmall);
208218
209219 if (v > thr) {
210- for (ii = 0; ii < 3; ii++)
211- for (jj = 0; jj < 4; jj++)
212- for (jb = 0; jb < nb; jb++)
213- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
214- libtetrabz_dbldelta2(nb, ej2, w2);
220+ for (jj = 0; jj < 3; jj++)
221+ for (jb = 0; jb < nb; jb++) ej2[jb][jj] = 0.0;
215222 for (ii = 0; ii < 4; ii++)
216223 for (jj = 0; jj < 3; jj++)
217224 for (jb = 0; jb < nb; jb++)
218- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
225+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
226+ libtetrabz_dbldelta2(nb, ej2, w2);
227+ for (ii = 0; ii < 4; ii++)
228+ for (jb = 0; jb < nb; jb++)
229+ for (jj = 0; jj < 3; jj++)
230+ w1[ib][indx[ii]][jb] += v * tsmall[ii][jj] * w2[jb][jj];
219231 }
220232 }
221233 else {
@@ -223,10 +235,10 @@ void dbldelta(
223235 }
224236 }
225237 for (ii = 0; ii < 20; ii++)
226- for (jj = 0; jj < 4; jj++)
227- for (ib = 0; ib < nb; ib++)
238+ for (ib = 0; ib < nb; ib++)
239+ for (jj = 0; jj < 4; jj++)
228240 for (jb = 0; jb < nb; jb++)
229- wght[ikv[it][ii]][ib][jb] += w1[ib][jb][jj] * wlsm[jj][ii];
241+ wght[ikv[it][ii]][ib][jb] += wlsm[ii][jj] * w1[ib][jj][jb];
230242 }
231243 for (ik = 0; ik < nk; ik++)
232244 for (ib = 0; ib < nb; ib++)
--- a/src/c/libtetrabz_dblstep.c
+++ b/src/c/libtetrabz_dblstep.c
@@ -22,6 +22,7 @@
2222 */
2323 #include "libtetrabz_common.h"
2424 #include <math.h>
25+#include <stdlib.h>
2526
2627 void libtetrabz_dblstep2(
2728 int nb,
@@ -42,7 +43,7 @@ void libtetrabz_dblstep2(
4243 for (ii = 0; ii < 3; ii++)
4344 w1[ib][ii] = 0.0;
4445
45- for (ii = 0; ii < 3; ii++) e[ii] = - ei[ii] + ej[ii][ib];
46+ for (ii = 0; ii < 3; ii++) e[ii] = - ei[ii] + ej[ib][ii];
4647 eig_sort(4, e, indx);
4748
4849 if (fabs(e[0]) < thr && fabs(e[3]) < thr) {
@@ -57,41 +58,41 @@ void libtetrabz_dblstep2(
5758 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
5859 for (ii = 0; ii < 4; ii++)
5960 for (jj = 0; jj < 4; jj++)
60- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
61+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
6162 }
6263 else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
6364
6465 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
6566 for (ii = 0; ii < 4; ii++)
6667 for (jj = 0; jj < 4; jj++)
67- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
68+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
6869
6970 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
7071 for (ii = 0; ii < 4; ii++)
7172 for (jj = 0; jj < 4; jj++)
72- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
73+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
7374
7475 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
7576 for (ii = 0; ii < 4; ii++)
7677 for (jj = 0; jj < 4; jj++)
77- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
78+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
7879 }
7980 else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
8081
8182 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
8283 for (ii = 0; ii < 4; ii++)
8384 for (jj = 0; jj < 4; jj++)
84- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
85+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
8586
8687 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
8788 for (ii = 0; ii < 4; ii++)
8889 for (jj = 0; jj < 4; jj++)
89- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
90+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
9091
9192 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
9293 for (ii = 0; ii < 4; ii++)
9394 for (jj = 0; jj < 4; jj++)
94- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
95+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
9596 }
9697 else if (e[3] <= 0.0) {
9798 for (ii = 0; ii < 4; ii++)
@@ -114,23 +115,45 @@ void dblstep(
114115 Main SUBROUTINE for Theta(- E1) * Theta(E1 - E2)
115116 */
116117 int it, ik, ib, ii, jj, jb, **ikv, indx[4];
117- double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], ***w1, **w2, v, tsmall[4][4], thr;
118+ double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], *** w1, ** w2, v, tsmall[4][4], thr;
118119
119- ikv = (int**)calloc(6 * nk, sizeof(int*));
120- ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
120+ ikv = (int**)malloc(6 * nk * sizeof(int*));
121+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
121122 for (ik = 0; ik < 6 * nk; ik++) {
122123 ikv[ik] = ikv[0] + ik * 20;
123124 }
124125
125- ei1 = (double**)calloc(4, sizeof(double*));
126- ej1 = (double**)calloc(4, sizeof(double*));
127- ei1[0] = (double*)calloc(4 * nb, sizeof(double));
128- ej1[0] = (double*)calloc(4 * nb, sizeof(double));
126+ ei1 = (double**)malloc(4 * sizeof(double*));
127+ ej1 = (double**)malloc(4 * sizeof(double*));
128+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
129+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
129130 for (ii = 0; ii < 4; ii++) {
130131 ei1[ii] = ei1[0] + ii * nb;
131132 ej1[ii] = ej1[0] + ii * nb;
132133 }
133134
135+ w1 = (double***)malloc(nb * sizeof(double**));
136+ w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
137+ w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
138+ for (ib = 0; ib < nb; ib++) {
139+ 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;
142+ }
143+ }
144+
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+
151+ w2 = (double**)malloc(nb * sizeof(double*));
152+ w2[0] = (double*)malloc(nb * 4 * sizeof(double));
153+ for (ib = 0; ib < nb; ib++) {
154+ w2[ib] = w2[0] + ib * 4;
155+ }
156+
134157 libtetrabz_initialize(ng, bvec, wlsm, ikv);
135158
136159 for (ik = 0; ik < nk; ik++)
@@ -142,23 +165,25 @@ void dblstep(
142165
143166 for(it = 0; it < 6*nk; it++) {
144167
145- for (ii = 0; ii < 4; ii++) {
168+ for (ii = 0; ii < 4; ii++)
146169 for (ib = 0; ib < nb; ib++) {
147170 ei1[ii][ib] = 0.0;
148171 ej1[ii][ib] = 0.0;
149- for (jj = 0; jj < 20; jj++) {
150- ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
151- ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
172+ }
173+ for (jj = 0; jj < 20; jj++) {
174+ for (ii = 0; ii < 4; ii++) {
175+ 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];
152178 }
153179 }
154180 }
155-
156181
157182 for (ib = 0; ib < nb; ib++) {
158183
159- for (jb = 0; jb < nb; jb++)
160- for (ii = 0; ii < 4; ii++)
161- w1[ib][jb][ii] = 0.0;
184+ for (ii = 0; ii < 4; ii++)
185+ for (jb = 0; jb < nb; jb++)
186+ w1[ib][ii][jb] = 0.0;
162187
163188 for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
164189 eig_sort(4, e, indx);
@@ -168,17 +193,21 @@ void dblstep(
168193 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
169194
170195 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;
199+ }
171200 for (ii = 0; ii < 4; ii++)
172201 for (jj = 0; jj < 4; jj++) {
173- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
202+ ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
174203 for (jb = 0; jb < nb; jb++)
175- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
204+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
176205 }
177206 libtetrabz_dblstep2(nb, ei2, ej2, w2);
178207 for (ii = 0; ii < 4; ii++)
179- for (jj = 0; jj < 4; jj++)
180- for (jb = 0; jb < nb; jb++)
181- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
208+ 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];
182211 }
183212 }
184213 else if (e[1] <= 0.0 && 0.0 < e[2]) {
@@ -186,49 +215,61 @@ void dblstep(
186215 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
187216
188217 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;
221+ }
189222 for (ii = 0; ii < 4; ii++)
190223 for (jj = 0; jj < 4; jj++) {
191- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
224+ ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
192225 for (jb = 0; jb < nb; jb++)
193- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
226+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
194227 }
195228 libtetrabz_dblstep2(nb, ei2, ej2, w2);
196229 for (ii = 0; ii < 4; ii++)
197- for (jj = 0; jj < 4; jj++)
198- for (jb = 0; jb < nb; jb++)
199- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
230+ 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];
200233 }
201234
202235 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
203236
204237 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;
241+ }
205242 for (ii = 0; ii < 4; ii++)
206243 for (jj = 0; jj < 4; jj++) {
207- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
244+ ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
208245 for (jb = 0; jb < nb; jb++)
209- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
246+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
210247 }
211248 libtetrabz_dblstep2(nb, ei2, ej2, w2);
212249 for (ii = 0; ii < 4; ii++)
213- for (jj = 0; jj < 4; jj++)
214- for (jb = 0; jb < nb; jb++)
215- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
250+ 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];
216253 }
217254
218255 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
219256
220257 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;
261+ }
221262 for (ii = 0; ii < 4; ii++)
222263 for (jj = 0; jj < 4; jj++) {
223- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
264+ ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
224265 for (jb = 0; jb < nb; jb++)
225- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
266+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
226267 }
227268 libtetrabz_dblstep2(nb, ei2, ej2, w2);
228269 for (ii = 0; ii < 4; ii++)
229- for (jj = 0; jj < 4; jj++)
230- for (jb = 0; jb < nb; jb++)
231- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
270+ 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];
232273 }
233274 }
234275 else if (e[2] <= 0.0 && 0.0 < e[3]) {
@@ -236,61 +277,73 @@ void dblstep(
236277 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
237278
238279 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;
283+ }
239284 for (ii = 0; ii < 4; ii++)
240285 for (jj = 0; jj < 4; jj++) {
241- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
286+ ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
242287 for (jb = 0; jb < nb; jb++)
243- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
288+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
244289 }
245290 libtetrabz_dblstep2(nb, ei2, ej2, w2);
246291 for (ii = 0; ii < 4; ii++)
247- for (jj = 0; jj < 4; jj++)
248- for (jb = 0; jb < nb; jb++)
249- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
292+ 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];
250295 }
251296
252297 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
253298
254299 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;
303+ }
255304 for (ii = 0; ii < 4; ii++)
256305 for (jj = 0; jj < 4; jj++) {
257- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
306+ ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
258307 for (jb = 0; jb < nb; jb++)
259- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
308+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
260309 }
261310 libtetrabz_dblstep2(nb, ei2, ej2, w2);
262311 for (ii = 0; ii < 4; ii++)
263- for (jj = 0; jj < 4; jj++)
264- for (jb = 0; jb < nb; jb++)
265- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
312+ 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];
266315 }
267316
268317 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
269318
270319 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;
323+ }
271324 for (ii = 0; ii < 4; ii++)
272325 for (jj = 0; jj < 4; jj++) {
273- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
326+ ei2[jj] += ej1[indx[ii]][ib] * tsmall[ii][jj];
274327 for (jb = 0; jb < nb; jb++)
275- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
328+ ej2[jb][jj] += ej1[indx[ii]][jb] * tsmall[ii][jj];
276329 }
277330 libtetrabz_dblstep2(nb, ei2, ej2, w2);
278331 for (ii = 0; ii < 4; ii++)
279- for (jj = 0; jj < 4; jj++)
280- for (jb = 0; jb < nb; jb++)
281- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
332+ 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];
282335 }
283336 }
284337 else if (e[3] <= 0.0) {
285338 for (ii = 0; ii < 4; ii++) {
286339 ei2[ii] = ej1[ii][ib];
287340 for (jb = 0; jb < nb; jb++)
288- ej2[ii][jb] = ej1[ii][jb];
341+ ej2[jb][ii] = ej1[ii][jb];
289342 }
290343 libtetrabz_dblstep2(nb, ei2, ej2, w2);
291344 for (ii = 0; ii < 4; ii++)
292345 for (jb = 0; jb < nb; jb++)
293- w1[ib][jb][ii] += v * w2[jb][ii];
346+ w1[ib][ii][jb] += w2[jb][ii];
294347 }
295348 else {
296349 continue;
@@ -306,4 +359,18 @@ void dblstep(
306359 for (ib = 0; ib < nb; ib++)
307360 for (jb = 0; jb < nb; jb++)
308361 wght[ik][ib][jb] /= (6.0 * (double) nk);
362+
363+ free(ikv[0]);
364+ free(ikv);
365+ free(ei1[0]);
366+ free(ei1);
367+ free(ej1[0]);
368+ free(ej1);
369+ free(ej2[0]);
370+ free(ej2);
371+ free(w1[0][0]);
372+ free(w1[0]);
373+ free(w1);
374+ free(w2[0]);
375+ free(w2);
309376 }
--- a/src/c/libtetrabz_dos.c
+++ b/src/c/libtetrabz_dos.c
@@ -36,23 +36,23 @@ void dos(
3636 Compute Dos : Delta(E - E1)
3737 */
3838 int it, ik, ib, ii, jj, ie, ** ikv, indx[4];
39- double wlsm[4][20], ** ei1, e[4], *** w1, v, tsmall[3][4];
39+ double wlsm[20][4], ** ei1, e[4], *** w1, v, tsmall[4][3];
4040
41- ikv = (int**)calloc(6 * nk, sizeof(int*));
42- ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
41+ ikv = (int**)malloc(6 * nk * sizeof(int*));
42+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
4343 for (ik = 0; ik < 6 * nk; ik++) {
4444 ikv[ik] = ikv[0] + ik * 20;
4545 }
4646
47- ei1 = (double**)calloc(4, sizeof(double*));
48- ei1[0] = (double*)calloc(4 * nb, sizeof(double));
47+ ei1 = (double**)malloc(4 * sizeof(double*));
48+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
4949 for (ii = 0; ii < 4; ii++) {
5050 ei1[ii] = ei1[0] + ii * nb;
5151 }
5252
53- w1 = (double***)calloc(nb, sizeof(double**));
54- w1[0] = (double**)calloc(nb * ne, sizeof(double*));
55- w1[0][0] = (double*)calloc(nb * ne * 4, sizeof(double));
53+ w1 = (double***)malloc(nb * sizeof(double**));
54+ w1[0] = (double**)malloc(nb * ne * sizeof(double*));
55+ w1[0][0] = (double*)malloc(nb * ne * 4 * sizeof(double));
5656 for (ib = 0; ib < nb; ib++) {
5757 w1[ii] = w1[0] + ib * ne;
5858 for (ie = 0; ie < ne; ie++) {
@@ -69,13 +69,13 @@ void dos(
6969
7070 for (it = 0; it < 6 * nk; it++) {
7171
72- for (ii = 0; ii < 4; ii++) {
73- for (ib = 0; ib < nb; ib++) {
72+ for (ii = 0; ii < 4; ii++)
73+ for (ib = 0; ib < nb; ib++)
7474 ei1[ii][ib] = 0.0;
75- for (jj = 0; jj < 20; jj++)
76- ei1[ii][ib] += wlsm[ii][jj] * eig[ikv[it][jj]][ib];
77- }
78- }
75+ for (jj = 0; jj < 20; jj++)
76+ for (ii = 0; ii < 4; ii++)
77+ for (ib = 0; ib < nb; ib++)
78+ ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
7979
8080 for (ib = 0; ib < nb; ib++)
8181 for (ie = 0; ie < ne; ie++)
@@ -94,7 +94,7 @@ void dos(
9494 libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
9595 for (ii = 0; ii < 4; ii++)
9696 for (jj = 0; jj < 3; jj++)
97- w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
97+ w1[ib][ie][ii] += v * tsmall[ii][jj] / 3.0;
9898
9999 }
100100 else if ((e[1] < e0[ie] && e0[ie] <= e[2]) || (e[1] <= e0[ie] && e0[ie] < e[2])) {
@@ -102,29 +102,28 @@ void dos(
102102 libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
103103 for (ii = 0; ii < 4; ii++)
104104 for (jj = 0; jj < 3; jj++)
105- w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
105+ w1[ib][ie][ii] += v * tsmall[ii][jj] / 3.0;
106106
107107 libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
108108 for (ii = 0; ii < 4; ii++)
109109 for (jj = 0; jj < 3; jj++)
110- w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
110+ w1[ib][ie][ii] += v * tsmall[ii][jj] / 3.0;
111111 }
112112 else if ((e[2] < e0[ie] && e0[ie] <= e[3]) || (e[2] <= e0[ie] && e0[ie] < e[3])) {
113113
114114 libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
115115 for (ii = 0; ii < 4; ii++)
116116 for (jj = 0; jj < 3; jj++)
117- w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
117+ w1[ib][ie][ii] += v * tsmall[ii][jj] / 3.0;
118118 }
119119 else {
120120 continue;
121121 }
122122 for (ii = 0; ii < 20; ii++)
123- for (jj = 0; jj < 4; jj++)
124- for (ib = 0; ib < nb; ib++)
125- for (ie = 0; ie < ne; ie++)
126- wght[ikv[it][ii]][ib][ie] += w1[ib][ie][jj] * wlsm[jj][ii];
127-
123+ for (ib = 0; ib < nb; ib++)
124+ for (ie = 0; ie < ne; ie++)
125+ for (jj = 0; jj < 4; jj++)
126+ wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
128127 }
129128 }
130129 }
@@ -157,23 +156,23 @@ void intdos(
157156 */
158157
159158 int it, ik, ib, ii, jj, ie, ** ikv, indx[4];
160- double wlsm[4][20], ** ei1, e[4], *** w1, v, tsmall[4][4];
159+ double wlsm[20][4], ** ei1, e[4], *** w1, v, tsmall[4][4];
161160
162- ikv = (int**)calloc(6 * nk, sizeof(int*));
163- ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
161+ ikv = (int**)malloc(6 * nk * sizeof(int*));
162+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
164163 for (ik = 0; ik < 6 * nk; ik++) {
165164 ikv[ik] = ikv[0] + ik * 20;
166165 }
167166
168- ei1 = (double**)calloc(4, sizeof(double*));
169- ei1[0] = (double*)calloc(4 * nb, sizeof(double));
167+ ei1 = (double**)malloc(4 * sizeof(double*));
168+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
170169 for (ii = 0; ii < 4; ii++) {
171170 ei1[ii] = ei1[0] + ii * nb;
172171 }
173172
174- w1 = (double***)calloc(nb, sizeof(double**));
175- w1[0] = (double**)calloc(nb * ne, sizeof(double*));
176- w1[0][0] = (double*)calloc(nb * ne * 4, sizeof(double));
173+ w1 = (double***)malloc(nb * sizeof(double**));
174+ w1[0] = (double**)malloc(nb * ne * sizeof(double*));
175+ w1[0][0] = (double*)malloc(nb * ne * 4 * sizeof(double));
177176 for (ib = 0; ib < nb; ib++) {
178177 w1[ii] = w1[0] + ib * ne;
179178 for (ie = 0; ie < ne; ie++) {
@@ -190,13 +189,13 @@ void intdos(
190189
191190 for (it = 0; it < 6 * nk; it++) {
192191
193- for (ii = 0; ii < 4; ii++) {
194- for (ib = 0; ib < nb; ib++) {
192+ for (ii = 0; ii < 4; ii++)
193+ for (ib = 0; ib < nb; ib++)
195194 ei1[ii][ib] = 0.0;
196- for (jj = 0; jj < 20; jj++)
197- ei1[ii][ib] += wlsm[ii][jj] * eig[ikv[it][jj]][ib];
198- }
199- }
195+ for (jj = 0; jj < 20; jj++)
196+ for (ii = 0; ii < 4; ii++)
197+ for (ib = 0; ib < nb; ib++)
198+ ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
200199
201200 for (ib = 0; ib < nb; ib++)
202201 for (ie = 0; ie < ne; ie++)
@@ -214,24 +213,24 @@ void intdos(
214213 libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
215214 for (ii = 0; ii < 4; ii++)
216215 for (jj = 0; jj < 4; jj++)
217- w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
216+ w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
218217 }
219218 else if ((e[1] <= e0[ie] && e0[ie] < e[2]) || (e[1] < e0[ie] && e0[ie] <= e[2])) {
220219
221220 libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
222221 for (ii = 0; ii < 4; ii++)
223222 for (jj = 0; jj < 4; jj++)
224- w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
223+ w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
225224
226225 libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
227226 for (ii = 0; ii < 4; ii++)
228227 for (jj = 0; jj < 4; jj++)
229- w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
228+ w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
230229
231230 libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
232231 for (ii = 0; ii < 4; ii++)
233232 for (jj = 0; jj < 4; jj++)
234- w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
233+ w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
235234
236235 }
237236 else if ((e[2] <= e0[ie] && e0[ie] < e[3]) || (e[2] < e0[ie] && e0[ie] <= e[3])) {
@@ -239,17 +238,17 @@ void intdos(
239238 libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
240239 for (ii = 0; ii < 4; ii++)
241240 for (jj = 0; jj < 4; jj++)
242- w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
241+ w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
243242
244243 libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
245244 for (ii = 0; ii < 4; ii++)
246245 for (jj = 0; jj < 4; jj++)
247- w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
246+ w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
248247
249248 libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
250249 for (ii = 0; ii < 4; ii++)
251250 for (jj = 0; jj < 4; jj++)
252- w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
251+ w1[ib][ie][ii] += v * tsmall[ii][jj] * 0.25;
253252
254253 }
255254 else if (e[3] <= e0[ie]) {
@@ -263,10 +262,10 @@ void intdos(
263262 }
264263
265264 for (ii = 0; ii < 20; ii++)
266- for (jj = 0; jj < 4; jj++)
267- for (ib = 0; ib < nb; ib++)
268- for (ie = 0; ie < ne; ie++)
269- wght[ikv[it][ii]][ib][ie] += w1[ib][ie][jj] * wlsm[jj][ii];
265+ for (ib = 0; ib < nb; ib++)
266+ for (ie = 0; ie < ne; ie++)
267+ for (jj = 0; jj < 4; jj++)
268+ wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
270269 }
271270 for (ik = 0; ik < nk; ik++)
272271 for (ib = 0; ib < nb; ib++)
--- a/src/c/libtetrabz_fermigr.c
+++ b/src/c/libtetrabz_fermigr.c
@@ -21,6 +21,7 @@
2121 // SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
2222 //
2323 #include "libtetrabz_common.h"
24+#include <stdlib.h>
2425
2526 void libtetrabz_fermigr3(
2627 int ne,
@@ -28,41 +29,41 @@ void libtetrabz_fermigr3(
2829 double de[4],
2930 double **w1
3031 ) {
31- int ii, jj, ie, indx[3];
32- double e[4], tsmall[3][4], v;
32+ int i4, j3, ie, indx[3];
33+ double e[4], tsmall[4][3], v;
3334
34- for (ii = 0; ii < 3; ii++) e[ii] = de[ii];
35- eig_sort(3, e, indx);
35+ for (i4 = 0; i4 < 4; i4++) e[i4] = de[i4];
36+ eig_sort(4, e, indx);
3637
3738 for (ie = 0; ie < ne; ie++) {
3839
39- for (ii = 0; ii < 4; ii++) w1[ie][ii] = 0.0;
40+ for (i4 = 0; i4 < 4; i4++) w1[i4][ie] = 0.0;
4041
4142 if (e[0] < e0[ie] && e0[ie] <= e[1]) {
4243
4344 libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
44- for (ii = 0; ii < 4; ii++)
45- for (jj = 0; jj < 3; jj++)
46- w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
45+ for (i4 = 0; i4 < 4; i4++)
46+ for (j3 = 0; j3 < 3; j3++)
47+ w1[i4][ie] += v * tsmall[i4][j3] / 3.0;
4748 }
4849 else if (e[1] < e0[ie] && e0[ie] <= e[2]) {
4950
5051 libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
51- for (ii = 0; ii < 4; ii++)
52- for (jj = 0; jj < 3; jj++)
53- w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
52+ for (i4 = 0; i4 < 4; i4++)
53+ for (j3 = 0; j3 < 3; j3++)
54+ w1[i4][ie] += v * tsmall[i4][j3] / 3.0;
5455
5556 libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
56- for (ii = 0; ii < 4; ii++)
57- for (jj = 0; jj < 3; jj++)
58- w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
57+ for (i4 = 0; i4 < 4; i4++)
58+ for (j3 = 0; j3 < 3; j3++)
59+ w1[i4][ie] += v * tsmall[i4][j3] / 3.0;
5960 }
6061 else if (e[2] < e0[ie] && e0[ie] < e[3]) {
6162
6263 libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
63- for (ii = 0; ii < 4; ii++)
64- for (jj = 0; jj < 3; jj++)
65- w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
64+ for (i4 = 0; i4 < 4; i4++)
65+ for (j3 = 0; j3 < 3; j3++)
66+ w1[i4][ie] += v * tsmall[i4][j3] / 3.0;
6667 }
6768 }
6869 }
@@ -75,16 +76,24 @@ void libtetrabz_fermigr2(
7576 double **ej1,
7677 double ***w1
7778 ) {
78- int ib, ii, jj, ie, indx[4];
79+ int ib, i4, j4, ie, indx[4];
7980 double e[4], tsmall[4][4], v, de[4], thr, **w2;
8081
82+ w2 = (double**)malloc(4 * sizeof(double*));
83+ w2[0] = (double*)malloc(4 * ne * sizeof(double));
84+ for (i4 = 0; i4 < 4; i4++) {
85+ w2[i4] = w2[0] + i4 * ne;
86+ }
87+
88+ thr = 1.0e-8;
89+
8190 for (ib = 0; ib < nb; ib++) {
8291
83- for (ie = 0; ie < ne; ie++)
84- for (ii = 0; ii < 4; ii++)
85- w1[ib][ie][ii] = 0.0;
92+ for (i4 = 0; i4 < 4; i4++)
93+ for (ie = 0; ie < ne; ie++)
94+ w1[ib][i4][ie] = 0.0;
8695
87- for (ii = 0; ii < 4; ii++) e[ii] = -ej1[ii][ib];
96+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[i4][ib];
8897 eig_sort(4, e, indx);
8998
9099 if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
@@ -92,16 +101,15 @@ void libtetrabz_fermigr2(
92101 libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
93102
94103 if (v > thr) {
95- for (ii = 0; ii < 4; ii++) {
96- de[ii] = 0.0;
97- for (jj = 0; jj < 4; jj++)
98- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
99- }
104+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
105+ for (i4 = 0; i4 < 4; i4++)
106+ for (j4 = 0; j4 < 4; j4++)
107+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
100108 libtetrabz_fermigr3(ne, e0, de, w2);
101- for (ii = 0; ii < 4; ii++)
102- for (jj = 0; jj < 4; jj++)
109+ for (i4 = 0; i4 < 4; i4++)
110+ for (j4 = 0; j4 < 4; j4++)
103111 for (ie = 0; ie < ne; ie++)
104- w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
112+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
105113 }
106114 }
107115 else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
@@ -109,46 +117,43 @@ void libtetrabz_fermigr2(
109117 libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
110118
111119 if (v > thr) {
112- for (ii = 0; ii < 4; ii++) {
113- de[ii] = 0.0;
114- for (jj = 0; jj < 4; jj++)
115- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
116- }
120+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
121+ for (i4 = 0; i4 < 4; i4++)
122+ for (j4 = 0; j4 < 4; j4++)
123+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
117124 libtetrabz_fermigr3(ne, e0, de, w2);
118- for (ii = 0; ii < 4; ii++)
119- for (jj = 0; jj < 4; jj++)
125+ for (i4 = 0; i4 < 4; i4++)
126+ for (j4 = 0; j4 < 4; j4++)
120127 for (ie = 0; ie < ne; ie++)
121- w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
128+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
122129 }
123130
124131 libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
125132
126133 if (v > thr) {
127- for (ii = 0; ii < 4; ii++) {
128- de[ii] = 0.0;
129- for (jj = 0; jj < 4; jj++)
130- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
131- }
134+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
135+ for (i4 = 0; i4 < 4; i4++)
136+ for (j4 = 0; j4 < 4; j4++)
137+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
132138 libtetrabz_fermigr3(ne, e0, de, w2);
133- for (ii = 0; ii < 4; ii++)
134- for (jj = 0; jj < 4; jj++)
139+ for (i4 = 0; i4 < 4; i4++)
140+ for (j4 = 0; j4 < 4; j4++)
135141 for (ie = 0; ie < ne; ie++)
136- w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
142+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
137143 }
138144
139145 libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
140146
141147 if (v > thr) {
142- for (ii = 0; ii < 4; ii++) {
143- de[ii] = 0.0;
144- for (jj = 0; jj < 4; jj++)
145- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
146- }
148+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
149+ for (i4 = 0; i4 < 4; i4++)
150+ for (j4 = 0; j4 < 4; j4++)
151+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
147152 libtetrabz_fermigr3(ne, e0, de, w2);
148- for (ii = 0; ii < 4; ii++)
149- for (jj = 0; jj < 4; jj++)
153+ for (i4 = 0; i4 < 4; i4++)
154+ for (j4 = 0; j4 < 4; j4++)
150155 for (ie = 0; ie < ne; ie++)
151- w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
156+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
152157 }
153158 }
154159 else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
@@ -156,57 +161,58 @@ void libtetrabz_fermigr2(
156161 libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
157162
158163 if (v > thr) {
159- for (ii = 0; ii < 4; ii++) {
160- de[ii] = 0.0;
161- for (jj = 0; jj < 4; jj++)
162- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
163- }
164+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
165+ for (i4 = 0; i4 < 4; i4++)
166+ for (j4 = 0; j4 < 4; j4++)
167+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
164168 libtetrabz_fermigr3(ne, e0, de, w2);
165- for (ii = 0; ii < 4; ii++)
166- for (jj = 0; jj < 4; jj++)
169+ for (i4 = 0; i4 < 4; i4++)
170+ for (j4 = 0; j4 < 4; j4++)
167171 for (ie = 0; ie < ne; ie++)
168- w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
172+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
169173 }
170174
171175 libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
172176
173177 if (v > thr) {
174- for (ii = 0; ii < 4; ii++) {
175- de[ii] = 0.0;
176- for (jj = 0; jj < 4; jj++)
177- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
178- }
178+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
179+ for (i4 = 0; i4 < 4; i4++)
180+ for (j4 = 0; j4 < 4; j4++)
181+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
179182 libtetrabz_fermigr3(ne, e0, de, w2);
180- for (ii = 0; ii < 4; ii++)
181- for (jj = 0; jj < 4; jj++)
183+ for (i4 = 0; i4 < 4; i4++)
184+ for (j4 = 0; j4 < 4; j4++)
182185 for (ie = 0; ie < ne; ie++)
183- w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
186+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
184187 }
185188
186189 libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
187190
188191 if (v > thr) {
189- for (ii = 0; ii < 4; ii++) {
190- de[ii] = 0.0;
191- for (jj = 0; jj < 4; jj++)
192- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
193- }
192+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
193+ for (i4 = 0; i4 < 4; i4++)
194+ for (j4 = 0; j4 < 4; j4++)
195+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
194196 libtetrabz_fermigr3(ne, e0, de, w2);
195- for (ii = 0; ii < 4; ii++)
196- for (jj = 0; jj < 4; jj++)
197+ for (i4 = 0; i4 < 4; i4++)
198+ for (j4 = 0; j4 < 4; j4++)
197199 for (ie = 0; ie < ne; ie++)
198- w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
200+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
199201 }
200202 }
201203 else if (e[3] <= 0.0) {
202- for (ii = 0; ii < 4; ii++)
203- de[ii] = ej1[ii][ib] - ei1[ii];
204+ for (i4 = 0; i4 < 4; i4++) {
205+ de[i4] = tsmall[i4][j4] * (ej1[i4][ib] - ei1[i4]);
206+ }
204207 libtetrabz_fermigr3(ne, e0, de, w2);
205- for (ii = 0; ii < 4; ii++)
208+ for (i4 = 0; i4 < 4; i4++)
206209 for (ie = 0; ie < ne; ie++)
207- w1[ib][ie][ii] += w2[ie][ii];
210+ w1[ib][i4][ie] += w2[i4][ie];
208211 }
209212 }
213+
214+ free(w2[0]);
215+ free(w2);
210216 }
211217
212218 void fermigr(
@@ -224,8 +230,52 @@ void fermigr(
224230 /*
225231 Main SUBROUTINE for Fermi's Golden rule : Theta(- E1) * Theta(E2) * Delta(E2 - E1 - w)
226232 */
227- int it, ik, ib, ii, jj, jb, ie, **ikv, indx[4];
228- double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], ****w1, ***w2, v, tsmall[4][4], thr;
233+ int it, ik, ib, i20, i4, j4, jb, ie, **ikv, indx[4];
234+ double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], **** w1, *** w2, v, tsmall[4][4], thr;
235+
236+ ikv = (int**)malloc(6 * nk * sizeof(int*));
237+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
238+ for (it = 0; it < 6 * nk; it++) {
239+ ikv[it] = ikv[0] + it * 20;
240+ }
241+
242+ ei1 = (double**)malloc(4 * sizeof(double*));
243+ ej1 = (double**)malloc(4 * sizeof(double*));
244+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
245+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
246+ for (i4 = 0; i4 < 4; i4++) {
247+ ei1[i4] = ei1[0] + i4 * nb;
248+ ej1[i4] = ej1[0] + i4 * nb;
249+ }
250+
251+ ej2 = (double**)malloc(nb * sizeof(double*));
252+ ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
253+ for (ib = 0; ib < nb; ib++)
254+ ej2[ib] = ej2[0] + ib * 4;
255+
256+ w1 = (double****)malloc(nb * sizeof(double***));
257+ w1[0] = (double***)malloc(nb * 4 * sizeof(double**));
258+ w1[0][0] = (double**)malloc(nb * 4 * nb * sizeof(double*));
259+ w1[0][0][0] = (double*)malloc(nb * 4 * nb * ne * sizeof(double));
260+ for (ib = 0; ib < nb; ib++) {
261+ w1[ib] = w1[0] + ib * 4;
262+ for (i4 = 0; i4 < 4; i4++) {
263+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
264+ for (jb = 0; jb < nb; jb++) {
265+ w1[ib][i4][jb] = w1[0][0][0] + ib * 4 * nb * ne + i4 * nb * ne + jb * ne;
266+ }
267+ }
268+ }
269+
270+ w2 = (double***)malloc(nb * sizeof(double**));
271+ w2[0] = (double**)malloc(nb * 4 * sizeof(double*));
272+ w2[0][0] = (double*)malloc(nb * 4 * ne * sizeof(double));
273+ for (ib = 0; ib < nb; ib++) {
274+ w2[ib] = w2[0] + ib * 4;
275+ for (i4 = 0; i4 < 4; i4++) {
276+ w2[ib][i4] = w2[0][0] + ib * 4 * ne + i4 * ne;
277+ }
278+ }
229279
230280 libtetrabz_initialize(ng, bvec, wlsm, ikv);
231281
@@ -239,17 +289,16 @@ void fermigr(
239289
240290 for(it = 0; it < 6*nk; it++) {
241291
242- for (ii = 0; ii < 4; ii++) {
292+ for (i4 = 0; i4 < 4; i4++)
243293 for (ib = 0; ib < nb; ib++) {
244- ei1[ii][ib] = 0.0;
245- ej1[ii][ib] = 0.0;
294+ ei1[i4][ib] = 0.0;
295+ ej1[i4][ib] = 0.0;
246296 }
247- }
248- for (jj = 0; jj < 20; jj++) {
249- for (ib = 0; ib < nb; ib++) {
250- for (ii = 0; ii < 4; ii++) {
251- ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
252- ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
297+ for (i20 = 0; i20 < 20; i20++) {
298+ for (i4 = 0; i4 < 4; i4++) {
299+ for (ib = 0; ib < nb; ib++) {
300+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
301+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
253302 }
254303 }
255304 }
@@ -257,11 +306,11 @@ void fermigr(
257306 for (ib = 0; ib < nb; ib++) {
258307
259308 for (jb = 0; jb < nb; jb++)
260- for (ii = 0; ii < 4; ii++)
309+ for (i4 = 0; i4 < 4; i4++)
261310 for (ie = 0; ie < ne; ie++)
262- w1[ib][jb][ie][ii] = 0.0;
311+ w1[ib][jb][ie][i4] = 0.0;
263312
264- for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
313+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
265314 eig_sort(4, e, indx);
266315
267316 if (e[0] <= 0.0 && 0.0 < e[1]) {
@@ -269,21 +318,22 @@ void fermigr(
269318 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
270319
271320 if (v > thr) {
272- for (ii = 0; ii < 4; ii++) {
273- ei2[ii] = 0.0;
274- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
275- for (jj = 0; jj < 4; jj++) {
276- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
321+ for (j4 = 0; j4 < 4; j4++) {
322+ ei2[j4] = 0.0;
323+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
324+ }
325+ for (i4 = 0; i4 < 4; i4++)
326+ for (j4 = 0; j4 < 4; j4++) {
327+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
277328 for (jb = 0; jb < nb; jb++)
278- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
329+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
279330 }
280- }
281331 libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
282- for (ii = 0; ii < 4; ii++)
283- for (jj = 0; jj < 4; jj++)
284- for (jb = 0; jb < nb; jb++)
332+ for (i4 = 0; i4 < 4; i4++)
333+ for (jb = 0; jb < nb; jb++)
334+ for (j4 = 0; j4 < 4; j4++)
285335 for (ie = 0; ie < ne; ie++)
286- w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
336+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
287337 }
288338 }
289339 else if (e[1] <= 0.0 && 0.0 < e[2]) {
@@ -291,61 +341,64 @@ void fermigr(
291341 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
292342
293343 if (v > thr) {
294- for (ii = 0; ii < 4; ii++) {
295- ei2[ii] = 0.0;
296- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
297- for (jj = 0; jj < 4; jj++) {
298- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
344+ for (j4 = 0; j4 < 4; j4++) {
345+ ei2[j4] = 0.0;
346+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
347+ }
348+ for (i4 = 0; i4 < 4; i4++)
349+ for (j4 = 0; j4 < 4; j4++) {
350+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
299351 for (jb = 0; jb < nb; jb++)
300- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
352+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
301353 }
302- }
303354 libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
304- for (ii = 0; ii < 4; ii++)
305- for (jj = 0; jj < 4; jj++)
306- for (jb = 0; jb < nb; jb++)
355+ for (i4 = 0; i4 < 4; i4++)
356+ for (jb = 0; jb < nb; jb++)
357+ for (j4 = 0; j4 < 4; j4++)
307358 for (ie = 0; ie < ne; ie++)
308- w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
359+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
309360 }
310361
311362 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
312363
313364 if (v > thr) {
314- for (ii = 0; ii < 4; ii++) {
315- ei2[ii] = 0.0;
316- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
317- for (jj = 0; jj < 4; jj++) {
318- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
365+ for (j4 = 0; j4 < 4; j4++) {
366+ ei2[j4] = 0.0;
367+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
368+ }
369+ for (i4 = 0; i4 < 4; i4++)
370+ for (j4 = 0; j4 < 4; j4++) {
371+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
319372 for (jb = 0; jb < nb; jb++)
320- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
373+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
321374 }
322- }
323375 libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
324- for (ii = 0; ii < 4; ii++)
325- for (jj = 0; jj < 4; jj++)
326- for (jb = 0; jb < nb; jb++)
376+ for (i4 = 0; i4 < 4; i4++)
377+ for (jb = 0; jb < nb; jb++)
378+ for (j4 = 0; j4 < 4; j4++)
327379 for (ie = 0; ie < ne; ie++)
328- w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
380+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
329381 }
330382
331383 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
332384
333385 if (v > thr) {
334- for (ii = 0; ii < 4; ii++) {
335- ei2[ii] = 0.0;
336- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
337- for (jj = 0; jj < 4; jj++) {
338- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
386+ for (j4 = 0; j4 < 4; j4++) {
387+ ei2[j4] = 0.0;
388+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
389+ }
390+ for (i4 = 0; i4 < 4; i4++)
391+ for (j4 = 0; j4 < 4; j4++) {
392+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
339393 for (jb = 0; jb < nb; jb++)
340- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
394+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
341395 }
342- }
343396 libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
344- for (ii = 0; ii < 4; ii++)
345- for (jj = 0; jj < 4; jj++)
346- for (jb = 0; jb < nb; jb++)
397+ for (i4 = 0; i4 < 4; i4++)
398+ for (jb = 0; jb < nb; jb++)
399+ for (j4 = 0; j4 < 4; j4++)
347400 for (ie = 0; ie < ne; ie++)
348- w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
401+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
349402 }
350403 }
351404 else if (e[2] <= 0.0 && 0.0 < e[3]) {
@@ -353,87 +406,89 @@ void fermigr(
353406 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
354407
355408 if (v > thr) {
356- for (ii = 0; ii < 4; ii++) {
357- ei2[ii] = 0.0;
358- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
359- for (jj = 0; jj < 4; jj++) {
360- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
409+ for (j4 = 0; j4 < 4; j4++) {
410+ ei2[j4] = 0.0;
411+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
412+ }
413+ for (i4 = 0; i4 < 4; i4++)
414+ for (j4 = 0; j4 < 4; j4++) {
415+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
361416 for (jb = 0; jb < nb; jb++)
362- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
417+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
363418 }
364- }
365419 libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
366- for (ii = 0; ii < 4; ii++)
367- for (jj = 0; jj < 4; jj++)
368- for (jb = 0; jb < nb; jb++)
420+ for (i4 = 0; i4 < 4; i4++)
421+ for (jb = 0; jb < nb; jb++)
422+ for (j4 = 0; j4 < 4; j4++)
369423 for (ie = 0; ie < ne; ie++)
370- w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
424+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
371425 }
372426
373427 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
374428
375429 if (v > thr) {
376- for (ii = 0; ii < 4; ii++) {
377- ei2[ii] = 0.0;
378- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
379- for (jj = 0; jj < 4; jj++) {
380- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
430+ for (j4 = 0; j4 < 4; j4++) {
431+ ei2[j4] = 0.0;
432+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
433+ }
434+ for (i4 = 0; i4 < 4; i4++)
435+ for (j4 = 0; j4 < 4; j4++) {
436+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
381437 for (jb = 0; jb < nb; jb++)
382- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
438+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
383439 }
384- }
385440 libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
386- for (ii = 0; ii < 4; ii++)
387- for (jj = 0; jj < 4; jj++)
388- for (jb = 0; jb < nb; jb++)
441+ for (i4 = 0; i4 < 4; i4++)
442+ for (jb = 0; jb < nb; jb++)
443+ for (j4 = 0; j4 < 4; j4++)
389444 for (ie = 0; ie < ne; ie++)
390- w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
445+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
391446 }
392447
393448 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
394449
395450 if (v > thr) {
396- for (ii = 0; ii < 4; ii++) {
397- ei2[ii] = 0.0;
398- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
399- for (jj = 0; jj < 4; jj++) {
400- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
451+ for (j4 = 0; j4 < 4; j4++) {
452+ ei2[j4] = 0.0;
453+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
454+ }
455+ for (i4 = 0; i4 < 4; i4++)
456+ for (j4 = 0; j4 < 4; j4++) {
457+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
401458 for (jb = 0; jb < nb; jb++)
402- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
459+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
403460 }
404- }
405461 libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
406- for (ii = 0; ii < 4; ii++)
407- for (jj = 0; jj < 4; jj++)
408- for (jb = 0; jb < nb; jb++)
462+ for (i4 = 0; i4 < 4; i4++)
463+ for (jb = 0; jb < nb; jb++)
464+ for (j4 = 0; j4 < 4; j4++)
409465 for (ie = 0; ie < ne; ie++)
410- w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
466+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
411467 }
412468 }
413469 else if (e[3] <= 0.0) {
414470
415- for (ii = 0; ii < 4; ii++) {
416- ei2[ii] = ej1[ii][ib];
417- for (jb = 0; jb < nb; jb++)
418- ej2[ii][jb] = ej1[ii][jb];
419- }
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+ }
420476 libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
421-
422- for (ii = 0; ii < 4; ii++)
477+ for (i4 = 0; i4 < 4; i4++)
423478 for (jb = 0; jb < nb; jb++)
424479 for (ie = 0; ie < ne; ie++)
425- w1[ib][jb][ie][ii] += w2[jb][ie][ii];
480+ w1[ib][i4][jb][ie] += w2[jb][i4][ie];
426481 }
427482 else {
428483 continue;
429484 }
430485 }
431- for (ii = 0; ii < 20; ii++)
432- for (jj = 0; jj < 4; jj++)
433- for (ib = 0; ib < nb; ib++)
486+ for (i4 = 0; i4 < 20; i4++)
487+ for (ib = 0; ib < nb; ib++)
488+ for (j4 = 0; j4 < 4; j4++)
434489 for (jb = 0; jb < nb; jb++)
435490 for (ie = 0; ie < ne; ie++)
436- wght[ikv[it][ii]][ib][jb][ie] += w1[ib][jb][ie][jj] * wlsm[jj][ii];
491+ wght[ikv[it][i20]][ib][jb][ie] += wlsm[i20][j4] * w1[ib][j4][jb][ie];
437492 }
438493 for (ik = 0; ik < nk; ik++)
439494 for (ib = 0; ib < nb; ib++)
--- a/src/c/libtetrabz_occ.c
+++ b/src/c/libtetrabz_occ.c
@@ -36,24 +36,24 @@ void occ(
3636 Main SUBROUTINE for occupation : Theta(EF - E1)
3737 */
3838 int it, ik, ib, ii, jj, ** ikv, indx[4];
39- double wlsm[4][20], ** ei1, e[4], ** w1, v, tsmall[4][4];
39+ double wlsm[20][4], ** ei1, e[4], ** w1, v, tsmall[4][4];
4040
41- ikv = (int**)calloc(6 * nk, sizeof(int*));
42- ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
41+ ikv = (int**)malloc(6 * nk * sizeof(int*));
42+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
4343 for (ik = 0; ik < 6 * nk; ik++) {
4444 ikv[ik] = ikv[0] + ik * 20;
4545 }
4646
47- ei1 = (double**)calloc(4, sizeof(double*));
48- ei1[0] = (double*)calloc(4 * nb, sizeof(double));
47+ ei1 = (double**)malloc(4 * sizeof(double*));
48+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
4949 for (ii = 0; ii < 4; ii++) {
5050 ei1[ii] = ei1[0] + ii * nb;
5151 }
5252
53- w1 = (double**)calloc(nb, sizeof(double*));
54- w1[0] = (double*)calloc(nb * 4, sizeof(double));
53+ w1 = (double**)malloc(nb * sizeof(double*));
54+ w1[0] = (double*)malloc(nb * 4 * sizeof(double));
5555 for (ib = 0; ib < nb; ib++) {
56- w1[ii] = w1[0] + ib * 4;
56+ w1[ib] = w1[0] + ib * 4;
5757 }
5858
5959 libtetrabz_initialize(ng, bvec, wlsm, ikv);
@@ -64,13 +64,13 @@ void occ(
6464
6565 for (it = 0; it < 6 * nk; it++) {
6666
67- for (ii = 0; ii < 4; ii++) {
68- for (ib = 0; ib < nb; ib++) {
67+ for (ii = 0; ii < 4; ii++)
68+ for (ib = 0; ib < nb; ib++)
6969 ei1[ii][ib] = 0.0;
70- for (jj = 0; jj < 20; jj++)
71- ei1[ii][ib] += wlsm[ii][jj] * eig[ikv[it][jj]][ib];
72- }
73- }
70+ for (jj = 0; jj < 20; jj++)
71+ for (ii = 0; ii < 4; ii++)
72+ for (ib = 0; ib < nb; ib++)
73+ ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
7474
7575 for (ib = 0; ib < nb; ib++)
7676 for (ii = 0; ii < 4; ii++)
@@ -85,54 +85,54 @@ void occ(
8585 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
8686 for (ii = 0; ii < 4; ii++)
8787 for (jj = 0; jj < 4; jj++)
88- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
88+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
8989 }
9090 else if (e[1] <= 0.0 && 0.0 < e[2]) {
9191
9292 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
9393 for (ii = 0; ii < 4; ii++)
9494 for (jj = 0; jj < 4; jj++)
95- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
95+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
9696
9797 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
9898 for (ii = 0; ii < 4; ii++)
9999 for (jj = 0; jj < 4; jj++)
100- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
100+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
101101
102102 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
103103 for (ii = 0; ii < 4; ii++)
104104 for (jj = 0; jj < 4; jj++)
105- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
105+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
106106 }
107107 else if (e[2] <= 0.0 && 0.0 < e[3]) {
108108
109109 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
110110 for (ii = 0; ii < 4; ii++)
111111 for (jj = 0; jj < 4; jj++)
112- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
112+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
113113
114114 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
115115 for (ii = 0; ii < 4; ii++)
116116 for (jj = 0; jj < 4; jj++)
117- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
117+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
118118
119119 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
120120 for (ii = 0; ii < 4; ii++)
121121 for (jj = 0; jj < 4; jj++)
122- w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
122+ w1[ib][ii] += v * tsmall[ii][jj] * 0.25;
123123 }
124124 else if (e[3] <= 0.0) {
125125 for (ii = 0; ii < 4; ii++)
126- w1[ib][ii] = 0.25;
126+ w1[ib][ii] += 0.25;
127127 }
128128 else {
129129 continue;
130130 }
131131 }
132132 for (ii = 0; ii < 20; ii++)
133- for (jj = 0; jj < 4; jj++)
134- for (ib = 0; ib < nb; ib++)
135- wght[ikv[it][ii]][ib] += w1[ib][jj] * wlsm[jj][ii];
133+ for (ib = 0; ib < nb; ib++)
134+ for (jj = 0; jj < 4; jj++)
135+ wght[ikv[it][ii]][ib] += wlsm[ii][jj] * w1[ib][jj];
136136 }
137137 for (ik = 0; ik < nk; ik++)
138138 for (ib = 0; ib < nb; ib++)
@@ -160,9 +160,15 @@ void fermieng(
160160 /*
161161 Calculate Fermi energy
162162 */
163- int maxiter, ik, ib;
163+ int maxiter, ik, ib, iter;
164164 double eps, elw, eup, **eig2, sumkmid;
165165
166+ eig2 = (double**)malloc(nk * sizeof(double*));
167+ eig2[0] = (double*)malloc(nk * nb * sizeof(double));
168+ for (ik = 0; ik < nk; ik++) {
169+ eig2[ik] = eig2[0] + ik * nb;
170+ }
171+
166172 maxiter = 300;
167173 eps = 1.0e-10;
168174
@@ -177,7 +183,7 @@ void fermieng(
177183 /*
178184 Bisection method
179185 */
180- for (*iteration = 0; *iteration < maxiter; *iteration++) {
186+ for (iter = 0; iter < maxiter; iter++) {
181187
182188 *ef = (eup + elw) * 0.5;
183189 /*
@@ -197,11 +203,15 @@ void fermieng(
197203 /*
198204 convergence check
199205 */
200- if (fabs(sumkmid - nelec) < eps)
206+ if (fabs(sumkmid - nelec) < eps) {
207+ *iteration = iter;
208+ free(eig2);
201209 return;
210+ }
202211 else if (sumkmid < nelec)
203212 elw = *ef;
204213 else
205214 eup = *ef;
206215 }
216+ free(eig2);
207217 }
--- a/src/c/libtetrabz_polcmplx.c
+++ b/src/c/libtetrabz_polcmplx.c
@@ -22,463 +22,287 @@
2222 */
2323 #include "libtetrabz_common.h"
2424 #include <math.h>
25+#include <stdlib.h>
2526
26-void polcmplx(
27- int *ng,
28- int nk,
29- int nb,
30- int ne,
31- double **bvec,
32- double **eig1,
33- double **eig2,
34- double **e0,
35- double *****wght) {
27+void libtetrabz_polcmplx_1234(
28+ double g1,
29+ double g2,
30+ double g3,
31+ double g4,
32+ double* w
33+) {
3634 /*
37- Main SUBROUTINE for Polarization (Imaginary axis) : Theta(- E1) * Theta(E2) / (E2 - E1 - iw)
35+ 1, Different each other
36+ */
37+ double w2, w3, w4;
38+ /*
39+ Real
40+ */
41+ w2 = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1)) +
42+ (g2 * g2 - 3.0) * g2 * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
43+ w2 = -2.0 * (g2 * g2 - 1.0) + w2 / (g2 - g1);
44+ w2 = w2 / (g2 - g1);
45+ w3 = 2.0 * (3.0 * g3 * g3 - 1.0) * (atan(g3) - atan(g1)) +
46+ (g3 * g3 - 3.0) * g3 * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
47+ w3 = -2.0 * (g3 * g3 - 1.0) + w3 / (g3 - g1);
48+ w3 = w3 / (g3 - g1);
49+ w4 = 2.0 * (3.0 * g4 * g4 - 1.0) * (atan(g4) - atan(g1)) +
50+ (g4 * g4 - 3.0) * g4 * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
51+ w4 = -2.0 * (g4 * g4 - 1.0) + w4 / (g4 - g1);
52+ w4 = w4 / (g4 - g1);
53+ w2 = (w2 - w3) / (g2 - g3);
54+ w4 = (w4 - w3) / (g4 - g3);
55+ w[0] = (w4 - w2) / (2.0 * (g4 - g2));
56+ /*
57+ Imaginal
3858 */
39- int it, ik, ie, ib, ii, jj, kk, jb, **ikv, indx[4];
40- double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], *****w1, ****w2, v, tsmall[4][4], thr;
59+ w2 = 2.0 * (3.0 - g2 * g2) * g2 * (atan(g2) - atan(g1)) +
60+ (3.0 * g2 * g2 - 1.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
61+ w2 = 4.0 * g2 - w2 / (g2 - g1);
62+ w2 = w2 / (g2 - g1);
63+ w3 = 2.0 * (3.0 - g3 * g3) * g3 * (atan(g3) - atan(g1)) +
64+ (3.0 * g3 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
65+ w3 = 4.0 * g3 - w3 / (g3 - g1);
66+ w3 = w3 / (g3 - g1);
67+ w4 = 2.0 * (3.0 - g4 * g4) * g4 * (atan(g4) - atan(g1)) +
68+ (3.0 * g4 * g4 - 1.0) * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
69+ w4 = 4.0 * g4 - w4 / (g4 - g1);
70+ w4 = w4 / (g4 - g1);
71+ w2 = (w2 - w3) / (g2 - g3);
72+ w4 = (w4 - w3) / (g4 - g3);
73+ w[1] = (w4 - w2) / (2.0 * (g4 - g2));
74+}
4175
42- libtetrabz_initialize(ng, bvec, wlsm, ikv);
4376
44- for (ik = 0; ik < nk; ik++)
45- for (ib = 0; ib < nb; ib++)
46- for (jb = 0; jb < nb; jb++)
47- for (ie = 0; ie < ne; ie++)
48- for (jj = 0; jj < 2; jj++)
49- wght[ik][ib][jb][ie][jj] = 0.0;
77+void libtetrabz_polcmplx_1231(
78+ double g1,
79+ double g2,
80+ double g3,
81+ double* w
82+) {
83+ /*
84+ 2, g4 = g1
85+ */
86+ double w2, w3;
87+ /*
88+ Real
89+ */
90+ w2 = 2.0 * (-1.0 + 3.0 * g2 * g2) * (atan(g2) - atan(g1))
91+ + g2 * (-3.0 + g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
92+ w2 = 2.0 * (1.0 - g2 * g2) + w2 / (g2 - g1);
93+ w2 = -g1 + w2 / (g2 - g1);
94+ w2 = w2 / (g2 - g1);
95+ w3 = 2.0 * (-1.0 + 3.0 * g3 * g3) * (atan(g3) - atan(g1))
96+ + g3 * (-3.0 + g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
97+ w3 = 2.0 * (1 - g3 * g3) + w3 / (g3 - g1);
98+ w3 = -g1 + w3 / (g3 - g1);
99+ w3 = w3 / (g3 - g1);
100+ w[0] = (w3 - w2) / (2.0 * (g3 - g2));
101+ /*
102+ Imaginal
103+ */
104+ w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
105+ (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
106+ w2 = 4.0 * g2 - w2 / (g2 - g1);
107+ w2 = 1 + w2 / (g2 - g1);
108+ w2 = w2 / (g2 - g1);
109+ w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
110+ (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
111+ w3 = 4.0 * g3 - w3 / (g3 - g1);
112+ w3 = 1 + w3 / (g3 - g1);
113+ w3 = w3 / (g3 - g1);
114+ w[1] = (w3 - w2) / (2.0 * (g3 - g2));
115+}
50116
51- thr = 1.0e-8;
52117
53- for (it = 0; it < 6 * nk; it++) {
118+void libtetrabz_polcmplx_1233(
119+ double g1,
120+ double g2,
121+ double g3,
122+ double* w
123+) {
124+ /*
125+ 3, g4 = g3
126+ */
127+ double w2, w3;
128+ /*
129+ Real
130+ */
131+ w2 = 2.0 * (1.0 - 3.0 * g2 * g2) * (atan(g2) - atan(g1)) +
132+ g2 * (3.0 - g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
133+ w2 = 2.0 * (1 - g2 * g2) - w2 / (g2 - g1);
134+ w2 = w2 / (g2 - g1);
135+ w3 = 2.0 * (1.0 - 3.0 * g3 * g3) * (atan(g3) - atan(g1)) +
136+ g3 * (3.0 - g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
137+ w3 = 2.0 * (1 - g3 * g3) - w3 / (g3 - g1);
138+ w3 = w3 / (g3 - g1);
139+ w2 = (w3 - w2) / (g3 - g2);
140+ w3 = 4.0 * (1.0 - 3.0 * g1 * g3) * (atan(g3) - atan(g1))
141+ + (3.0 * g1 + 3.0 * g3 - 3.0 * g1 * g3 * g3 + g3 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
142+ w3 = -4.0 * (1.0 - g1 * g1) + w3 / (g3 - g1);
143+ w3 = 4.0 * g1 + w3 / (g3 - g1);
144+ w3 = w3 / (g3 - g1);
145+ w[0] = (w3 - w2) / (2.0 * (g3 - g2));
146+ /*
147+ Imaginal
148+ */
149+ w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
150+ (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
151+ w2 = 4.0 * g2 - w2 / (g2 - g1);
152+ w2 = w2 / (g2 - g1);
153+ w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
154+ (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
155+ w3 = 4.0 * g3 - w3 / (g3 - g1);
156+ w3 = w3 / (g3 - g1);
157+ w2 = (w3 - w2) / (g3 - g2);
158+ w3 = (3.0 * g1 - 3.0 * g1 * g3 * g3 + 3.0 * g3 + g3 * g3 * g3) * (atan(g3) - atan(g1))
159+ + (3.0 * g1 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
160+ w3 = w3 / (g3 - g1) - 4.0 * g1;
161+ w3 = w3 / (g3 - g1) - 2.0;
162+ w3 = (2.0 * w3) / (g3 - g1);
163+ w[1] = (w3 - w2) / (2.0 * (g3 - g2));
164+}
54165
55- for (ii = 0; ii < 4; ii++) {
56- for (ib = 0; ib < nb; ib++) {
57- ei1[ii][ib] = 0.0;
58- ej1[ii][ib] = 0.0;
59- }
60- }
61- for (jj = 0; jj < 20; jj++) {
62- for (ib = 0; ib < nb; ib++) {
63- for (ii = 0; ii < 4; ii++) {
64- ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
65- ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
66- }
67- }
68- }
166+void libtetrabz_polcmplx_1221(
167+ double g1,
168+ double g2,
169+ double* w
170+) {
171+ /*
172+ 4, g4 = g1 and g3 = g2
173+ */
174+ /*
175+ Real
176+ */
177+ w[0] = -2.0 * (-1.0 + 2.0 * g1 * g2 + g2 * g2) * (atan(g2) - atan(g1))
178+ + (g1 + 2.0 * g2 - g1 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
179+ w[0] = 2.0 * (-1.0 + g1 * g1) + w[0] / (g2 - g1);
180+ w[0] = 3.0 * g1 + w[0] / (g2 - g1);
181+ w[0] = 2.0 + (3.0 * w[0]) / (g2 - g1);
182+ w[0] = w[0] / (2.0 * (g2 - g1));
183+ /*
184+ Imaginal
185+ */
186+ w[1] = 2.0 * (g1 + 2.0 * g2 - g1 * g2 * g2) * (atan(g2) - atan(g1))
187+ + (-1.0 + 2.0 * g1 * g2 + g2 * g2) * log((1 + g2 * g2) / (1 + g1 * g1));
188+ w[1] = -4.0 * g1 + w[1] / (g2 - g1);
189+ w[1] = -3.0 + w[1] / (g2 - g1);
190+ w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
191+}
69192
70- for (ib = 0; ib < nb; ib++) {
71193
72- for (jb = 0; jb < nb; jb++)
73- for (ii = 0; ii < 4; ii++)
74- for (ie = 0; ie < ne; ie++)
75- for (jj = 0; jj < 2; jj++)
76- w1[ib][jb][ie][ii][jj] = 0.0;
194+void libtetrabz_polcmplx_1222(
195+ double g1,
196+ double g2,
197+ double* w
198+) {
199+ /*
200+ 5, g4 = g3 = g2
201+ */
202+ /*
203+ Real
204+ */
205+ w[0] = 2.0 * (-1.0 + g1 * g1 + 2.0 * g1 * g2) * (atan(g2) - atan(g1))
206+ + (-2.0 * g1 - g2 + g1 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
207+ w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
208+ w[0] = g1 - w[0] / (g2 - g1);
209+ w[0] = 1.0 - (3.0 * w[0]) / (g2 - g1);
210+ w[0] = w[0] / (2.0 * (g2 - g1));
211+ /*
212+ Imaginal
213+ */
214+ w[1] = 2.0 * (-2.0 * g1 - g2 + g1 * g1 * g2) * (atan(g2) - atan(g1))
215+ + (1.0 - g1 * g1 - 2.0 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
216+ w[1] = 4.0 * g1 + w[1] / (g2 - g1);
217+ w[1] = 1.0 + w[1] / (g2 - g1);
218+ w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
219+}
77220
78- for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
79- eig_sort(4, e, indx);
80221
81- if (e[0] <= 0.0 && 0.0 < e[1]) {
222+void libtetrabz_polcmplx_1211(
223+ double g1,
224+ double g2,
225+ double* w
226+) {
227+ /*
228+ 6, g4 = g3 = g1
229+ */
230+ /*
231+ Real
232+ */
233+ w[0] = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1))
234+ + g2 * (g2 * g2 - 3.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
235+ w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
236+ w[0] = -5.0 * g1 + w[0] / (g2 - g1);
237+ w[0] = -11.0 + (3.0 * w[0]) / (g2 - g1);
238+ w[0] = w[0] / (6.0 * (g2 - g1));
239+ /*
240+ Imaginal
241+ */
242+ w[1] = 2.0 * g2 * (-3.0 + g2 * g2) * (atan(g2) - atan(g1))
243+ + (1.0 - 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
244+ w[1] = 4.0 * g2 + w[1] / (g2 - g1);
245+ w[1] = 1.0 + w[1] / (g2 - g1);
246+ w[1] = w[1] / (2.0 * (g2 - g1) * (g2 - g1));
247+}
82248
83- libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
249+void libtetrabz_polcmplx3(
250+ int ne,
251+ double** e0,
252+ double de[4],
253+ double*** w1
254+) {
255+ /*
256+ Tetrahedron method for delta(om - ep + e)
257+ */
258+ int i4, ir, indx[4], ie;
259+ double e[4], x[4], thr, w2[4][2];
84260
85- if (v > thr) {
86- for (ii = 0; ii < 4; ii++) {
87- ei2[ii] = 0.0;
88- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
89- for (jj = 0; jj < 4; jj++) {
90- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
91- for (jb = 0; jb < nb; jb++)
92- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
93- }
94- }
95- libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
96- for (ii = 0; ii < 4; ii++)
97- for (jj = 0; jj < 4; jj++)
98- for (jb = 0; jb < nb; jb++)
99- for (ie = 0; ie < ne; ie++)
100- for (kk = 0; kk < 4; kk++)
101- w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
102- }
103- }
104- else if (e[1] <= 0.0 && 0.0 < e[2]) {
261+ for (i4 = 0; i4 < 3; i4++) e[i4] = de[i4];
262+ eig_sort(4, e, indx);
105263
106- libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
264+ for (ie = 0; ie < ne; ie++) {
265+ /*
266+ I am not sure which one is better.
267+ The former is more stable.
268+ The latter is more accurate ?
269+ */
270+ for (i4 = 0; i4 < 4; i4++)
271+ for (ir = 0; ir < 2; ir++)
272+ w1[i4][ie][ir] = 0.25 / (de[i4] + e0[ie][ir]);
107273
108- if (v > thr) {
109- for (ii = 0; ii < 4; ii++) {
110- ei2[ii] = 0.0;
111- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
112- for (jj = 0; jj < 4; jj++) {
113- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
114- for (jb = 0; jb < nb; jb++)
115- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
116- }
117- }
118- libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
119- for (ii = 0; ii < 4; ii++)
120- for (jj = 0; jj < 4; jj++)
121- for (jb = 0; jb < nb; jb++)
122- for (ie = 0; ie < ne; ie++)
123- for (kk = 0; kk < 4; kk++)
124- w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
125- }
274+ continue;
126275
127- libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
276+ for (i4 = 0; i4 < 4; i4++)
277+ x[i4] = (e[i4] + e0[ie][0]) / e0[ie][1];
278+ /* thr = maxval(de(1:4)) * 1d-3*/
279+ thr = fabs(x[0]);
280+ for (i4 = 0; i4 < 4; i4++)
281+ if (thr < fabs(x[i4])) thr = fabs(x[i4]);
282+ thr = fmax(1.0e-3, thr * 1.0e-2);
128283
129- if (v > thr) {
130- for (ii = 0; ii < 4; ii++) {
131- ei2[ii] = 0.0;
132- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
133- for (jj = 0; jj < 4; jj++) {
134- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
135- for (jb = 0; jb < nb; jb++)
136- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
137- }
284+ if (fabs(x[3] - x[2]) < thr) {
285+ if (fabs(x[3] - x[1]) < thr) {
286+ if (fabs(x[3] - x[0]) < thr) {
287+ /*
288+ e[3] = e[2] = e[1] = e[0]
289+ */
290+ w2[3][0] = 0.25 * x[3] / (1.0 + x[3] * x[3]);
291+ w2[3][1] = 0.25 / (1.0 + x[3] * x[3]);
292+ for (ir = 0; ir < 2; ir++) {
293+ w2[2][ir] = w2[3][ir];
294+ w2[1][ir] = w2[3][ir];
295+ w2[0][ir] = w2[3][ir];
138296 }
139- libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
140- for (ii = 0; ii < 4; ii++)
141- for (jj = 0; jj < 4; jj++)
142- for (jb = 0; jb < nb; jb++)
143- for (ie = 0; ie < ne; ie++)
144- for (kk = 0; kk < 4; kk++)
145- w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
146297 }
147-
148- libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
149-
150- if (v > thr) {
151- for (ii = 0; ii < 4; ii++) {
152- ei2[ii] = 0.0;
153- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
154- for (jj = 0; jj < 4; jj++) {
155- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
156- for (jb = 0; jb < nb; jb++)
157- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
158- }
159- }
160- libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
161- for (ii = 0; ii < 4; ii++)
162- for (jj = 0; jj < 4; jj++)
163- for (jb = 0; jb < nb; jb++)
164- for (ie = 0; ie < ne; ie++)
165- for (kk = 0; kk < 4; kk++)
166- w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
167- }
168- }
169- else if (e[2] <= 0.0 && 0.0 < e[3]) {
170-
171- libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
172-
173- if (v > thr) {
174- for (ii = 0; ii < 4; ii++) {
175- ei2[ii] = 0.0;
176- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
177- for (jj = 0; jj < 4; jj++) {
178- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
179- for (jb = 0; jb < nb; jb++)
180- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
181- }
182- }
183- libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
184- for (ii = 0; ii < 4; ii++)
185- for (jj = 0; jj < 4; jj++)
186- for (jb = 0; jb < nb; jb++)
187- for (ie = 0; ie < ne; ie++)
188- for (kk = 0; kk < 4; kk++)
189- w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
190- }
191-
192- libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
193-
194- if (v > thr) {
195- for (ii = 0; ii < 4; ii++) {
196- ei2[ii] = 0.0;
197- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
198- for (jj = 0; jj < 4; jj++) {
199- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
200- for (jb = 0; jb < nb; jb++)
201- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
202- }
203- }
204- libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
205- for (ii = 0; ii < 4; ii++)
206- for (jj = 0; jj < 4; jj++)
207- for (jb = 0; jb < nb; jb++)
208- for (ie = 0; ie < ne; ie++)
209- for (kk = 0; kk < 4; kk++)
210- w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
211- }
212-
213- libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
214-
215- if (v > thr) {
216- for (ii = 0; ii < 4; ii++) {
217- ei2[ii] = 0.0;
218- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
219- for (jj = 0; jj < 4; jj++) {
220- ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
221- for (jb = 0; jb < nb; jb++)
222- ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
223- }
224- }
225- libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
226- for (ii = 0; ii < 4; ii++)
227- for (jj = 0; jj < 4; jj++)
228- for (jb = 0; jb < nb; jb++)
229- for (ie = 0; ie < ne; ie++)
230- for (kk = 0; kk < 4; kk++)
231- w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
232- }
233- }
234- else if (e[3] <= 0.0) {
235-
236- for (ii = 0; ii < 4; ii++) {
237- ei2[ii] = ej1[ii][ib];
238- for (jb = 0; jb < nb; jb++) ej2[ii][jb] = ej1[ii][jb];
239- }
240- libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
241- for (ii = 0; ii < 4; ii++)
242- for (jb = 0; jb < nb; jb++)
243- for (ie = 0; ie < ne; ie++)
244- for (kk = 0; kk < 4; kk++)
245- w1[ib][jb][ie][kk][ii] += w2[jb][ie][kk][ii];
246- }
247- else {
248- continue;
249- }
250- }
251- for (ii = 0; ii < 20; ii++)
252- for (jj = 0; jj < 4; jj++)
253- for (ib = 0; ib < nb; ib++)
254- for (jb = 0; jb < nb; jb++)
255- for (ie = 0; ie < ne; ie++)
256- for (kk = 0; kk < 2; kk++)
257- wght[ikv[it][ii]][ib][jb][ie][kk] += w1[ib][jb][ie][kk][jj] * wlsm[jj][ii];
258- }
259- for (ik = 0; ik < nk; ik++)
260- for (ib = 0; ib < nb; ib++)
261- for (jb = 0; jb < nb; jb++)
262- for (ie = 0; ie < ne; ie++)
263- for (ii = 0; ii < 2; ii++)
264- wght[ik][ib][jb][ie][ii] /= (6.0 * (double) nk);
265-}
266-
267-void libtetrabz_polcmplx2(
268- int nb,
269- int ne,
270- double *e0,
271- double *ei1,
272- double **ej1,
273- double ****w1
274-) {
275- /*
276- Tetrahedron method for theta( - E2)
277- */
278- int ib, ii, jj, kk, ie, indx[4];
279- double e[4], tsmall[4][4], v, de[4], thr, ***w2;
280-
281- thr = 1.0e-8;
282-
283- for (ib = 0; ib < nb; ib++) {
284-
285- for (ie = 0; ie < ne; ie++)
286- for (ii = 0; ii < 4; ii++)
287- for (jj = 0; jj < 2; jj++)
288- w1[ib][ie][jj][ii] = 0.0;
289-
290- for (ii = 0; ii < 4; ii++) e[ii] = -ej1[ii][ib];
291- eig_sort(4, e, indx);
292-
293- if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
294-
295- libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
296-
297- if (v > thr) {
298- for (ii = 0; ii < 4; ii++) {
299- de[ii] = 0.0;
300- for (jj = 0; jj < 4; jj++)
301- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
302- }
303- libtetrabz_polcmplx3(ne, e0, de, w2);
304- for (ii = 0; ii < 4; ii++)
305- for (jj = 0; jj < 4; jj++)
306- for (ie = 0; ie < ne; ie++)
307- for (kk = 0; kk < 2; kk++)
308- w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
309- }
310- }
311- else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
312-
313- libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
314-
315- if (v > thr) {
316- for (ii = 0; ii < 4; ii++) {
317- de[ii] = 0.0;
318- for (jj = 0; jj < 4; jj++)
319- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
320- }
321- libtetrabz_polcmplx3(ne, e0, de, w2);
322- for (ii = 0; ii < 4; ii++)
323- for (jj = 0; jj < 4; jj++)
324- for (ie = 0; ie < ne; ie++)
325- for (kk = 0; kk < 2; kk++)
326- w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
327- }
328-
329- libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
330-
331- if (v > thr) {
332- for (ii = 0; ii < 4; ii++) {
333- de[ii] = 0.0;
334- for (jj = 0; jj < 4; jj++)
335- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
336- }
337- libtetrabz_polcmplx3(ne, e0, de, w2);
338- for (ii = 0; ii < 4; ii++)
339- for (jj = 0; jj < 4; jj++)
340- for (ie = 0; ie < ne; ie++)
341- for (kk = 0; kk < 2; kk++)
342- w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
343- }
344-
345- libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
346-
347- if (v > thr) {
348- for (ii = 0; ii < 4; ii++) {
349- de[ii] = 0.0;
350- for (jj = 0; jj < 4; jj++)
351- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
352- }
353- libtetrabz_polcmplx3(ne, e0, de, w2);
354- for (ii = 0; ii < 4; ii++)
355- for (jj = 0; jj < 4; jj++)
356- for (ie = 0; ie < ne; ie++)
357- for (kk = 0; kk < 2; kk++)
358- w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
359- }
360- }
361- else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
362-
363- libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
364-
365- if (v > thr) {
366- for (ii = 0; ii < 4; ii++) {
367- de[ii] = 0.0;
368- for (jj = 0; jj < 4; jj++)
369- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
370- }
371- libtetrabz_polcmplx3(ne, e0, de, w2);
372- for (ii = 0; ii < 4; ii++)
373- for (jj = 0; jj < 4; jj++)
374- for (ie = 0; ie < ne; ie++)
375- for (kk = 0; kk < 2; kk++)
376- w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
377- }
378-
379- libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
380-
381- if (v > thr) {
382- for (ii = 0; ii < 4; ii++) {
383- de[ii] = 0.0;
384- for (jj = 0; jj < 4; jj++)
385- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
386- }
387- libtetrabz_polcmplx3(ne, e0, de, w2);
388- for (ii = 0; ii < 4; ii++)
389- for (jj = 0; jj < 4; jj++)
390- for (ie = 0; ie < ne; ie++)
391- for (kk = 0; kk < 2; kk++)
392- w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
393- }
394-
395- libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
396-
397- if (v > thr) {
398- for (ii = 0; ii < 4; ii++) {
399- de[ii] = 0.0;
400- for (jj = 0; jj < 4; jj++)
401- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
402- }
403- libtetrabz_polcmplx3(ne, e0, de, w2);
404- for (ii = 0; ii < 4; ii++)
405- for (jj = 0; jj < 4; jj++)
406- for (ie = 0; ie < ne; ie++)
407- for (kk = 0; kk < 2; kk++)
408- w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
409- }
410- }
411- else if (e[3] <= 0.0) {
412- for (ii = 0; ii < 4; ii++) {
413- de[ii] = tsmall[ii][jj] * (ej1[ii][ib] - ei1[ii]);
414- }
415- libtetrabz_polcmplx3(ne, e0, de, w2);
416- for (ii = 0; ii < 4; ii++)
417- for (ie = 0; ie < ne; ie++)
418- for (kk = 0; kk < 2; kk++)
419- w1[ib][ie][kk][ii] += w2[ie][kk][ii];
420- }
421- }
422-}
423-
424-
425-void libtetrabz_polcmplx3(
426- int ne,
427- double **e0,
428- double de[4],
429- double ***w1
430-) {
431- /*
432- Tetrahedron method for delta(om - ep + e)
433- */
434- int ii, jj, indx[4], ie;
435- double e[4], x[4], thr, w2[4][2];
436-
437- for (ii = 0; ii < 3; ii++) e[ii] = de[ii];
438- eig_sort(4, e, indx);
439-
440- for (ie = 0; ie < ne; ie++) {
441- /*
442- I am not sure which one is better.
443- The former is more stable.
444- The latter is more accurate ?
445- */
446- for (ii = 0; ii < 4; ii++)
447- for (jj = 0; jj < 2; jj++)
448- w1[ie][jj][ii] = 0.25 / (de[ii] + e0[ie][jj]);
449-
450- continue;
451-
452- for (ii = 0; ii < 4; ii++)
453- x[ii] = (e[ii] + e0[ie][0]) / e0[ie][1];
454- /* thr = maxval(de(1:4)) * 1d-3*/
455- thr = fabs(x[0]);
456- for(ii=0;ii<4;ii++)
457- if(thr < fabs(x[ii])) thr = fabs(x[ii]);
458- thr = max(1.0e-3, thr * 1.0e-2);
459-
460- if (fabs(x[3] - x[2]) < thr) {
461- if (fabs(x[3] - x[1]) < thr) {
462- if (fabs(x[3] - x[0]) < thr) {
463- /*
464- e[3] = e[2] = e[1] = e[0]
465- */
466- w2[3][0] = 0.25 * x[3] / (1.0 + x[3] *x[3]);
467- w2[3][1] = 0.25 / (1.0 + x[3] * x[3]);
468- for (ii = 0; ii < 2; ii++) {
469- w2[2][ii] = w2[3][ii];
470- w2[1][ii] = w2[3][ii];
471- w2[0][ii] = w2[3][ii];
472- }
473- }
474- else {
475- /*
476- e[3] = e[2] = e[1]
477- */
478- libtetrabz_polcmplx_1211(x[3], x[0], w2[3]);
479- for (ii = 0; ii < 2; ii++) {
480- w2[2][ii] = w2[3][ii];
481- w2[1][ii] = w2[3][ii];
298+ else {
299+ /*
300+ e[3] = e[2] = e[1]
301+ */
302+ libtetrabz_polcmplx_1211(x[3], x[0], w2[3]);
303+ for (ir = 0; ir < 2; ir++) {
304+ w2[2][ir] = w2[3][ir];
305+ w2[1][ir] = w2[3][ir];
482306 }
483307 libtetrabz_polcmplx_1222(x[0], x[3], w2[0]);
484308 /*
@@ -495,9 +319,9 @@ void libtetrabz_polcmplx3(
495319 e[3] = e[2], e[1] = e[0]
496320 */
497321 libtetrabz_polcmplx_1221(x[3], x[1], w2[3]);
498- for (ii = 0; ii < 2; ii++) w2[2][ii] = w2[3][ii];
322+ for (ir = 0; ir < 2; ir++) w2[2][ir] = w2[3][ir];
499323 libtetrabz_polcmplx_1221(x[1], x[3], w2[1]);
500- for (ii = 0; ii < 2; ii++) w2[0][ii] = w2[1][ii];
324+ for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[1][ir];
501325 /*
502326 IF(ANY(w2(1:2,1:4) < 0.0)):
503327 WRITE(*,*) ie
@@ -511,7 +335,7 @@ void libtetrabz_polcmplx3(
511335 e[3] = e[2]
512336 */
513337 libtetrabz_polcmplx_1231(x[3], x[0], x[1], w2[3]);
514- for (ii = 0; ii < 2; ii++) w2[2][ii] = w2[3][ii];
338+ for (ir = 0; ir < 2; ir++) w2[2][ir] = w2[3][ir];
515339 libtetrabz_polcmplx_1233(x[1], x[0], x[3], w2[1]);
516340 libtetrabz_polcmplx_1233(x[0], x[1], x[3], w2[0]);
517341 /*
@@ -530,8 +354,8 @@ void libtetrabz_polcmplx3(
530354 */
531355 libtetrabz_polcmplx_1222(x[3], x[2], w2[3]);
532356 libtetrabz_polcmplx_1211(x[2], x[3], w2[2]);
533- for (ii = 0; ii < 2; ii++) w2[1][ii] = w2[2][ii];
534- for (ii = 0; ii < 2; ii++) w2[0][ii] = w2[2][ii];
357+ for (ir = 0; ir < 2; ir++) w2[1][ir] = w2[2][ir];
358+ for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[2][ir];
535359 /*
536360 IF(ANY(w2(1:2,1:4) < 0.0)):
537361 WRITE(*,*) ie
@@ -546,7 +370,7 @@ void libtetrabz_polcmplx3(
546370 */
547371 libtetrabz_polcmplx_1233(x[3], x[0], x[2], w2[3]);
548372 libtetrabz_polcmplx_1231(x[2], x[0], x[3], w2[2]);
549- for (ii = 0; ii < 2; ii++) w2[1][ii] = w2[2][ii];
373+ for (ir = 0; ir < 2; ir++) w2[1][ir] = w2[2][ir];
550374 libtetrabz_polcmplx_1233(x[0], x[3], x[2], w2[0]);
551375 /*
552376 IF(ANY(w2(1:2,1:4) < 0.0)):
@@ -564,7 +388,7 @@ void libtetrabz_polcmplx3(
564388 libtetrabz_polcmplx_1233(x[3], x[2], x[1], w2[3]);
565389 libtetrabz_polcmplx_1233(x[2], x[3], x[1], w2[2]);
566390 libtetrabz_polcmplx_1231(x[1], x[2], x[3], w2[1]);
567- for (ii = 0; ii < 2; ii++) w2[0][ii] = w2[1][ii];
391+ for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[1][ir];
568392 /*
569393 IF(ANY(w2(1:2,1:4) < 0.0)):
570394 WRITE(*,*) ie
@@ -589,9 +413,9 @@ void libtetrabz_polcmplx3(
589413 STOP "weighting"
590414 */
591415 }
592- for (ii = 0; ii < 4; ii++) {
593- w1[ie][indx[ii]][0] = w2[ii][0] / e0[ie][1];
594- w1[ie][indx[ii]][1] = w2[ii][1] / (-e0[ie][1]);
416+ for (i4 = 0; i4 < 4; i4++) {
417+ w1[indx[i4]][ie][0] = w2[i4][0] / e0[ie][1];
418+ w1[indx[i4]][ie][1] = w2[i4][1] / (-e0[ie][1]);
595419 }
596420 }
597421 }
@@ -599,225 +423,484 @@ void libtetrabz_polcmplx3(
599423 // Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
600424 // for 0<x<1, 0<y<1-x, 0<z<1-x-y
601425
602-
603-void libtetrabz_polcmplx_1234(
604- double g1,
605- double g2,
606- double g3,
607- double g4,
608- double *w
609- ) {
610- /*
611- 1, Different each other
612- */
613- double w2, w3, w4;
614- /*
615- Real
616- */
617- w2 = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1)) +
618- (g2 * g2 - 3.0) * g2 * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
619- w2 = -2.0 * (g2 * g2 - 1.0) + w2 / (g2 - g1);
620- w2 = w2 / (g2 - g1);
621- w3 = 2.0 * (3.0 * g3 * g3 - 1.0) * (atan(g3) - atan(g1)) +
622- (g3 * g3 - 3.0) * g3 * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
623- w3 = -2.0 * (g3 * g3 - 1.0) + w3 / (g3 - g1);
624- w3 = w3 / (g3 - g1);
625- w4 = 2.0 * (3.0 * g4 * g4 - 1.0) * (atan(g4) - atan(g1)) +
626- (g4 * g4 - 3.0) * g4 * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
627- w4 = -2.0 * (g4 * g4 - 1.0) + w4 / (g4 - g1);
628- w4 = w4 / (g4 - g1);
629- w2 = (w2 - w3) / (g2 - g3);
630- w4 = (w4 - w3) / (g4 - g3);
631- w[0] = (w4 - w2) / (2.0 * (g4 - g2));
426+void libtetrabz_polcmplx2(
427+ int nb,
428+ int ne,
429+ double* *e0,
430+ double* ei1,
431+ double** ej1,
432+ double**** w1
433+) {
632434 /*
633- Imaginal
435+ Tetrahedron method for theta( - E2)
634436 */
635- w2 = 2.0 * (3.0 - g2 * g2) * g2 * (atan(g2) - atan(g1)) +
636- (3.0 * g2 * g2 - 1.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
637- w2 = 4.0 * g2 - w2 / (g2 - g1);
638- w2 = w2 / (g2 - g1);
639- w3 = 2.0 * (3.0 - g3 * g3) * g3 * (atan(g3) - atan(g1)) +
640- (3.0 * g3 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
641- w3 = 4.0 * g3 - w3 / (g3 - g1);
642- w3 = w3 / (g3 - g1);
643- w4 = 2.0 * (3.0 - g4 * g4) * g4 * (atan(g4) - atan(g1)) +
644- (3.0 * g4 * g4 - 1.0) * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
645- w4 = 4.0 * g4 - w4 / (g4 - g1);
646- w4 = w4 / (g4 - g1);
647- w2 = (w2 - w3) / (g2 - g3);
648- w4 = (w4 - w3) / (g4 - g3);
649- w[1] = (w4 - w2) / (2.0 * (g4 - g2));
650-}
437+ int ib, i4, j4, ir, ie, indx[4];
438+ double e[4], tsmall[4][4], v, de[4], thr, *** w2;
439+
440+ w2 = (double***)malloc(4 * sizeof(double**));
441+ w2[0] = (double**)malloc(4 * ne * sizeof(double*));
442+ w2[0][0] = (double*)malloc(4 * ne * 2 * sizeof(double));
443+ for (i4 = 0; i4 < 4; i4++) {
444+ w2[i4] = w2[0] + i4 * ne;
445+ for (ie = 0; ie < ne; ie++) {
446+ w2[i4][ie] = w2[0][0] + i4 * ne * 2 + ie * 2;
447+ }
448+ }
651449
450+ thr = 1.0e-8;
652451
653-void libtetrabz_polcmplx_1231(
654- double g1,
655- double g2,
656- double g3,
657- double *w
658- ) {
659- /*
660- 2, g4 = g1
661- */
662- double w2, w3;
663- /*
664- Real
665- */
666- w2 = 2.0 * (-1.0 + 3.0 * g2 * g2) * (atan(g2) - atan(g1))
667- + g2 * (-3.0 + g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
668- w2 = 2.0 * (1.0 - g2 * g2) + w2 / (g2 - g1);
669- w2 = -g1 + w2 / (g2 - g1);
670- w2 = w2 / (g2 - g1);
671- w3 = 2.0 * (-1.0 + 3.0 * g3 * g3) * (atan(g3) - atan(g1))
672- + g3 * (-3.0 + g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
673- w3 = 2.0 * (1 - g3 * g3) + w3 / (g3 - g1);
674- w3 = -g1 + w3 / (g3 - g1);
675- w3 = w3 / (g3 - g1);
676- w[0] = (w3 - w2) / (2.0 * (g3 - g2));
677- /*
678- Imaginal
679- */
680- w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
681- (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
682- w2 = 4.0 * g2 - w2 / (g2 - g1);
683- w2 = 1 + w2 / (g2 - g1);
684- w2 = w2 / (g2 - g1);
685- w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
686- (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
687- w3 = 4.0 * g3 - w3 / (g3 - g1);
688- w3 = 1 + w3 / (g3 - g1);
689- w3 = w3 / (g3 - g1);
690- w[1] = (w3 - w2) / (2.0 * (g3 - g2));
691-}
452+ for (ib = 0; ib < nb; ib++) {
692453
454+ for (i4 = 0; i4 < 4; i4++)
455+ for (ie = 0; ie < ne; ie++)
456+ for (ir = 0; ir < 2; ir++)
457+ w1[ib][i4][ie][ir] = 0.0;
693458
694-void libtetrabz_polcmplx_1233(
695- double g1,
696- double g2,
697- double g3,
698- double *w
699- ) {
700- /*
701- 3, g4 = g3
702- */
703- double w2, w3;
704- /*
705- Real
706- */
707- w2 = 2.0 * (1.0 - 3.0 * g2 * g2) * (atan(g2) - atan(g1)) +
708- g2 * (3.0 - g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
709- w2 = 2.0 * (1 - g2 * g2) - w2 / (g2 - g1);
710- w2 = w2 / (g2 - g1);
711- w3 = 2.0 * (1.0 - 3.0 * g3 * g3) * (atan(g3) - atan(g1)) +
712- g3 * (3.0 - g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
713- w3 = 2.0 * (1 - g3 * g3) - w3 / (g3 - g1);
714- w3 = w3 / (g3 - g1);
715- w2 = (w3 - w2) / (g3 - g2);
716- w3 = 4.0 * (1.0 - 3.0 * g1 * g3) * (atan(g3) - atan(g1))
717- + (3.0 * g1 + 3.0 * g3 - 3.0 * g1 * g3 * g3 + g3 * g3*g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
718- w3 = -4.0 * (1.0 - g1 * g1) + w3 / (g3 - g1);
719- w3 = 4.0 * g1 + w3 / (g3 - g1);
720- w3 = w3 / (g3 - g1);
721- w[0] = (w3 - w2) / (2.0 * (g3 - g2));
722- /*
723- Imaginal
724- */
725- w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
726- (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
727- w2 = 4.0 * g2 - w2 / (g2 - g1);
728- w2 = w2 / (g2 - g1);
729- w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
730- (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
731- w3 = 4.0 * g3 - w3 / (g3 - g1);
732- w3 = w3 / (g3 - g1);
733- w2 = (w3 - w2) / (g3 - g2);
734- w3 = (3.0 * g1 - 3.0 * g1 * g3 * g3 + 3.0 * g3 + g3 * g3*g3) * (atan(g3) - atan(g1))
735- + (3.0 * g1 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
736- w3 = w3 / (g3 - g1) - 4.0 * g1;
737- w3 = w3 / (g3 - g1) - 2.0;
738- w3 = (2.0 * w3) / (g3 - g1);
739- w[1] = (w3 - w2) / (2.0 * (g3 - g2));
740-}
459+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[i4][ib];
460+ eig_sort(4, e, indx);
741461
742-void libtetrabz_polcmplx_1221(
743- double g1,
744- double g2,
745- double *w
746- ) {
747- /*
748- 4, g4 = g1 and g3 = g2
749- */
750- /*
751- Real
752- */
753- w[0] = -2.0 * (-1.0 + 2.0 * g1 * g2 + g2 * g2) * (atan(g2) - atan(g1))
754- + (g1 + 2.0 * g2 - g1 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
755- w[0] = 2.0 * (-1.0 + g1 * g1) + w[0] / (g2 - g1);
756- w[0] = 3.0 * g1 + w[0] / (g2 - g1);
757- w[0] = 2.0 + (3.0 * w[0]) / (g2 - g1);
758- w[0] = w[0] / (2.0 * (g2 - g1));
759- /*
760- Imaginal
761- */
762- w[1] = 2.0 * (g1 + 2.0 * g2 - g1 * g2 * g2) * (atan(g2) - atan(g1))
763- + (-1.0 + 2.0 * g1 * g2 + g2 * g2) * log((1 + g2 * g2) / (1 + g1 * g1));
764- w[1] = -4.0 * g1 + w[1] / (g2 - g1);
765- w[1] = -3.0 + w[1] / (g2 - g1);
766- w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
767-}
462+ if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
768463
464+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
769465
770-void libtetrabz_polcmplx_1222(
771- double g1,
772- double g2,
773- double *w
774- ) {
775- /*
776- 5, g4 = g3 = g2
777- */
778- /*
779- Real
780- */
781- w[0] = 2.0 * (-1.0 + g1 * g1 + 2.0 * g1 * g2) * (atan(g2) - atan(g1))
782- + (-2.0 * g1 - g2 + g1 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
783- w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
784- w[0] = g1 - w[0] / (g2 - g1);
785- w[0] = 1.0 - (3.0 * w[0]) / (g2 - g1);
786- w[0] = w[0] / (2.0 * (g2 - g1));
787- /*
788- Imaginal
789- */
790- w[1] = 2.0 * (-2.0 * g1 - g2 + g1 * g1 * g2) * (atan(g2) - atan(g1))
791- + (1.0 - g1 * g1 - 2.0 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
792- w[1] = 4.0 * g1 + w[1] / (g2 - g1);
793- w[1] = 1.0 + w[1] / (g2 - g1);
794- w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
466+ if (v > thr) {
467+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
468+ for (i4 = 0; i4 < 4; i4++)
469+ for (j4 = 0; j4 < 4; j4++)
470+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
471+ libtetrabz_polcmplx3(ne, e0, de, w2);
472+ for (i4 = 0; i4 < 4; i4++)
473+ for (j4 = 0; j4 < 4; j4++)
474+ for (ie = 0; ie < ne; ie++)
475+ for (ir = 0; ir < 2; ir++)
476+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
477+ }
478+ }
479+ else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
480+
481+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
482+
483+ if (v > thr) {
484+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
485+ for (i4 = 0; i4 < 4; i4++)
486+ for (j4 = 0; j4 < 4; j4++)
487+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
488+ libtetrabz_polcmplx3(ne, e0, de, w2);
489+ for (i4 = 0; i4 < 4; i4++)
490+ for (j4 = 0; j4 < 4; j4++)
491+ for (ie = 0; ie < ne; ie++)
492+ for (ir = 0; ir < 2; ir++)
493+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
494+ }
495+
496+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
497+
498+ if (v > thr) {
499+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
500+ for (i4 = 0; i4 < 4; i4++)
501+ for (j4 = 0; j4 < 4; j4++)
502+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
503+ libtetrabz_polcmplx3(ne, e0, de, w2);
504+ for (i4 = 0; i4 < 4; i4++)
505+ for (j4 = 0; j4 < 4; j4++)
506+ for (ie = 0; ie < ne; ie++)
507+ for (ir = 0; ir < 2; ir++)
508+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
509+ }
510+
511+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
512+
513+ if (v > thr) {
514+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
515+ for (i4 = 0; i4 < 4; i4++)
516+ for (j4 = 0; j4 < 4; j4++)
517+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
518+ libtetrabz_polcmplx3(ne, e0, de, w2);
519+ for (i4 = 0; i4 < 4; i4++)
520+ for (j4 = 0; j4 < 4; j4++)
521+ for (ie = 0; ie < ne; ie++)
522+ for (ir = 0; ir < 2; ir++)
523+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
524+ }
525+ }
526+ else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
527+
528+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
529+
530+ if (v > thr) {
531+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
532+ for (i4 = 0; i4 < 4; i4++)
533+ for (j4 = 0; j4 < 4; j4++)
534+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
535+ libtetrabz_polcmplx3(ne, e0, de, w2);
536+ for (i4 = 0; i4 < 4; i4++)
537+ for (j4 = 0; j4 < 4; j4++)
538+ for (ie = 0; ie < ne; ie++)
539+ for (ir = 0; ir < 2; ir++)
540+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
541+ }
542+
543+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
544+
545+ if (v > thr) {
546+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
547+ for (i4 = 0; i4 < 4; i4++)
548+ for (j4 = 0; j4 < 4; j4++)
549+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
550+ libtetrabz_polcmplx3(ne, e0, de, w2);
551+ for (i4 = 0; i4 < 4; i4++)
552+ for (j4 = 0; j4 < 4; j4++)
553+ for (ie = 0; ie < ne; ie++)
554+ for (ir = 0; ir < 2; ir++)
555+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
556+ }
557+
558+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
559+
560+ if (v > thr) {
561+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
562+ for (i4 = 0; i4 < 4; i4++)
563+ for (j4 = 0; j4 < 4; j4++)
564+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
565+ libtetrabz_polcmplx3(ne, e0, de, w2);
566+ for (i4 = 0; i4 < 4; i4++)
567+ for (j4 = 0; j4 < 4; j4++)
568+ for (ie = 0; ie < ne; ie++)
569+ for (ir = 0; ir < 2; ir++)
570+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
571+ }
572+ }
573+ else if (e[3] <= 0.0) {
574+ for (i4 = 0; i4 < 4; i4++) {
575+ de[i4] = tsmall[i4][j4] * (ej1[i4][ib] - ei1[i4]);
576+ }
577+ libtetrabz_polcmplx3(ne, e0, de, w2);
578+ for (i4 = 0; i4 < 4; i4++)
579+ for (ie = 0; ie < ne; ie++)
580+ for (ir = 0; ir < 2; ir++)
581+ w1[ib][i4][ie][ir] += w2[i4][ie][ir];
582+ }
583+ }
584+ free(w2[0][0]);
585+ free(w2[0]);
586+ free(w2);
795587 }
796588
797589
798-void libtetrabz_polcmplx_1211(
799- double g1,
800- double g2,
801- double *w
802-) {
803- /*
804- 6, g4 = g3 = g1
805- */
806- /*
807- Real
808- */
809- w[0] = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1))
810- + g2 * (g2 * g2 - 3.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
811- w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
812- w[0] = -5.0 * g1 + w[0] / (g2 - g1);
813- w[0] = -11.0 + (3.0 * w[0]) / (g2 - g1);
814- w[0] = w[0] / (6.0 * (g2 - g1));
590+void polcmplx(
591+ int *ng,
592+ int nk,
593+ int nb,
594+ int ne,
595+ double **bvec,
596+ double **eig1,
597+ double **eig2,
598+ double **e0,
599+ double *****wght) {
815600 /*
816- Imaginal
601+ Main SUBROUTINE for Polarization (Imaginary axis) : Theta(- E1) * Theta(E2) / (E2 - E1 - iw)
817602 */
818- w[1] = 2.0 * g2 * (-3.0 + g2 * g2) * (atan(g2) - atan(g1))
819- + (1.0 - 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
820- w[1] = 4.0 * g2 + w[1] / (g2 - g1);
821- w[1] = 1.0 + w[1] / (g2 - g1);
822- w[1] = w[1] / (2.0 * (g2 - g1) * (g2 - g1));
603+ int it, ik, ie, ib, i4, j4, ir, jb, **ikv, indx[4], i20;
604+ double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], ***** w1, **** w2, v, tsmall[4][4], thr;
605+
606+ ikv = (int**)malloc(6 * nk * sizeof(int*));
607+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
608+ for (it = 0; it < 6 * nk; it++) {
609+ ikv[it] = ikv[0] + it * 20;
610+ }
611+
612+ ei1 = (double**)malloc(4 * sizeof(double*));
613+ ej1 = (double**)malloc(4 * sizeof(double*));
614+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
615+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
616+ for (i4 = 0; i4 < 4; i4++) {
617+ ei1[i4] = ei1[0] + i4 * nb;
618+ ej1[i4] = ej1[0] + i4 * nb;
619+ }
620+
621+ ej2 = (double**)malloc(nb * sizeof(double*));
622+ ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
623+ for (ib = 0; ib < nb; ib++)
624+ ej2[ib] = ej2[0] + ib * 4;
625+
626+ w1 = (double*****)malloc(nb * sizeof(double****));
627+ w1[0] = (double****)malloc(nb * 4 * sizeof(double***));
628+ w1[0][0] = (double***)malloc(nb * 4 * nb * sizeof(double**));
629+ w1[0][0][0] = (double**)malloc(nb * 4 * nb * ne * sizeof(double*));
630+ w1[0][0][0][0] = (double*)malloc(nb * 4 * nb * ne * 2 * sizeof(double));
631+ for (ib = 0; ib < nb; ib++) {
632+ w1[ib] = w1[0] + ib * 4;
633+ for (i4 = 0; i4 < 4; i4++) {
634+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
635+ for (jb = 0; jb < nb; jb++) {
636+ w1[ib][i4][jb] = w1[0][0][0] + ib * 4 * nb * ne + i4 * nb * ne + jb * ne;
637+ for (ie = 0; ie < ne; ie++) {
638+ w1[ib][i4][jb][ie] = w1[0][0][0][0] + ib * 4 * nb * ne * 2 + i4 * nb * ne * 2 + jb * ne * 2 + ie * 2;
639+ }
640+ }
641+ }
642+ }
643+
644+ w2 = (double****)malloc(nb * sizeof(double***));
645+ w2[0] = (double***)malloc(nb * 4 * sizeof(double**));
646+ w2[0][0] = (double**)malloc(nb * 4 * ne * sizeof(double*));
647+ w2[0][0][0] = (double*)malloc(nb * 4 * ne * 2 * sizeof(double));
648+ for (ib = 0; ib < nb; ib++) {
649+ w2[ib] = w2[0] + ib * 4;
650+ for (i4 = 0; i4 < 4; i4++) {
651+ w2[ib][i4] = w2[0][0] + ib * 4 * ne + i4 * ne;
652+ for (ie = 0; ie < ne; ie++) {
653+ w2[ib][i4][ie] = w2[0][0][0] + ib * 4 * ne * 2 + i4 * ne * 2 + ie * 2;
654+ }
655+ }
656+ }
657+
658+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
659+
660+ for (ik = 0; ik < nk; ik++)
661+ for (ib = 0; ib < nb; ib++)
662+ for (jb = 0; jb < nb; jb++)
663+ for (ie = 0; ie < ne; ie++)
664+ for (ir = 0; ir < 2; ir++)
665+ wght[ik][ib][jb][ie][ir] = 0.0;
666+
667+ thr = 1.0e-8;
668+
669+ for (it = 0; it < 6 * nk; it++) {
670+
671+ for (i4 = 0; i4 < 4; i4++)
672+ for (ib = 0; ib < nb; ib++) {
673+ ei1[i4][ib] = 0.0;
674+ ej1[i4][ib] = 0.0;
675+ }
676+ for (i20 = 0; i20 < 20; i20++) {
677+ for (i4 = 0; i4 < 4; i4++) {
678+ for (ib = 0; ib < nb; ib++) {
679+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
680+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
681+ }
682+ }
683+ }
684+
685+ for (ib = 0; ib < nb; ib++) {
686+
687+ for (jb = 0; jb < nb; jb++)
688+ for (i4 = 0; i4 < 4; i4++)
689+ for (ie = 0; ie < ne; ie++)
690+ for (ir = 0; ir < 2; ir++)
691+ w1[ib][jb][ie][i4][ir] = 0.0;
692+
693+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
694+ eig_sort(4, e, indx);
695+
696+ if (e[0] <= 0.0 && 0.0 < e[1]) {
697+
698+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
699+
700+ if (v > thr) {
701+ for (j4 = 0; j4 < 4; j4++) {
702+ ei2[j4] = 0.0;
703+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
704+ }
705+ for (i4 = 0; i4 < 4; i4++)
706+ for (j4 = 0; j4 < 4; j4++) {
707+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
708+ for (jb = 0; jb < nb; jb++)
709+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
710+ }
711+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
712+ for (i4 = 0; i4 < 4; i4++)
713+ for (jb = 0; jb < nb; jb++)
714+ for (j4 = 0; j4 < 4; j4++)
715+ for (ie = 0; ie < ne; ie++)
716+ for (ir = 0; ir < 2; ir++)
717+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
718+ }
719+ }
720+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
721+
722+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
723+
724+ if (v > thr) {
725+ for (j4 = 0; j4 < 4; j4++) {
726+ ei2[j4] = 0.0;
727+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
728+ }
729+ for (i4 = 0; i4 < 4; i4++)
730+ for (j4 = 0; j4 < 4; j4++) {
731+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
732+ for (jb = 0; jb < nb; jb++)
733+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
734+ }
735+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
736+ for (i4 = 0; i4 < 4; i4++)
737+ for (jb = 0; jb < nb; jb++)
738+ for (j4 = 0; j4 < 4; j4++)
739+ for (ie = 0; ie < ne; ie++)
740+ for (ir = 0; ir < 2; ir++)
741+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
742+ }
743+
744+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
745+
746+ if (v > thr) {
747+ for (j4 = 0; j4 < 4; j4++) {
748+ ei2[j4] = 0.0;
749+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
750+ }
751+ for (i4 = 0; i4 < 4; i4++)
752+ for (j4 = 0; j4 < 4; j4++) {
753+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
754+ for (jb = 0; jb < nb; jb++)
755+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
756+ }
757+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
758+ for (i4 = 0; i4 < 4; i4++)
759+ for (jb = 0; jb < nb; jb++)
760+ for (j4 = 0; j4 < 4; j4++)
761+ for (ie = 0; ie < ne; ie++)
762+ for (ir = 0; ir < 2; ir++)
763+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
764+ }
765+
766+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
767+
768+ if (v > thr) {
769+ for (j4 = 0; j4 < 4; j4++) {
770+ ei2[j4] = 0.0;
771+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
772+ }
773+ for (i4 = 0; i4 < 4; i4++)
774+ for (j4 = 0; j4 < 4; j4++) {
775+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
776+ for (jb = 0; jb < nb; jb++)
777+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
778+ }
779+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
780+ for (i4 = 0; i4 < 4; i4++)
781+ for (jb = 0; jb < nb; jb++)
782+ for (j4 = 0; j4 < 4; j4++)
783+ for (ie = 0; ie < ne; ie++)
784+ for (ir = 0; ir < 2; ir++)
785+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
786+ }
787+ }
788+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
789+
790+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
791+
792+ if (v > thr) {
793+ for (j4 = 0; j4 < 4; j4++) {
794+ ei2[j4] = 0.0;
795+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
796+ }
797+ for (i4 = 0; i4 < 4; i4++)
798+ for (j4 = 0; j4 < 4; j4++) {
799+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
800+ for (jb = 0; jb < nb; jb++)
801+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
802+ }
803+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
804+ for (i4 = 0; i4 < 4; i4++)
805+ for (jb = 0; jb < nb; jb++)
806+ for (j4 = 0; j4 < 4; j4++)
807+ for (ie = 0; ie < ne; ie++)
808+ for (ir = 0; ir < 2; ir++)
809+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
810+ }
811+
812+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
813+
814+ if (v > thr) {
815+ for (j4 = 0; j4 < 4; j4++) {
816+ ei2[j4] = 0.0;
817+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
818+ }
819+ for (i4 = 0; i4 < 4; i4++)
820+ for (j4 = 0; j4 < 4; j4++) {
821+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
822+ for (jb = 0; jb < nb; jb++)
823+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
824+ }
825+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
826+ for (i4 = 0; i4 < 4; i4++)
827+ for (jb = 0; jb < nb; jb++)
828+ for (j4 = 0; j4 < 4; j4++)
829+ for (ie = 0; ie < ne; ie++)
830+ for (ir = 0; ir < 2; ir++)
831+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
832+ }
833+
834+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
835+
836+ if (v > thr) {
837+ for (j4 = 0; j4 < 4; j4++) {
838+ ei2[j4] = 0.0;
839+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
840+ }
841+ for (i4 = 0; i4 < 4; i4++)
842+ for (j4 = 0; j4 < 4; j4++) {
843+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
844+ for (jb = 0; jb < nb; jb++)
845+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
846+ }
847+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
848+ for (i4 = 0; i4 < 4; i4++)
849+ for (jb = 0; jb < nb; jb++)
850+ for (j4 = 0; j4 < 4; j4++)
851+ for (ie = 0; ie < ne; ie++)
852+ for (ir = 0; ir < 2; ir++)
853+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
854+ }
855+ }
856+ else if (e[3] <= 0.0) {
857+
858+ for (i4 = 0; i4 < 4; i4++) {
859+ ei2[i4] = ej1[i4][ib];
860+ for (jb = 0; jb < nb; jb++)
861+ ej2[jb][i4] = ej1[i4][jb];
862+ }
863+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
864+ for (i4 = 0; i4 < 4; i4++)
865+ for (jb = 0; jb < nb; jb++)
866+ for (ie = 0; ie < ne; ie++)
867+ for (ir = 0; ir < 2; ir++)
868+ w1[ib][i4][jb][ie][ir] += w2[jb][i4][ie][ir];
869+ }
870+ else {
871+ continue;
872+ }
873+ }
874+ for (i4 = 0; i4 < 20; i4++)
875+ for (ib = 0; ib < nb; ib++)
876+ for (j4 = 0; j4 < 4; j4++)
877+ for (jb = 0; jb < nb; jb++)
878+ for (ie = 0; ie < ne; ie++)
879+ for (ir = 0; ir < 2; ir++)
880+ wght[ikv[it][i20]][ib][jb][ie][ir] += wlsm[i20][j4] * w1[ib][j4][jb][ie][ir];
881+ }
882+ for (ik = 0; ik < nk; ik++)
883+ for (ib = 0; ib < nb; ib++)
884+ for (jb = 0; jb < nb; jb++)
885+ for (ie = 0; ie < ne; ie++)
886+ for (i4 = 0; i4 < 2; i4++)
887+ wght[ik][ib][jb][ie][i4] /= (6.0 * (double) nk);
888+
889+ free(ikv[0]);
890+ free(ikv);
891+ free(ei1[0]);
892+ free(ei1);
893+ free(ej1[0]);
894+ free(ej1);
895+ free(ej2[0]);
896+ free(ej2);
897+ free(w1[0][0][0][0]);
898+ free(w1[0][0][0]);
899+ free(w1[0][0]);
900+ free(w1[0]);
901+ free(w1);
902+ free(w2[0][0][0]);
903+ free(w2[0][0]);
904+ free(w2[0]);
905+ free(w2);
823906 }
--- a/src/c/libtetrabz_polstat.c
+++ b/src/c/libtetrabz_polstat.c
@@ -23,6 +23,7 @@
2323 #include "libtetrabz_common.h"
2424 #include <math.h>
2525 #include <stdio.h>
26+#include <stdlib.h>
2627
2728 /*
2829 Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
@@ -155,25 +156,25 @@ void libtetrabz_polstat3(
155156 /*
156157 Tetrahedron method for delta(om - ep + e)
157158 */
158- int ii, indx[4];
159+ int i4, indx[4];
159160 double thr, thr2, e[4], ln[4];
160161
161- for (ii = 0; ii < 3; ii++) e[ii] = de[ii];
162+ for (i4 = 0; i4 < 4; i4++) e[i4] = de[i4];
162163 eig_sort(4, e, indx);
163164
164165 thr = e[3] * 1.0e-3;
165166 thr2 = 1.0e-8;
166167
167- for(ii=0;ii<4;ii++){
168- if (e[ii] < thr2) {
169- if (ii == 3) {
168+ for(i4 =0; i4 <4; i4++){
169+ if (e[i4] < thr2) {
170+ if (i4 == 3) {
170171 printf(" Nesting # ");
171172 }
172- ln[ii] = 0.0;
173- e[ii] = 0.0;
173+ ln[i4] = 0.0;
174+ e[i4] = 0.0;
174175 }
175176 else{
176- ln[ii] = log(e[ii]);
177+ ln[i4] = log(e[i4]);
177178 }
178179 }
179180
@@ -233,55 +234,55 @@ void libtetrabz_polstat3(
233234 printf("weighting 4=3");
234235 }
235236 }
236- }
237- else if(fabs(e[2] - e[1]) < thr){
238- if (fabs(e[2] - e[0]) < thr) {
239- /*
240- e[2] = e[1] = e[0]
241- */
242- w1[indx[3]] = libtetrabz_polstat_1222(e[3], e[2], ln[3], ln[2]);
243- w1[indx[2]] = libtetrabz_polstat_1211(e[2], e[3], ln[2], ln[3]);
244- w1[indx[1]] = w1[indx[2]];
245- w1[indx[0]] = w1[indx[2]];
246-
247- if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
248- printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
249- printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
250- printf("weighting 3=2=1");
251- }
252- }
253- else {
254- /*
255- e[2] = e[1]
256- */
257- w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[0], e[2], ln[3], ln[0], ln[2]);
258- w1[indx[2]] = libtetrabz_polstat_1231(e[2], e[0], e[3], ln[2], ln[0], ln[3]);
259- w1[indx[1]] = w1[indx[2]];
260- w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[3], e[2], ln[0], ln[3], ln[2]);
237+ }
238+ else if (fabs(e[2] - e[1]) < thr) {
239+ if (fabs(e[2] - e[0]) < thr) {
240+ /*
241+ e[2] = e[1] = e[0]
242+ */
243+ w1[indx[3]] = libtetrabz_polstat_1222(e[3], e[2], ln[3], ln[2]);
244+ w1[indx[2]] = libtetrabz_polstat_1211(e[2], e[3], ln[2], ln[3]);
245+ w1[indx[1]] = w1[indx[2]];
246+ w1[indx[0]] = w1[indx[2]];
261247
262- if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
263- printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
264- printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
265- printf("weighting 3=2");
266- }
248+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
249+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
250+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
251+ printf("weighting 3=2=1");
267252 }
268253 }
269- else if(fabs(e[1] - e[0]) < thr) {
254+ else {
270255 /*
271- e[1] = e[0]
256+ e[2] = e[1]
272257 */
273- w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[2], e[1], ln[3], ln[2], ln[1]);
274- w1[indx[2]] = libtetrabz_polstat_1233(e[2], e[3], e[1], ln[2], ln[3], ln[1]);
275- w1[indx[1]] = libtetrabz_polstat_1231(e[1], e[2], e[3], ln[1], ln[2], ln[3]);
276- w1[indx[0]] = w1[indx[1]];
258+ w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[0], e[2], ln[3], ln[0], ln[2]);
259+ w1[indx[2]] = libtetrabz_polstat_1231(e[2], e[0], e[3], ln[2], ln[0], ln[3]);
260+ w1[indx[1]] = w1[indx[2]];
261+ w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[3], e[2], ln[0], ln[3], ln[2]);
277262
278263 if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
279264 printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
280265 printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
281- printf("weighting 2=1");
266+ printf("weighting 3=2");
282267 }
283268 }
284- else {
269+ }
270+ else if (fabs(e[1] - e[0]) < thr) {
271+ /*
272+ e[1] = e[0]
273+ */
274+ w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[2], e[1], ln[3], ln[2], ln[1]);
275+ w1[indx[2]] = libtetrabz_polstat_1233(e[2], e[3], e[1], ln[2], ln[3], ln[1]);
276+ w1[indx[1]] = libtetrabz_polstat_1231(e[1], e[2], e[3], ln[1], ln[2], ln[3]);
277+ w1[indx[0]] = w1[indx[1]];
278+
279+ if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
280+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
281+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
282+ printf("weighting 2=1");
283+ }
284+ }
285+ else {
285286 /*
286287 Different each other.
287288 */
@@ -310,14 +311,14 @@ void libtetrabz_polstat2(
310311 :param ej1:
311312 :return:
312313 */
313- int ii, jj, ib, indx[4];
314+ int i4, j4, ib, indx[4];
314315 double v, thr, e[4], de[4], w2[4], tsmall[4][4];
315316
316317 thr = 1.0e-8;
317318
318319 for (ib = 0; ib < nb; ib++) {
319320
320- for (ii = 0; ii < 4; ii++) e[ii] = -ej1[ii][ib];
321+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
321322 eig_sort(4, e, indx);
322323
323324 if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
@@ -325,15 +326,14 @@ void libtetrabz_polstat2(
325326 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
326327
327328 if (v > thr) {
328- for (ii = 0; ii < 4; ii++) {
329- de[ii] = 0.0;
330- for (jj = 0; jj < 4; jj++)
331- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
332- }
329+ for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
330+ for (i4 = 0; i4 < 4; i4++)
331+ for (j4 = 0; j4 < 4; j4++)
332+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
333333 libtetrabz_polstat3(de, w2);
334- for (ii = 0; ii < 4; ii++)
335- for (jj = 0; jj < 4; jj++)
336- w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
334+ for (i4 = 0; i4 < 4; i4++)
335+ for (j4 = 0; j4 < 4; j4++)
336+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
337337 }
338338 }
339339 else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
@@ -341,43 +341,40 @@ void libtetrabz_polstat2(
341341 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
342342
343343 if (v > thr) {
344- for (ii = 0; ii < 4; ii++) {
345- de[ii] = 0.0;
346- for (jj = 0; jj < 4; jj++)
347- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
348- }
344+ for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
345+ for (i4 = 0; i4 < 4; i4++)
346+ for (j4 = 0; j4 < 4; j4++)
347+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
349348 libtetrabz_polstat3(de, w2);
350- for (ii = 0; ii < 4; ii++)
351- for (jj = 0; jj < 4; jj++)
352- w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
349+ for (i4 = 0; i4 < 4; i4++)
350+ for (j4 = 0; j4 < 4; j4++)
351+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
353352 }
354353
355354 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
356355
357356 if (v > thr) {
358- for (ii = 0; ii < 4; ii++) {
359- de[ii] = 0.0;
360- for (jj = 0; jj < 4; jj++)
361- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
362- }
357+ for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
358+ for (i4 = 0; i4 < 4; i4++)
359+ for (j4 = 0; j4 < 4; j4++)
360+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
363361 libtetrabz_polstat3(de, w2);
364- for (ii = 0; ii < 4; ii++)
365- for (jj = 0; jj < 4; jj++)
366- w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
362+ for (i4 = 0; i4 < 4; i4++)
363+ for (j4 = 0; j4 < 4; j4++)
364+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
367365 }
368366
369367 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
370368
371369 if (v > thr) {
372- for (ii = 0; ii < 4; ii++) {
373- de[ii] = 0.0;
374- for (jj = 0; jj < 4; jj++)
375- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
376- }
370+ for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
371+ for (i4 = 0; i4 < 4; i4++)
372+ for (j4 = 0; j4 < 4; j4++)
373+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
377374 libtetrabz_polstat3(de, w2);
378- for (ii = 0; ii < 4; ii++)
379- for (jj = 0; jj < 4; jj++)
380- w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
375+ for (i4 = 0; i4 < 4; i4++)
376+ for (j4 = 0; j4 < 4; j4++)
377+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
381378 }
382379 }
383380 else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
@@ -385,50 +382,47 @@ void libtetrabz_polstat2(
385382 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
386383
387384 if (v > thr) {
388- for (ii = 0; ii < 4; ii++) {
389- de[ii] = 0.0;
390- for (jj = 0; jj < 4; jj++)
391- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
392- }
385+ for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
386+ for (i4 = 0; i4 < 4; i4++)
387+ for (j4 = 0; j4 < 4; j4++)
388+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
393389 libtetrabz_polstat3(de, w2);
394- for (ii = 0; ii < 4; ii++)
395- for (jj = 0; jj < 4; jj++)
396- w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
390+ for (i4 = 0; i4 < 4; i4++)
391+ for (j4 = 0; j4 < 4; j4++)
392+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
397393 }
398394
399395 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
400396
401397 if (v > thr) {
402- for (ii = 0; ii < 4; ii++) {
403- de[ii] = 0.0;
404- for (jj = 0; jj < 4; jj++)
405- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
406- }
398+ for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
399+ for (i4 = 0; i4 < 4; i4++)
400+ for (j4 = 0; j4 < 4; j4++)
401+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
407402 libtetrabz_polstat3(de, w2);
408- for (ii = 0; ii < 4; ii++)
409- for (jj = 0; jj < 4; jj++)
410- w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
403+ for (i4 = 0; i4 < 4; i4++)
404+ for (j4 = 0; j4 < 4; j4++)
405+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
411406 }
412407
413408 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
414409
415410 if (v > thr) {
416- for (ii = 0; ii < 4; ii++) {
417- de[ii] = 0.0;
418- for (jj = 0; jj < 4; jj++)
419- de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
420- }
411+ for (j4 = 0; j4 < 4; j4++)de[j4] = 0.0;
412+ for (i4 = 0; i4 < 4; i4++)
413+ for (j4 = 0; j4 < 4; j4++)
414+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
421415 libtetrabz_polstat3(de, w2);
422- for (ii = 0; ii < 4; ii++)
423- for (jj = 0; jj < 4; jj++)
424- w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
416+ for (i4 = 0; i4 < 4; i4++)
417+ for (j4 = 0; j4 < 4; j4++)
418+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
425419 }
426420 }
427421 else if (e[3] <= 0.0) {
428- for (ii = 0; ii < 4; ii++)
429- de[ii] = ej1[indx[ii]][ib] - ei1[indx[ii]];
422+ for (i4 = 0; i4 < 4; i4++)
423+ de[i4] = ej1[i4][ib] - ei1[i4];
430424 libtetrabz_polstat3(de, w2);
431- for (ii = 0; ii < 4; ii++) w1[ib][ii] += w2[ii];
425+ for (i4 = 0; i4 < 4; i4++) w1[ib][i4] += w2[i4];
432426 }
433427 }
434428 }
@@ -446,11 +440,48 @@ void polstat(
446440 /*
447441 Compute Static polarization function
448442 */
449- int it, ik, ib, ii, jj, jb, **ikv, indx[4];
450- double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], ***w1, **w2, v, tsmall[4][4], thr;
443+ int it, ik, ib, jb, i20, i4, j4, **ikv, indx[4];
444+ double wlsm[20][4], **ei1, **ej1, ei2[4], ** ej2, e[4], *** w1, ** w2, v, tsmall[4][4], thr;
451445
452446 thr = 1.0e-10;
453447
448+ ikv = (int**)malloc(6 * nk * sizeof(int*));
449+ ikv[0] = (int*)malloc(6 * nk * 20 * sizeof(int));
450+ for (ik = 0; ik < 6 * nk; ik++) {
451+ ikv[ik] = ikv[0] + ik * 20;
452+ }
453+
454+ ei1 = (double**)malloc(4 * sizeof(double*));
455+ ej1 = (double**)malloc(4 * sizeof(double*));
456+ ei1[0] = (double*)malloc(4 * nb * sizeof(double));
457+ ej1[0] = (double*)malloc(4 * nb * sizeof(double));
458+ for (i4 = 0; i4 < 4; i4++) {
459+ ei1[i4] = ei1[0] + i4 * nb;
460+ ej1[i4] = ej1[0] + i4 * nb;
461+ }
462+
463+ ej2 = (double**)malloc(nb * sizeof(double*));
464+ ej2[0] = (double*)malloc(nb * 4 * sizeof(double));
465+ for (ib = 0; ib < nb; ib++) {
466+ ej2[ib] = ej2[0] + ib * 4;
467+ }
468+
469+ w1 = (double***)malloc(nb * sizeof(double**));
470+ w1[0] = (double**)malloc(nb * 4 * sizeof(double*));
471+ w1[0][0] = (double*)malloc(nb * 4 * nb * sizeof(double));
472+ for (ib = 0; ib < nb; ib++) {
473+ w1[ib] = w1[0] + ib * 4;
474+ for (i4 = 0; i4 < 4; i4++) {
475+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
476+ }
477+ }
478+
479+ w2 = (double**)malloc(nb * sizeof(double*));
480+ w2[0] = (double*)malloc(nb * 4 * sizeof(double));
481+ for (ib = 0; ib < nb; ib++) {
482+ w2[ib] = w2[0] + ib * 4;
483+ }
484+
454485 libtetrabz_initialize(ng, bvec, wlsm, ikv);
455486
456487 for (ik = 0; ik < nk; ik++)
@@ -460,25 +491,27 @@ void polstat(
460491
461492 for(it = 0; it < 6*nk; it++) {
462493
463- for (ii = 0; ii < 4; ii++) {
494+ for (i4 = 0; i4 < 4; i4++)
464495 for (ib = 0; ib < nb; ib++) {
465- ei1[ii][ib] = 0.0;
466- ej1[ii][ib] = 0.0;
467- for (jj = 0; jj < 20; jj++) {
468- ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
469- ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
496+ ei1[i4][ib] = 0.0;
497+ ej1[i4][ib] = 0.0;
498+ }
499+ for (i20 = 0; i20 < 20; i20++) {
500+ for (i4 = 0; i4 < 4; i4++) {
501+ for (ib = 0; ib < nb; ib++) {
502+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
503+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
470504 }
471505 }
472506 }
473507
474- for (ib = 0; ib < nb; ib++)
475- for (jb = 0; jb < nb; jb++)
476- for (ii = 0; ii < 4; ii++)
477- w1[ib][jb][ii] = 0.0;
478-
479508 for (ib = 0; ib < nb; ib++) {
480509
481- for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
510+ for (i4 = 0; i4 < 4; i4++)
511+ for (jb = 0; jb < nb; jb++)
512+ w1[ib][i4][jb] = 0.0;
513+
514+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
482515 eig_sort(4, e, indx);
483516
484517 if (e[0] <= 0.0 && 0.0 < e[1]) {
@@ -486,17 +519,21 @@ void polstat(
486519 libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
487520
488521 if (v > thr) {
489- for (ii = 0; ii < 4; ii++)
490- for (jj = 0; jj < 4; jj++) {
491- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
522+ for (j4 = 0; j4 < 4; j4++) {
523+ ei2[j4] = 0.0;
524+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
525+ }
526+ for (i4 = 0; i4 < 4; i4++)
527+ for (j4 = 0; j4 < 4; j4++) {
528+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
492529 for (jb = 0; jb < nb; jb++)
493- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
530+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
494531 }
495532 libtetrabz_polstat2(nb, ei2, ej2, w2);
496- for (ii = 0; ii < 4; ii++)
497- for (jj = 0; jj < 4; jj++)
498- for (jb = 0; jb < nb; jb++)
499- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
533+ for (i4 = 0; i4 < 4; i4++)
534+ for (jb = 0; jb < nb; jb++)
535+ for (j4 = 0; j4 < 4; j4++)
536+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
500537 }
501538 }
502539 else if (e[1] <= 0.0 && 0.0 < e[2]) {
@@ -504,49 +541,61 @@ void polstat(
504541 libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
505542
506543 if (v > thr) {
507- for (ii = 0; ii < 4; ii++)
508- for (jj = 0; jj < 4; jj++) {
509- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
544+ for (j4 = 0; j4 < 4; j4++) {
545+ ei2[j4] = 0.0;
546+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
547+ }
548+ for (i4 = 0; i4 < 4; i4++)
549+ for (j4 = 0; j4 < 4; j4++) {
550+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
510551 for (jb = 0; jb < nb; jb++)
511- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
552+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
512553 }
513554 libtetrabz_polstat2(nb, ei2, ej2, w2);
514- for (ii = 0; ii < 4; ii++)
515- for (jj = 0; jj < 4; jj++)
516- for (jb = 0; jb < nb; jb++)
517- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
555+ for (i4 = 0; i4 < 4; i4++)
556+ for (jb = 0; jb < nb; jb++)
557+ for (j4 = 0; j4 < 4; j4++)
558+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
518559 }
519560
520561 libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
521562
522563 if (v > thr) {
523- for (ii = 0; ii < 4; ii++)
524- for (jj = 0; jj < 4; jj++) {
525- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
564+ for (j4 = 0; j4 < 4; j4++) {
565+ ei2[j4] = 0.0;
566+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
567+ }
568+ for (i4 = 0; i4 < 4; i4++)
569+ for (j4 = 0; j4 < 4; j4++) {
570+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
526571 for (jb = 0; jb < nb; jb++)
527- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
572+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
528573 }
529574 libtetrabz_polstat2(nb, ei2, ej2, w2);
530- for (ii = 0; ii < 4; ii++)
531- for (jj = 0; jj < 4; jj++)
532- for (jb = 0; jb < nb; jb++)
533- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
575+ for (i4 = 0; i4 < 4; i4++)
576+ for (jb = 0; jb < nb; jb++)
577+ for (j4 = 0; j4 < 4; j4++)
578+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
534579 }
535580
536581 libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
537582
538583 if (v > thr) {
539- for (ii = 0; ii < 4; ii++)
540- for (jj = 0; jj < 4; jj++) {
541- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
584+ for (j4 = 0; j4 < 4; j4++) {
585+ ei2[j4] = 0.0;
586+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
587+ }
588+ for (i4 = 0; i4 < 4; i4++)
589+ for (j4 = 0; j4 < 4; j4++) {
590+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
542591 for (jb = 0; jb < nb; jb++)
543- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
592+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
544593 }
545594 libtetrabz_polstat2(nb, ei2, ej2, w2);
546- for (ii = 0; ii < 4; ii++)
547- for (jj = 0; jj < 4; jj++)
548- for (jb = 0; jb < nb; jb++)
549- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
595+ for (i4 = 0; i4 < 4; i4++)
596+ for (jb = 0; jb < nb; jb++)
597+ for (j4 = 0; j4 < 4; j4++)
598+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
550599 }
551600 }
552601 else if (e[2] <= 0.0 && 0.0 < e[3]) {
@@ -554,74 +603,100 @@ void polstat(
554603 libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
555604
556605 if (v > thr) {
557- for (ii = 0; ii < 4; ii++)
558- for (jj = 0; jj < 4; jj++) {
559- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
606+ for (j4 = 0; j4 < 4; j4++) {
607+ ei2[j4] = 0.0;
608+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
609+ }
610+ for (i4 = 0; i4 < 4; i4++)
611+ for (j4 = 0; j4 < 4; j4++) {
612+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
560613 for (jb = 0; jb < nb; jb++)
561- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
614+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
562615 }
563616 libtetrabz_polstat2(nb, ei2, ej2, w2);
564- for (ii = 0; ii < 4; ii++)
565- for (jj = 0; jj < 4; jj++)
566- for (jb = 0; jb < nb; jb++)
567- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
617+ for (i4 = 0; i4 < 4; i4++)
618+ for (jb = 0; jb < nb; jb++)
619+ for (j4 = 0; j4 < 4; j4++)
620+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
568621 }
569622
570623 libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
571624
572625 if (v > thr) {
573- for (ii = 0; ii < 4; ii++)
574- for (jj = 0; jj < 4; jj++) {
575- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
626+ for (j4 = 0; j4 < 4; j4++) {
627+ ei2[j4] = 0.0;
628+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
629+ }
630+ for (i4 = 0; i4 < 4; i4++)
631+ for (j4 = 0; j4 < 4; j4++) {
632+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
576633 for (jb = 0; jb < nb; jb++)
577- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
634+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
578635 }
579636 libtetrabz_polstat2(nb, ei2, ej2, w2);
580- for (ii = 0; ii < 4; ii++)
581- for (jj = 0; jj < 4; jj++)
582- for (jb = 0; jb < nb; jb++)
583- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
637+ for (i4 = 0; i4 < 4; i4++)
638+ for (jb = 0; jb < nb; jb++)
639+ for (j4 = 0; j4 < 4; j4++)
640+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
584641 }
585642
586643 libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
587644
588645 if (v > thr) {
589- for (ii = 0; ii < 4; ii++)
590- for (jj = 0; jj < 4; jj++) {
591- ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
646+ for (j4 = 0; j4 < 4; j4++) {
647+ ei2[j4] = 0.0;
648+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
649+ }
650+ for (i4 = 0; i4 < 4; i4++)
651+ for (j4 = 0; j4 < 4; j4++) {
652+ ei2[j4] += ej1[indx[i4]][ib] * tsmall[i4][j4];
592653 for (jb = 0; jb < nb; jb++)
593- ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
654+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
594655 }
595656 libtetrabz_polstat2(nb, ei2, ej2, w2);
596- for (ii = 0; ii < 4; ii++)
597- for (jj = 0; jj < 4; jj++)
598- for (jb = 0; jb < nb; jb++)
599- w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
657+ for (i4 = 0; i4 < 4; i4++)
658+ for (jb = 0; jb < nb; jb++)
659+ for (j4 = 0; j4 < 4; j4++)
660+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
600661 }
601- }
662+ }
602663 else if (e[3] <= 0.0) {
603- for (ii = 0; ii < 4; ii++) {
604- ei2[ii] = ej1[ii][ib];
664+ for (i4 = 0; i4 < 4; i4++) {
665+ ei2[i4] = ej1[i4][ib];
605666 for (jb = 0; jb < nb; jb++)
606- ej2[ii][jb] = ej1[ii][jb];
667+ ej2[jb][i4] = ej1[i4][jb];
607668 }
608669 libtetrabz_polstat2(nb, ei2, ej2, w2);
609- for (ii = 0; ii < 4; ii++)
670+ for (i4 = 0; i4 < 4; i4++)
610671 for (jb = 0; jb < nb; jb++)
611- w1[ib][jb][ii] += v * w2[jb][ii];
672+ w1[ib][i4][jb] += w2[jb][i4];
612673 }
613674 else {
614675 continue;
615676 }
616677 }
617- for (ii = 0; ii < 20; ii++)
618- for (jj = 0; jj < 4; jj++)
619- for (ib = 0; ib < nb; ib++)
678+ for (i20 = 0; i20 < 20; i20++)
679+ for (ib = 0; ib < nb; ib++)
680+ for (i4 = 0; i4 < 4; i4++)
620681 for (jb = 0; jb < nb; jb++)
621- wght[ikv[it][ii]][ib][jb] += w1[ib][jb][jj] * wlsm[jj][ii];
682+ wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
622683 }
623684 for (ik = 0; ik < nk; ik++)
624685 for (ib = 0; ib < nb; ib++)
625686 for (jb = 0; jb < nb; jb++)
626687 wght[ik][ib][jb] /= (6.0 * (double) nk);
688+
689+ free(ikv[0]);
690+ free(ikv);
691+ free(ei1[0]);
692+ free(ei1);
693+ free(ej1[0]);
694+ free(ej1);
695+ free(ej2[0]);
696+ free(ej2);
697+ free(w1[0][0]);
698+ free(w1[0]);
699+ free(w1);
700+ free(w2[0]);
701+ free(w2);
627702 }
Show on old repository browser