Changeset 42389 for trunk/Ohana/src/opihi/lib.data/svdcmp.c

Timestamp:

Feb 8, 2023, 12:17:35 PM (3 years ago)

Author:

eugene

Message:

merge from eam_branches/ipp-20220316. fixes for more pendantic gcc; add opihi memory stats; string vector improvements; use named macros for fixed string lengths

Location:

trunk/Ohana

Files:

• . (modified) (1 prop)
• src/opihi/lib.data/svdcmp.c (modified) (10 diffs)

trunk/Ohana/src/opihi/lib.data/svdcmp.c

@@ -23 +23 @@
   l = nm = 0;
   ALLOCATE (rv1, float, Nx);
+
   for (i = 0; i < Nx; i++) {
     l = i + 1;
@@ -29 +30 @@
     if (i < Ny) {
       for (k = i; k < Ny; k++) scale += fabs(a[k*Nx + i]);
-      if (scale) {
+      if (scale > 0.0) {
         for (k = i; k < Ny; k++) {
           a[k*Nx + i] /= scale;
@@ -35 +36 @@
         }
         f = a[i*Nx + i];
-        g  = -FSIGN(sqrt(s), f);
+        g = sqrt(s);
+        if (f > 0.0) { g = -g; }
         h = f*g - s;
         a[i*Nx + i] = f-g;
         if (i != Nx - 1) {
           for (j = l; j < Nx; j++) {
-            for (s = 0.0, k = i; k < Ny; k++) s += a[k*Nx + i]*a[k*Nx + j];
+            s = 0.0;
+            for (k = i; k < Ny; k++) { s += a[k*Nx + i]*a[k*Nx + j]; }
             f = s/h;
-            for (k = i; k < Ny; k++) a[k*Nx + j] += f*a[k*Nx + i];
+            for (k = i; k < Ny; k++) { a[k*Nx + j] += f*a[k*Nx + i]; }
           }
         }
-        for (k = i; k < Ny; k++) a[k*Nx + i] *= scale;
+        for (k = i; k < Ny; k++) { a[k*Nx + i] *= scale; }
       }
     }
@@ -51 +54 @@
     g = s = scale = 0.0;
     if ((i < Ny) && (i != (Nx - 1))) {
-      for (k = l; k < Nx; k++) scale += fabs(a[i*Nx + k]);
-      if (scale) {
+      for (k = l; k < Nx; k++) { scale += fabs(a[i*Nx + k]); }
+      if (scale > 0.0) {
         for (k = l; k < Nx; k++) {
           a[i*Nx + k] /= scale;
@@ -58 +61 @@
         }
         f = a[i*Nx + l];
-        g = -FSIGN(sqrt(s), f);
+        g = sqrt(s);
+        if (f > 0.0) { g = -g; }
         h = f*g - s;
         a[i*Nx + l] = f-g;
-        for (k = l; k < Nx; k++) rv1[k] = a[i*Nx + k]/h;
+        for (k = l; k < Nx; k++) { rv1[k] = a[i*Nx + k]/h; }
         if (i != Ny - 1) {
           for (j = l; j < Ny; j++) {
-            for (s = 0.0, k = l; k < Nx; k++) s += a[j*Nx + k]*a[i*Nx + k];
-            for (k = l; k < Nx; k++) a[j*Nx + k] += s*rv1[k];
+            s = 0.0;
+            for (k = l; k < Nx; k++) { s += a[j*Nx + k]*a[i*Nx + k];}
+            for (k = l; k < Nx; k++) { a[j*Nx + k] += s*rv1[k]; } 
           }
         }
-        for (k = l; k < Nx; k++) a[i*Nx + k] *= scale;
+        for (k = l; k < Nx; k++) { a[i*Nx + k] *= scale; }
       }
     }
@@ -76 +81 @@
   for (i = Nx - 1; i >= 0; i--) {
     if (i < Nx - 1) {
-      if (g) {
+      if (g != 0.0) {
         /* isn't l == n to start?? */
         for (j = l; j < Nx; j++) v[j*Nx + i] = (a[i*Nx + j]/a[i*Nx + l])/g;
         for (j = l; j < Nx; j++) {
-          for (s = 0.0, k = l; k < Nx; k++) s += a[i*Nx + k]*v[k*Nx + j];
+          s = 0.0;
+          for (k = l; k < Nx; k++) s += a[i*Nx + k]*v[k*Nx + j];
           for (k = l; k < Nx; k++) v[k*Nx + j] += s*v[k*Nx + i];
         }
       }
-      for (j = l; j < Nx; j++) v[i*Nx + j] = v[j*Nx + i] = 0.0;
+      for (j = l; j < Nx; j++) { v[i*Nx + j] = v[j*Nx + i] = 0.0; }
     }
     v[i*Nx + i] = 1.0;
@@ -93 +99 @@
     l = i + 1;
     g = w[i];
-    if (i < Nx - 1)
-      for (j = l; j < Nx; j++) a[i*Nx + j] = 0.0;
-    if (g) {
+    if (i < Nx - 1) { for (j = l; j < Nx; j++) {a[i*Nx + j] = 0.0; }}
+    if (g != 0.0) {
       g = 1.0/g;
       if (i != Nx - 1) {
         for (j = l; j < Nx; j++) {
-          for (s = 0.0, k = l; k < Ny; k++) s += a[k*Nx + i]*a[k*Nx + j];
+          s = 0.0;
+          for (k = l; k < Ny; k++) { s += a[k*Nx + i]*a[k*Nx + j];}
           f = (s/a[i*Nx + i])*g;
-          for (k = i; k < Ny; k++) a[k*Nx + j] += f*a[k*Nx + i];
-        }
-      }
-      for (j = i; j < Ny; j++) a[j*Nx + i] *= g;
+          for (k = i; k < Ny; k++) { a[k*Nx + j] += f*a[k*Nx + i]; }
+        }
+      }
+      for (j = i; j < Ny; j++) { a[j*Nx + i] *= g; }
     } else {
-      for (j = i; j < Ny; j++) a[j*Nx + i] = 0.0;
-    }
-    ++a[i*Nx + i];
+      for (j = i; j < Ny; j++) { a[j*Nx + i] = 0.0; }
+    }
+    a[i*Nx + i] ++;
   }
 
@@ -126 +132 @@
         c = 0.0;
         s = 1.0;
-        for (i = l; i < k; i++) {
+        for (i = l; i <= k; i++) {
           f = s*rv1[i];
           if (fabs(f)+anorm != anorm) {
@@ -160 +166 @@
       f = ((y-z)*(y+z) + (g-h)*(g+h))/(2.0*h*y);
       g = hypot (f, 1.0);
-      f = ((x-z)*(x+z) + h*((y/(f+FSIGN(g, f)))-h))/x;
+      
+      float tmp = g;
+      if (f < 0.0) { tmp = - tmp; }
+
+      f = ((x-z)*(x+z) + h*((y/(f+tmp))-h))/x;
+
       c = s = 1.0;
-      for (j = l; j < nm; j++) {
+      for (j = l; j <= nm; j++) {
         i = j+1;
         g = rv1[i];
@@ -184 +195 @@
         z = hypot(f, h);
         w[j] = z;
-        if (z) {
+        if (z != 0.0) {
           z = 1.0/z;
           c = f*z;