Index: /branches/eam_branches/ipp-20220316/Ohana/src/opihi/lib.data/svdcmp.c
===================================================================
--- /branches/eam_branches/ipp-20220316/Ohana/src/opihi/lib.data/svdcmp.c	(revision 42162)
+++ /branches/eam_branches/ipp-20220316/Ohana/src/opihi/lib.data/svdcmp.c	(revision 42163)
@@ -23,4 +23,5 @@
   l = nm = 0;
   ALLOCATE (rv1, float, Nx);
+
   for (i = 0; i < Nx; i++) {
     l = i + 1;
@@ -29,5 +30,5 @@
     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,15 +36,17 @@
 	}
 	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,6 +54,6 @@
     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,15 +61,17 @@
 	}
 	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,13 +81,14 @@
   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,20 +99,20 @@
     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,5 +132,5 @@
 	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,7 +166,12 @@
       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,5 +195,5 @@
 	z = hypot(f, h);
 	w[j] = z;
-	if (z) {
+	if (z != 0.0) {
 	  z = 1.0/z;
 	  c = f*z;