Index: tnt_linalg.h
===================================================================
--- tnt_linalg.h	(revision 37575)
+++ tnt_linalg.h	(revision 37576)
@@ -11,7 +11,6 @@
 namespace TNT {
 namespace Linear_Algebra {
 
-	using namespace std;
 
 /**
    <P>
@@ -442,12 +441,12 @@
 
       Real f = Real(0.0);
       Real tst1 = Real(0.0);
-      Real eps = pow(2.0,-52.0);
+      Real eps = std::pow(2.0,-52.0);
       for (int l = 0; l < n; l++) {
 
          // Find small subdiagonal element
 
-         tst1 = max(tst1,abs(d[l]) + abs(e[l]));
+         tst1 = std::max(tst1,abs(d[l]) + abs(e[l]));
          int m = l;
 
         // Original while-loop from Java code
@@ -677,7 +676,7 @@
       int n = nn-1;
       int low = 0;
       int high = nn-1;
-      Real eps = pow(2.0,-52.0);
+      Real eps = std::pow(2.0,-52.0);
       Real exshift = Real(0.0);
       Real p=0,q=0,r=0,s=0,z=0,t,w,x,y;
 
@@ -689,7 +688,7 @@
             d[i] = H[i][i];
             e[i] = Real(0.0);
          }
-         for (int j = max(i-Real(1.0)); j < nn; j++) {
+         for (int j = std::max(i-Real(1.0)); j < nn; j++) {
             norm = norm + abs(H[i][j]);
          }
       }
@@ -920,7 +919,7 @@
 
                   // Column modification
 
-                  for (int i = 0; i <= min(n,k+3); i++) {
+                  for (int i = 0; i <= std::min(n,k+3); i++) {
                      p = x * H[i][k] + y * H[i][k+1];
                      if (notlast) {
                         p = p + z * H[i][k+2];
@@ -1069,7 +1068,7 @@
 
                   // Overflow control
 
-                  t = max(abs(H[i][n-1]),abs(H[i][n]));
+                  t = std::max(abs(H[i][n-1]),abs(H[i][n]));
                   if ((eps * t) * t > 1) {
                      for (int j = i; j <= n; j++) {
                         H[j][n-1] = H[j][n-1] / t;
@@ -1096,7 +1095,7 @@
       for (int j = nn-1; j >= low; j--) {
          for (int i = low; i <= high; i++) {
             z = Real(0.0);
-            for (int k = low; k <= min(j,high); k++) {
+            for (int k = low; k <= std::min(j,high); k++) {
                z = z + V[i][k] * H[k][j];
             }
             V[i][j] = z;
@@ -1330,7 +1329,7 @@
 
             // Most of the time is spent in the following dot product.
 
-            int kmax = min(i,j);
+            int kmax = std::min(i,j);
             double s = Real(0.0);
             for (int k = 0; k < kmax; k++) {
                s += LUrowi[k]*LUcolj[k];
@@ -1881,8 +1880,8 @@
 
       m = Arg.num_rows();
       n = Arg.num_cols();
-      int nu = min(m,n);
-      s = Vector<Real>(min(m+1,n));
+      int nu = std::min(m,n);
+      s = Vector<Real>(std::min(m+1,n));
       U = Matrix<Real>(m, nu, Real(0));
       V = Matrix<Real>(n,n);
       Vector<Real> e(n);
@@ -1895,9 +1894,9 @@
       // Reduce A to bidiagonal form, storing the diagonal elements
       // in s and the super-diagonal elements in e.
 
-      int nct = min(m-1,n);
-      int nrt = max(0,min(n-2,m));
-      for (k = 0; k < max(nct,nrt); k++) {
+      int nct = std::min(m-1,n);
+      int nrt = std::max(0,std::min(n-2,m));
+      for (k = 0; k < std::max(nct,nrt); k++) {
          if (k < nct) {
 
             // Compute the transformation for the k-th column and
@@ -1999,7 +1998,7 @@
 
       // Set up the final bidiagonal matrix or order p.
 
-      int p = min(n,m+1);
+      int p = std::min(n,m+1);
       if (nct < n) {
          s[nct] = A[nct][nct];
       }
@@ -2075,7 +2074,7 @@
 
       int pp = p-1;
       int iter = 0;
-      double eps = pow(2.0,-52.0);
+      double eps = std::pow(2.0,-52.0);
       while (p > 0) {
          int k=0;
 		 int kase=0;
@@ -2185,7 +2184,7 @@
 
                // Calculate the shift.
 
-               double scale = max(max(max(max(
+               double scale = std::max(std::max(std::max(std::max(
                        abs(s[p-1]),abs(s[p-2])),abs(e[p-2])),
                        abs(s[k])),abs(e[k]));
                double sp = s[p-1]/scale;
@@ -2294,7 +2293,7 @@
 
    void getU (Matrix<Real> &A)
    {
-   	  int minm = min(m+1,n);
+          int minm = std::min(m+1,n);
 
 	  A = Matrix<Real>(m, minm);
 
@@ -2332,13 +2331,13 @@
       }
    }
 
-   /** Two norm  (max(S)) */
+   /** Two norm  (std::max(S)) */
 
    double norm2 () {
       return s[0];
    }
 
-   /** Two norm of condition number (max(S)/min(S)) */
+   /** Two norm of condition number (std::max(S)/std::min(S)) */
 
    double cond () {
       return s[0]/s[min(m,n)-1];
@@ -2350,8 +2349,8 @@
 
    int rank ()
    {
-      double eps = pow(2.0,-52.0);
-      double tol = max(m,n)*s[0]*eps;
+      double eps = std::pow(2.0,-52.0);
+      double tol = std::max(m,n)*s[0]*eps;
       int r = 0;
       for (int i = 0; i < s.dim(); i++) {
          if (s[i] > tol) {
Index: jama_eig.h
===================================================================
--- jama_eig.h	(revision 37575)
+++ jama_eig.h	(revision 37576)
@@ -12,8 +12,6 @@
 #include <cmath>
 // for abs() below
 
-using namespace TNT;
-using namespace std;
 
 namespace JAMA
 {
@@ -120,7 +118,7 @@
          Real scale = 0.0;
          Real h = 0.0;
          for (int k = 0; k < i; k++) {
-            scale = scale + abs(d[k]);
+            scale = scale + std::abs(d[k]);
          }
          if (scale == 0.0) {
             e[i] = d[i-1];
@@ -231,17 +229,17 @@
 
       Real f = 0.0;
       Real tst1 = 0.0;
-      Real eps = pow(2.0,-52.0);
+      Real eps = std::pow(2.0,-52.0);
       for (int l = 0; l < n; l++) {
 
          // Find small subdiagonal element
 
-         tst1 = max(tst1,abs(d[l]) + abs(e[l]));
+         tst1 = std::max(tst1,std::abs(d[l]) + std::abs(e[l]));
          int m = l;
 
         // Original while-loop from Java code
          while (m < n) {
-            if (abs(e[m]) <= eps*tst1) {
+            if (std::abs(e[m]) <= eps*tst1) {
                break;
             }
             m++;
@@ -309,7 +307,7 @@
 
                // Check for convergence.
 
-            } while (abs(e[l]) > eps*tst1);
+            } while (std::abs(e[l]) > eps*tst1);
          }
          d[l] = d[l] + f;
          e[l] = 0.0;
@@ -356,7 +354,7 @@
 
          Real scale = 0.0;
          for (int i = m; i <= high; i++) {
-            scale = scale + abs(H[i][m-1]);
+            scale = scale + std::abs(H[i][m-1]);
          }
          if (scale != 0.0) {
 
@@ -437,7 +435,7 @@
    Real cdivr, cdivi;
    void cdiv(Real xr, Real xi, Real yr, Real yi) {
       Real r,d;
-      if (abs(yr) > abs(yi)) {
+      if (std::abs(yr) > std::abs(yi)) {
          r = yi/yr;
          d = yr + r*yi;
          cdivr = (xr + r*xi)/d;
@@ -466,7 +464,7 @@
       int n = nn-1;
       int low = 0;
       int high = nn-1;
-      Real eps = pow(2.0,-52.0);
+      Real eps = std::pow(2.0,-52.0);
       Real exshift = 0.0;
       Real p=0,q=0,r=0,s=0,z=0,t,w,x,y;
 
@@ -478,8 +476,8 @@
             d[i] = H[i][i];
             e[i] = 0.0;
          }
-         for (int j = max(i-1,0); j < nn; j++) {
-            norm = norm + abs(H[i][j]);
+         for (int j = std::max(i-1,0); j < nn; j++) {
+            norm = norm + std::abs(H[i][j]);
          }
       }
 
@@ -492,11 +490,11 @@
 
          int l = n;
          while (l > low) {
-            s = abs(H[l-1][l-1]) + abs(H[l][l]);
+            s = std::abs(H[l-1][l-1]) + std::abs(H[l][l]);
             if (s == 0.0) {
                s = norm;
             }
-            if (abs(H[l][l-1]) < eps * s) {
+            if (std::abs(H[l][l-1]) < eps * s) {
                break;
             }
             l--;
@@ -518,7 +516,7 @@
             w = H[n][n-1] * H[n-1][n];
             p = (H[n-1][n-1] - H[n][n]) / 2.0;
             q = p * p + w;
-            z = sqrt(abs(q));
+            z = sqrt(std::abs(q));
             H[n][n] = H[n][n] + exshift;
             H[n-1][n-1] = H[n-1][n-1] + exshift;
             x = H[n][n];
@@ -539,7 +537,7 @@
                e[n-1] = 0.0;
                e[n] = 0.0;
                x = H[n][n-1];
-               s = abs(x) + abs(z);
+               s = std::abs(x) + std::abs(z);
                p = x / s;
                q = z / s;
                r = sqrt(p * p+q * q);
@@ -602,7 +600,7 @@
                for (int i = low; i <= n; i++) {
                   H[i][i] -= x;
                }
-               s = abs(H[n][n-1]) + abs(H[n-1][n-2]);
+               s = std::abs(H[n][n-1]) + std::abs(H[n-1][n-2]);
                x = y = 0.75 * s;
                w = -0.4375 * s * s;
             }
@@ -638,16 +636,16 @@
                p = (r * s - w) / H[m+1][m] + H[m][m+1];
                q = H[m+1][m+1] - z - r - s;
                r = H[m+2][m+1];
-               s = abs(p) + abs(q) + abs(r);
+               s = std::abs(p) + std::abs(q) + std::abs(r);
                p = p / s;
                q = q / s;
                r = r / s;
                if (m == l) {
                   break;
                }
-               if (abs(H[m][m-1]) * (abs(q) + abs(r)) <
-                  eps * (abs(p) * (abs(H[m-1][m-1]) + abs(z) +
-                  abs(H[m+1][m+1])))) {
+               if (std::abs(H[m][m-1]) * (std::abs(q) + std::abs(r)) <
+                  eps * (std::abs(p) * (std::abs(H[m-1][m-1]) + std::abs(z) +
+                  std::abs(H[m+1][m+1])))) {
                      break;
                }
                m--;
@@ -668,7 +666,7 @@
                   p = H[k][k-1];
                   q = H[k+1][k-1];
                   r = (notlast ? H[k+2][k-1] : 0.0);
-                  x = abs(p) + abs(q) + abs(r);
+                  x = std::abs(p) + std::abs(q) + std::abs(r);
                   if (x != 0.0) {
                      p = p / x;
                      q = q / x;
@@ -709,7 +707,7 @@
 
                   // Column modification
 
-                  for (int i = 0; i <= min(n,k+3); i++) {
+                  for (int i = 0; i <= std::min(n,k+3); i++) {
                      p = x * H[i][k] + y * H[i][k+1];
                      if (notlast) {
                         p = p + z * H[i][k+2];
@@ -776,7 +774,7 @@
                      q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
                      t = (x * s - z * r) / q;
                      H[i][n] = t;
-                     if (abs(x) > abs(z)) {
+                     if (std::abs(x) > std::abs(z)) {
                         H[i+1][n] = (-r - w * t) / x;
                      } else {
                         H[i+1][n] = (-s - y * t) / z;
@@ -785,7 +783,7 @@
 
                   // Overflow control
 
-                  t = abs(H[i][n]);
+                  t = std::abs(H[i][n]);
                   if ((eps * t) * t > 1) {
                      for (int j = i; j <= n; j++) {
                         H[j][n] = H[j][n] / t;
@@ -801,7 +799,7 @@
 
             // Last vector component imaginary so matrix is triangular
 
-            if (abs(H[n][n-1]) > abs(H[n-1][n])) {
+            if (std::abs(H[n][n-1]) > std::abs(H[n-1][n])) {
                H[n-1][n-1] = q / H[n][n-1];
                H[n-1][n] = -(H[n][n] - p) / H[n][n-1];
             } else {
@@ -840,13 +838,13 @@
                      vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
                      vi = (d[i] - p) * 2.0 * q;
                      if ((vr == 0.0) && (vi == 0.0)) {
-                        vr = eps * norm * (abs(w) + abs(q) +
-                        abs(x) + abs(y) + abs(z));
+                        vr = eps * norm * (std::abs(w) + std::abs(q) +
+                        std::abs(x) + std::abs(y) + std::abs(z));
                      }
                      cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
                      H[i][n-1] = cdivr;
                      H[i][n] = cdivi;
-                     if (abs(x) > (abs(z) + abs(q))) {
+                     if (std::abs(x) > (std::abs(z) + std::abs(q))) {
                         H[i+1][n-1] = (-ra - w * H[i][n-1] + q * H[i][n]) / x;
                         H[i+1][n] = (-sa - w * H[i][n] - q * H[i][n-1]) / x;
                      } else {
@@ -858,7 +856,7 @@
 
                   // Overflow control
 
-                  t = max(abs(H[i][n-1]),abs(H[i][n]));
+                  t = std::max(std::abs(H[i][n-1]),std::abs(H[i][n]));
                   if ((eps * t) * t > 1) {
                      for (int j = i; j <= n; j++) {
                         H[j][n-1] = H[j][n-1] / t;
@@ -885,7 +883,7 @@
       for (int j = nn-1; j >= low; j--) {
          for (int i = low; i <= high; i++) {
             z = 0.0;
-            for (int k = low; k <= min(j,high); k++) {
+            for (int k = low; k <= std::min(j,high); k++) {
                z = z + V[i][k] * H[k][j];
             }
             V[i][j] = z;
Index: jama_lu.h
===================================================================
--- jama_lu.h	(revision 37575)
+++ jama_lu.h	(revision 37576)
@@ -9,9 +9,6 @@
 #define NEAR_ZERO(val, epsilon)(((val) > -epsilon) && ((val) < epsilon))
 
 
-using namespace TNT;
-using namespace std;
-
 namespace JAMA
 {
 
@@ -34,17 +31,17 @@
 
 
    /* Array for internal storage of decomposition.  */
-   Array2D<Real>  LU_;
+   TNT::Array2D<Real>  LU_;
    int m, n, pivsign;
-   Array1D<int> piv;
+   TNT::Array1D<int> piv;
 
 
-   Array2D<Real> permute_copy(const Array2D<Real> &A,
-   			const Array1D<int> &pivot, int j0, int j1)
+   TNT::Array2D<Real> permute_copy(const TNT::Array2D<Real> &A,
+   			const TNT::Array1D<int> &pivot, int j0, int j1)
 	{
 		int pivotlength = pivot.dim();
 
-		Array2D<Real> X(pivotlength, j1-j0+1);
+		TNT::Array2D<Real> X(pivotlength, j1-j0+1);
 
 
          for (int i = 0; i < pivotlength; i++)
@@ -54,14 +51,14 @@
 		return X;
 	}
 
-   Array1D<Real> permute_copy(const Array1D<Real> &A,
-   		const Array1D<int> &pivot)
+   TNT::Array1D<Real> permute_copy(const TNT::Array1D<Real> &A,
+   		const TNT::Array1D<int> &pivot)
 	{
 		int pivotlength = pivot.dim();
 		if (pivotlength != A.dim())
-			return Array1D<Real>();
+			return TNT::Array1D<Real>();
 
-		Array1D<Real> x(pivotlength);
+		TNT::Array1D<Real> x(pivotlength);
 
 
          for (int i = 0; i < pivotlength; i++)
@@ -78,7 +75,7 @@
    @return     LU Decomposition object to access L, U and piv.
    */
 
-    LU (const Array2D<Real> &A) : LU_(A.copy()), m(A.dim1()), n(A.dim2()),
+    LU (const TNT::Array2D<Real> &A) : LU_(A.copy()), m(A.dim1()), n(A.dim2()),
 		piv(A.dim1())
 
 	{
@@ -91,7 +88,7 @@
       }
       pivsign = 1;
       Real *LUrowi = 0;;
-      Array1D<Real> LUcolj(m);
+      TNT::Array1D<Real> LUcolj(m);
 
       // Outer loop.
 
@@ -110,7 +107,7 @@
 
             // Most of the time is spent in the following dot product.
 
-            int kmax = min(i,j);
+            int kmax = std::min(i,j);
             double s = 0.0;
             for (int k = 0; k < kmax; k++) {
                s += LUrowi[k]*LUcolj[k];
@@ -123,7 +120,7 @@
 
          int p = j;
          for (int i = j+1; i < m; i++) {
-            if (abs(LUcolj[i]) > abs(LUcolj[p])) {
+            if (std::abs(LUcolj[i]) > std::abs(LUcolj[p])) {
                p = i;
             }
          }
@@ -168,8 +165,8 @@
    @return     L
    */
 
-   Array2D<Real> getL () {
-      Array2D<Real> L_(m,n);
+   TNT::Array2D<Real> getL () {
+      TNT::Array2D<Real> L_(m,n);
       for (int i = 0; i < m; i++) {
          for (int j = 0; j < n; j++) {
             if (i > j) {
@@ -188,8 +185,8 @@
    @return     U portion of LU factorization.
    */
 
-   Array2D<Real> getU () {
-   	  Array2D<Real> U_(n,n);
+   TNT::Array2D<Real> getU () {
+   	  TNT::Array2D<Real> U_(n,n);
       for (int i = 0; i < n; i++) {
          for (int j = 0; j < n; j++) {
             if (i <= j) {
@@ -206,7 +203,7 @@
    @return     piv
    */
 
-   Array1D<int> getPivot () {
+   TNT::Array1D<int> getPivot () {
       return piv;
    }
 
@@ -232,23 +229,23 @@
    					0x0 (null) array.
    */
 
-   Array2D<Real> solve (const Array2D<Real> &B)
+   TNT::Array2D<Real> solve (const TNT::Array2D<Real> &B)
    {
 
 	  /* Dimensions: A is mxn, X is nxk, B is mxk */
 
       if (B.dim1() != m) {
-	  	return Array2D<Real>(0,0);
+	  	return TNT::Array2D<Real>(0,0);
       }
       if (!isNonsingular()) {
-        return Array2D<Real>(0,0);
+        return TNT::Array2D<Real>(0,0);
       }
 
       // Copy right hand side with pivoting
       int nx = B.dim2();
 
 
-	  Array2D<Real> X = permute_copy(B, piv, 0, nx-1);
+	  TNT::Array2D<Real> X = permute_copy(B, piv, 0, nx-1);
 
       // Solve L*Y = B(piv,:)
       for (int k = 0; k < n; k++) {
@@ -276,26 +273,26 @@
    /** Solve A*x = b, where x and b are vectors of length equal
    		to the number of rows in A.
 
-   @param  b   a vector (Array1D> of length equal to the first dimension
+   @param  b   a vector (TNT::Array1D> of length equal to the first dimension
    						of A.
-   @return x a vector (Array1D> so that L*U*x = b(piv), if B is nonconformant,
+   @return x a vector (TNT::Array1D> so that L*U*x = b(piv), if B is nonconformant,
    					returns 0x0 (null) array.
    */
 
-   Array1D<Real> solve (const Array1D<Real> &b)
+   TNT::Array1D<Real> solve (const TNT::Array1D<Real> &b)
    {
 
 	  /* Dimensions: A is mxn, X is nxk, B is mxk */
 
       if (b.dim1() != m) {
-	  	return Array1D<Real>();
+	  	return TNT::Array1D<Real>();
       }
       if (!isNonsingular()) {
-        return Array1D<Real>();
+        return TNT::Array1D<Real>();
       }
 
 
-	  Array1D<Real> x = permute_copy(b, piv);
+	  TNT::Array1D<Real> x = permute_copy(b, piv);
 
       // Solve L*Y = B(piv)
       for (int k = 0; k < n; k++) {
Index: jama_cholesky.h
===================================================================
--- jama_cholesky.h	(revision 37575)
+++ jama_cholesky.h	(revision 37576)
@@ -8,8 +8,6 @@
 namespace JAMA
 {
 
-using namespace TNT;
-
 /**
    <P>
    For a symmetric, positive definite matrix A, this function
@@ -21,8 +19,8 @@
 
    <p>Typical usage looks like:
    <pre>
-	Array2D<double> A(n,n);
-	Array2D<double> L;
+	TNT::Array2D<double> A(n,n);
+	TNT::Array2D<double> L;
 
 	 ...
 
@@ -46,16 +44,16 @@
 template <class Real>
 class Cholesky
 {
-	Array2D<Real> L_;		// lower triangular factor
+	TNT::Array2D<Real> L_;		// lower triangular factor
 	int isspd;				// 1 if matrix to be factored was SPD
 
 public:
 
 	Cholesky();
-	Cholesky(const Array2D<Real> &A);
-	Array2D<Real> getL() const;
-	Array1D<Real> solve(const Array1D<Real> &B);
-	Array2D<Real> solve(const Array2D<Real> &B);
+	Cholesky(const TNT::Array2D<Real> &A);
+	TNT::Array2D<Real> getL() const;
+	TNT::Array1D<Real> solve(const TNT::Array1D<Real> &B);
+	TNT::Array2D<Real> solve(const TNT::Array2D<Real> &B);
 	int is_spd() const;
 
 };
@@ -77,7 +75,7 @@
 	@return the lower triangular factor, L, such that L*L'=A.
 */
 template <class Real>
-Array2D<Real> Cholesky<Real>::getL() const
+TNT::Array2D<Real> Cholesky<Real>::getL() const
 {
 	return L_;
 }
@@ -89,7 +87,7 @@
 	evalutate true (1) then the factorizaiton was successful.
 */
 template <class Real>
-Cholesky<Real>::Cholesky(const Array2D<Real> &A)
+Cholesky<Real>::Cholesky(const TNT::Array2D<Real> &A)
 {
 
 
@@ -100,11 +98,11 @@
 
 	if (m != n)
 	{
-		L_ = Array2D<Real>(0,0);
+		L_ = TNT::Array2D<Real>(0,0);
 		return;
 	}
 
-	L_ = Array2D<Real>(n,n);
+	L_ = TNT::Array2D<Real>(n,n);
 
 
       // Main loop.
@@ -143,14 +141,14 @@
    						array is returned.
 */
 template <class Real>
-Array1D<Real> Cholesky<Real>::solve(const Array1D<Real> &b)
+TNT::Array1D<Real> Cholesky<Real>::solve(const TNT::Array1D<Real> &b)
 {
 	int n = L_.dim1();
 	if (b.dim1() != n)
-		return Array1D<Real>();
+		return TNT::Array1D<Real>();
 
 
-	Array1D<Real> x = b.copy();
+	TNT::Array1D<Real> x = b.copy();
 
 
       // Solve L*y = b;
@@ -185,14 +183,14 @@
    						array is returned.
 */
 template <class Real>
-Array2D<Real> Cholesky<Real>::solve(const Array2D<Real> &B)
+TNT::Array2D<Real> Cholesky<Real>::solve(const TNT::Array2D<Real> &B)
 {
 	int n = L_.dim1();
 	if (B.dim1() != n)
-		return Array2D<Real>();
+		return TNT::Array2D<Real>();
 
 
-	Array2D<Real> X = B.copy();
+	TNT::Array2D<Real> X = B.copy();
 	int nx = B.dim2();
 
 // Cleve's original code
Index: jama_svd.h
===================================================================
--- jama_svd.h	(revision 37575)
+++ jama_svd.h	(revision 37576)
@@ -13,8 +13,6 @@
 #include <cmath>
 // for abs() below
 
-using namespace TNT;
-using namespace std;
 
 namespace JAMA
 {
@@ -40,25 +38,25 @@
 {
 
 
-	Array2D<Real> U, V;
-	Array1D<Real> s;
+	TNT::Array2D<Real> U, V;
+	TNT::Array1D<Real> s;
 	int m, n;
 
   public:
 
 
-   SVD (const Array2D<Real> &Arg) {
+   SVD (const TNT::Array2D<Real> &Arg) {
 
 
       m = Arg.dim1();
       n = Arg.dim2();
-      int nu = min(m,n);
-      s = Array1D<Real>(min(m+1,n));
-      U = Array2D<Real>(m, nu, Real(0));
-      V = Array2D<Real>(n,n);
-      Array1D<Real> e(n);
-      Array1D<Real> work(m);
-	  Array2D<Real> A(Arg.copy());
+      int nu = std::min(m,n);
+      s = TNT::Array1D<Real>(std::min(m+1,n));
+      U = TNT::Array2D<Real>(m, nu, Real(0));
+      V = TNT::Array2D<Real>(n,n);
+      TNT::Array1D<Real> e(n);
+      TNT::Array1D<Real> work(m);
+	  TNT::Array2D<Real> A(Arg.copy());
       int wantu = 1;  					/* boolean */
       int wantv = 1;  					/* boolean */
 	  int i=0, j=0, k=0;
@@ -66,9 +64,9 @@
       // Reduce A to bidiagonal form, storing the diagonal elements
       // in s and the super-diagonal elements in e.
 
-      int nct = min(m-1,n);
-      int nrt = max(0,min(n-2,m));
-      for (k = 0; k < max(nct,nrt); k++) {
+      int nct = std::min(m-1,n);
+      int nrt = std::max(0,std::min(n-2,m));
+      for (k = 0; k < std::max(nct,nrt); k++) {
          if (k < nct) {
 
             // Compute the transformation for the k-th column and
@@ -170,7 +168,7 @@
 
       // Set up the final bidiagonal matrix or order p.
 
-      int p = min(n,m+1);
+      int p = std::min(n,m+1);
       if (nct < n) {
          s[nct] = A[nct][nct];
       }
@@ -246,7 +244,7 @@
 
       int pp = p-1;
       int iter = 0;
-      Real eps(pow(2.0,-52.0));
+      Real eps(std::pow(2.0,-52.0));
       while (p > 0) {
          int k=0;
 		 int kase=0;
@@ -267,7 +265,7 @@
             if (k == -1) {
                break;
             }
-            if (abs(e[k]) <= eps*(abs(s[k]) + abs(s[k+1]))) {
+            if (std::abs(e[k]) <= eps*(std::abs(s[k]) + std::abs(s[k+1]))) {
                e[k] = 0.0;
                break;
             }
@@ -280,9 +278,9 @@
                if (ks == k) {
                   break;
                }
-               Real t( (ks != p ? abs(e[ks]) : 0.) +
-                          (ks != k+1 ? abs(e[ks-1]) : 0.));
-               if (abs(s[ks]) <= eps*t)  {
+               Real t( (ks != p ? std::abs(e[ks]) : 0.) +
+                          (ks != k+1 ? std::abs(e[ks-1]) : 0.));
+               if (std::abs(s[ks]) <= eps*t)  {
                   s[ks] = 0.0;
                   break;
                }
@@ -356,9 +354,9 @@
 
                // Calculate the shift.
 
-               Real scale = max(max(max(max(
-                       abs(s[p-1]),abs(s[p-2])),abs(e[p-2])),
-                       abs(s[k])),abs(e[k]));
+               Real scale = std::max(std::max(std::max(std::max(
+                       std::abs(s[p-1]),std::abs(s[p-2])),std::abs(e[p-2])),
+                       std::abs(s[k])),std::abs(e[k]));
                Real sp = s[p-1]/scale;
                Real spm1 = s[p-2]/scale;
                Real epm1 = e[p-2]/scale;
@@ -463,11 +461,11 @@
    }
 
 
-   void getU (Array2D<Real> &A)
+   void getU (TNT::Array2D<Real> &A)
    {
-   	  int minm = min(m+1,n);
+   	  int minm = std::min(m+1,n);
 
-	  A = Array2D<Real>(m, minm);
+	  A = TNT::Array2D<Real>(m, minm);
 
 	  for (int i=0; i<m; i++)
 	  	for (int j=0; j<minm; j++)
@@ -477,14 +475,14 @@
 
    /* Return the right singular vectors */
 
-   void getV (Array2D<Real> &A)
+   void getV (TNT::Array2D<Real> &A)
    {
    	  A = V;
    }
 
    /** Return the one-dimensional array of singular values */
 
-   void getSingularValues (Array1D<Real> &x)
+   void getSingularValues (TNT::Array1D<Real> &x)
    {
       x = s;
    }
@@ -493,8 +491,8 @@
    @return     S
    */
 
-   void getS (Array2D<Real> &A) {
-   	  A = Array2D<Real>(n,n);
+   void getS (TNT::Array2D<Real> &A) {
+   	  A = TNT::Array2D<Real>(n,n);
       for (int i = 0; i < n; i++) {
          for (int j = 0; j < n; j++) {
             A[i][j] = 0.0;
@@ -503,16 +501,16 @@
       }
    }
 
-   /** Two norm  (max(S)) */
+   /** Two norm  (std::max(S)) */
 
    Real norm2 () {
       return s[0];
    }
 
-   /** Two norm of condition number (max(S)/min(S)) */
+   /** Two norm of condition number (std::max(S)/std::min(S)) */
 
    Real cond () {
-      return s[0]/s[min(m,n)-1];
+      return s[0]/s[std::min(m,n)-1];
    }
 
    /** Effective numerical matrix rank
@@ -521,8 +519,8 @@
 
    int rank ()
    {
-      Real eps = pow(2.0,-52.0);
-      Real tol = max(m,n)*s[0]*eps;
+      Real eps = std::pow(2.0,-52.0);
+      Real tol = std::max(m,n)*s[0]*eps;
       int r = 0;
       for (int i = 0; i < s.dim(); i++) {
          if (s[i] > tol) {
