OSDN Git Service

Backup
authorMitsuaki Kawamura <kawamitsuaki@gmail.com>
Fri, 18 Mar 2022 02:56:36 +0000 (11:56 +0900)
committerMitsuaki Kawamura <kawamitsuaki@gmail.com>
Fri, 18 Mar 2022 02:56:36 +0000 (11:56 +0900)
src/c/libtetrabz_dbldelta.c

index ef53fb9..3924731 100644 (file)
@@ -23,6 +23,7 @@
 #include "libtetrabz_common.h"
 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
 
 void libtetrabz_dbldelta2(
     int nb,
@@ -32,15 +33,19 @@ void libtetrabz_dbldelta2(
   /*
   2nd step of tetrahedron method.
   */
-  int ii, ib, indx[4];
-  double a10, a20, a02, a12, v, e[3];
+  int i3, ib, indx[3];
+  double a10, a20, a02, a12, v, e[3], e_abs;
 
   for (ib = 0; ib < nb; ib++) {
 
-    for (ii = 0; ii < 3; ii++) e[ii] = ej[ib][ii];
+    e_abs = 0.0;
+    for (i3 = 0; i3 < 3; i3++) {
+      e[i3] = ej[ib][i3];
+      if (e_abs < fabs(e[i3])) e_abs = fabs(e[i3]);
+    }
     eig_sort(3, e, indx);
 
-    if (e[2] < 1.0e-10) {
+    if (e_abs < 1.0e-10) {
       printf("Nesting ##\n");
     }
 
@@ -67,8 +72,8 @@ void libtetrabz_dbldelta2(
       w[ib][indx[2]] = v * (2.0 - a02 - a12);
     }
     else {
-      for (ii = 0; ii < 3; ii++)
-        w[ib][ii] = 0.0;
+      for (i3 = 0; i3 < 3; i3++)
+        w[ib][i3] = 0.0;
     }
   }
 }
@@ -123,7 +128,7 @@ void dbldelta(
   }
 
   w2 = (double**)malloc(nb * sizeof(double*));
-  w2[0] = (double*)malloc(nb*3 * sizeof(double));
+  w2[0] = (double*)malloc(nb * 3 * sizeof(double));
   for (ib = 0; ib < nb; ib++) {
     w2[ib] = w2[0] + ib * 3;
   }
@@ -160,16 +165,16 @@ void dbldelta(
       for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
       eig_sort(4, e, indx);
 
-      if (e[0] < 0.0 && e[0] <= e[1]) {
+      if (e[0] < 0.0 && 0.0 <= e[1]) {
 
         libtetrabz_triangle_a1(e, 0.0, &v, tsmall);
 
         if (v > thr) {
-          for (j3 = 0; j3 < 3; j3++)
-            for (jb = 0; jb < nb; jb++) ej2[jb][j3] = 0.0;
+          for (jb = 0; jb < nb; jb++)
+            for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
           for (i4 = 0; i4 < 4; i4++)
-            for (j3 = 0; j3 < 3; j3++)
-              for (jb = 0; jb < nb; jb++)
+            for (jb = 0; jb < nb; jb++)
+              for (j3 = 0; j3 < 3; j3++)
                 ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
           libtetrabz_dbldelta2(nb, ej2, w2);
           for (i4 = 0; i4 < 4; i4++)
@@ -183,11 +188,11 @@ void dbldelta(
         libtetrabz_triangle_b1(e, 0.0, &v, tsmall);
 
         if (v > thr) {
-          for (j3 = 0; j3 < 3; j3++)
-            for (jb = 0; jb < nb; jb++) ej2[jb][j3] = 0.0;
+          for (jb = 0; jb < nb; jb++)
+            for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
           for (i4 = 0; i4 < 4; i4++)
-            for (j3 = 0; j3 < 3; j3++)
-              for (jb = 0; jb < nb; jb++)
+            for (jb = 0; jb < nb; jb++)
+              for (j3 = 0; j3 < 3; j3++)
                 ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
           libtetrabz_dbldelta2(nb, ej2, w2);
           for (i4 = 0; i4 < 4; i4++)
@@ -199,11 +204,11 @@ void dbldelta(
         libtetrabz_triangle_b2(e, 0.0, &v, tsmall);
 
         if (v > thr) {
-          for (j3 = 0; j3 < 3; j3++)
-            for (jb = 0; jb < nb; jb++) ej2[jb][j3] = 0.0;
+          for (jb = 0; jb < nb; jb++)
+            for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
           for (i4 = 0; i4 < 4; i4++)
-            for (j3 = 0; j3 < 3; j3++)
-              for (jb = 0; jb < nb; jb++)
+            for (jb = 0; jb < nb; jb++)
+              for (j3 = 0; j3 < 3; j3++)
                 ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
           libtetrabz_dbldelta2(nb, ej2, w2);
           for (i4 = 0; i4 < 4; i4++)
@@ -217,11 +222,11 @@ void dbldelta(
         libtetrabz_triangle_c1(e, 0.0, &v, tsmall);
 
         if (v > thr) {
-          for (j3 = 0; j3 < 3; j3++)
-            for (jb = 0; jb < nb; jb++) ej2[jb][j3] = 0.0;
+          for (jb = 0; jb < nb; jb++)
+            for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
           for (i4 = 0; i4 < 4; i4++)
-            for (j3 = 0; j3 < 3; j3++)
-              for (jb = 0; jb < nb; jb++)
+            for (jb = 0; jb < nb; jb++)
+              for (j3 = 0; j3 < 3; j3++)
                 ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
           libtetrabz_dbldelta2(nb, ej2, w2);
           for (i4 = 0; i4 < 4; i4++)