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; /**

@@ -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(min(m+1,n)); + int nu = std::min(m,n); + s = Vector(std::min(m+1,n)); U = Matrix(m, nu, Real(0)); V = Matrix(n,n); Vector 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 &A) { - int minm = min(m+1,n); + int minm = std::min(m+1,n); A = Matrix(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 // 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 LU_; + TNT::Array2D LU_; int m, n, pivsign; - Array1D piv; + TNT::Array1D piv; - Array2D permute_copy(const Array2D &A, - const Array1D &pivot, int j0, int j1) + TNT::Array2D permute_copy(const TNT::Array2D &A, + const TNT::Array1D &pivot, int j0, int j1) { int pivotlength = pivot.dim(); - Array2D X(pivotlength, j1-j0+1); + TNT::Array2D X(pivotlength, j1-j0+1); for (int i = 0; i < pivotlength; i++) @@ -54,14 +51,14 @@ return X; } - Array1D permute_copy(const Array1D &A, - const Array1D &pivot) + TNT::Array1D permute_copy(const TNT::Array1D &A, + const TNT::Array1D &pivot) { int pivotlength = pivot.dim(); if (pivotlength != A.dim()) - return Array1D(); + return TNT::Array1D(); - Array1D x(pivotlength); + TNT::Array1D 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 &A) : LU_(A.copy()), m(A.dim1()), n(A.dim2()), + LU (const TNT::Array2D &A) : LU_(A.copy()), m(A.dim1()), n(A.dim2()), piv(A.dim1()) { @@ -91,7 +88,7 @@ } pivsign = 1; Real *LUrowi = 0;; - Array1D LUcolj(m); + TNT::Array1D 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 getL () { - Array2D L_(m,n); + TNT::Array2D getL () { + TNT::Array2D 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 getU () { - Array2D U_(n,n); + TNT::Array2D getU () { + TNT::Array2D 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 getPivot () { + TNT::Array1D getPivot () { return piv; } @@ -232,23 +229,23 @@ 0x0 (null) array. */ - Array2D solve (const Array2D &B) + TNT::Array2D solve (const TNT::Array2D &B) { /* Dimensions: A is mxn, X is nxk, B is mxk */ if (B.dim1() != m) { - return Array2D(0,0); + return TNT::Array2D(0,0); } if (!isNonsingular()) { - return Array2D(0,0); + return TNT::Array2D(0,0); } // Copy right hand side with pivoting int nx = B.dim2(); - Array2D X = permute_copy(B, piv, 0, nx-1); + TNT::Array2D 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 solve (const Array1D &b) + TNT::Array1D solve (const TNT::Array1D &b) { /* Dimensions: A is mxn, X is nxk, B is mxk */ if (b.dim1() != m) { - return Array1D(); + return TNT::Array1D(); } if (!isNonsingular()) { - return Array1D(); + return TNT::Array1D(); } - Array1D x = permute_copy(b, piv); + TNT::Array1D 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; - /**

For a symmetric, positive definite matrix A, this function @@ -21,8 +19,8 @@

Typical usage looks like: 

-	Array2D A(n,n);
-	Array2D L;
+	TNT::Array2D A(n,n);
+	TNT::Array2D L;
 
 	 ...
 
@@ -46,16 +44,16 @@
 template 
 class Cholesky
 {
-	Array2D L_;		// lower triangular factor
+	TNT::Array2D L_;		// lower triangular factor
 	int isspd;				// 1 if matrix to be factored was SPD
 
 public:
 
 	Cholesky();
-	Cholesky(const Array2D &A);
-	Array2D getL() const;
-	Array1D solve(const Array1D &B);
-	Array2D solve(const Array2D &B);
+	Cholesky(const TNT::Array2D &A);
+	TNT::Array2D getL() const;
+	TNT::Array1D solve(const TNT::Array1D &B);
+	TNT::Array2D solve(const TNT::Array2D &B);
 	int is_spd() const;
 
 };
@@ -77,7 +75,7 @@
 	@return the lower triangular factor, L, such that L*L'=A.
 */
 template 
-Array2D Cholesky::getL() const
+TNT::Array2D Cholesky::getL() const
 {
 	return L_;
 }
@@ -89,7 +87,7 @@
 	evalutate true (1) then the factorizaiton was successful.
 */
 template 
-Cholesky::Cholesky(const Array2D &A)
+Cholesky::Cholesky(const TNT::Array2D &A)
 {
 
 
@@ -100,11 +98,11 @@
 
 	if (m != n)
 	{
-		L_ = Array2D(0,0);
+		L_ = TNT::Array2D(0,0);
 		return;
 	}
 
-	L_ = Array2D(n,n);
+	L_ = TNT::Array2D(n,n);
 
 
       // Main loop.
@@ -143,14 +141,14 @@
    						array is returned.
 */
 template 
-Array1D Cholesky::solve(const Array1D &b)
+TNT::Array1D Cholesky::solve(const TNT::Array1D &b)
 {
 	int n = L_.dim1();
 	if (b.dim1() != n)
-		return Array1D();
+		return TNT::Array1D();
 
 
-	Array1D x = b.copy();
+	TNT::Array1D x = b.copy();
 
 
       // Solve L*y = b;
@@ -185,14 +183,14 @@
    						array is returned.
 */
 template 
-Array2D Cholesky::solve(const Array2D &B)
+TNT::Array2D Cholesky::solve(const TNT::Array2D &B)
 {
 	int n = L_.dim1();
 	if (B.dim1() != n)
-		return Array2D();
+		return TNT::Array2D();
 
 
-	Array2D X = B.copy();
+	TNT::Array2D 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 
 // for abs() below
 
-using namespace TNT;
-using namespace std;
 
 namespace JAMA
 {
@@ -40,25 +38,25 @@
 {
 
 
-	Array2D U, V;
-	Array1D s;
+	TNT::Array2D U, V;
+	TNT::Array1D s;
 	int m, n;
 
   public:
 
 
-   SVD (const Array2D &Arg) {
+   SVD (const TNT::Array2D &Arg) {
 
 
       m = Arg.dim1();
       n = Arg.dim2();
-      int nu = min(m,n);
-      s = Array1D(min(m+1,n));
-      U = Array2D(m, nu, Real(0));
-      V = Array2D(n,n);
-      Array1D e(n);
-      Array1D work(m);
-	  Array2D A(Arg.copy());
+      int nu = std::min(m,n);
+      s = TNT::Array1D(std::min(m+1,n));
+      U = TNT::Array2D(m, nu, Real(0));
+      V = TNT::Array2D(n,n);
+      TNT::Array1D e(n);
+      TNT::Array1D work(m);
+	  TNT::Array2D 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 &A)
+   void getU (TNT::Array2D &A)
    {
-   	  int minm = min(m+1,n);
+   	  int minm = std::min(m+1,n);
 
-	  A = Array2D(m, minm);
+	  A = TNT::Array2D(m, minm);
 
 	  for (int i=0; i &A)
+   void getV (TNT::Array2D &A)
    {
    	  A = V;
    }
 
    /** Return the one-dimensional array of singular values */
 
-   void getSingularValues (Array1D &x)
+   void getSingularValues (TNT::Array1D &x)
    {
       x = s;
    }
@@ -493,8 +491,8 @@
    @return     S
    */
 
-   void getS (Array2D &A) {
-   	  A = Array2D(n,n);
+   void getS (TNT::Array2D &A) {
+   	  A = TNT::Array2D(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) {