Index: jama_eig.h =================================================================== --- jama_eig.h (revision 20) +++ jama_eig.h (revision 49) @@ -896,6 +896,8 @@ public: + Eigenvalue() {} // Added by NWM to help in wrapping + /** Check for symmetry, then construct the eigenvalue decomposition @param A Square real (non-complex) matrix */ Index: tnt_vec.h =================================================================== --- tnt_vec.h (revision 20) +++ tnt_vec.h (revision 49) @@ -150,8 +150,8 @@ iterator begin() { return v_;} iterator end() { return v_ + n_; } - const iterator begin() const { return v_;} - const iterator end() const { return v_ + n_; } + iterator begin() const { return v_;} + iterator end() const { return v_ + n_; } // destructor Index: jama_cholesky.h =================================================================== --- jama_cholesky.h (revision 20) +++ jama_cholesky.h (revision 49) @@ -110,10 +110,10 @@ // Main loop. for (int j = 0; j < n; j++) { - double d = 0.0; + Real d(0.0); for (int k = 0; k < j; k++) { - Real s = 0.0; + Real s(0.0); for (int i = 0; i < k; i++) { s += L_[k][i]*L_[j][i]; Index: tnt_sparse_matrix_csr.h =================================================================== --- tnt_sparse_matrix_csr.h (revision 20) +++ tnt_sparse_matrix_csr.h (revision 49) @@ -59,7 +59,7 @@ Sparse_Matrix_CompRow(const Sparse_Matrix_CompRow &S); Sparse_Matrix_CompRow(int M, int N, int nz, const T *val, - const int *r, const int *c); + int *r, int *c); @@ -93,7 +93,7 @@ */ template Sparse_Matrix_CompRow::Sparse_Matrix_CompRow(int M, int N, int nz, - const T *val, const int *r, const int *c) : val_(nz,val), + const T *val, int *r, int *c) : val_(nz,val), rowptr_(M, r), colind_(nz, c), dim1_(M), dim2_(N) {} Index: jama_svd.h =================================================================== --- jama_svd.h (revision 20) +++ jama_svd.h (revision 49) @@ -94,7 +94,7 @@ // Apply the transformation. - double t = 0; + Real t(0.0); for (i = k; i < m; i++) { t += A[i][k]*A[i][j]; } @@ -150,7 +150,7 @@ } } for (j = k+1; j < n; j++) { - double t = -e[j]/e[k+1]; + Real t(-e[j]/e[k+1]); for (i = k+1; i < m; i++) { A[i][j] += t*work[i]; } @@ -194,7 +194,7 @@ for (k = nct-1; k >= 0; k--) { if (s[k] != 0.0) { for (j = k+1; j < nu; j++) { - double t = 0; + Real t(0.0); for (i = k; i < m; i++) { t += U[i][k]*U[i][j]; } @@ -225,7 +225,7 @@ for (k = n-1; k >= 0; k--) { if ((k < nrt) & (e[k] != 0.0)) { for (j = k+1; j < nu; j++) { - double t = 0; + Real t(0.0); for (i = k+1; i < n; i++) { t += V[i][k]*V[i][j]; } @@ -246,7 +246,7 @@ int pp = p-1; int iter = 0; - double eps = pow(2.0,-52.0); + Real eps(1e-10); //pow(2.0,-52.0)); while (p > 0) { int k=0; int kase=0; @@ -280,8 +280,8 @@ if (ks == k) { break; } - double t = (ks != p ? abs(e[ks]) : 0.) + - (ks != k+1 ? abs(e[ks-1]) : 0.); + Real t( (ks != p ? abs(e[ks]) : 0.) + + (ks != k+1 ? abs(e[ks-1]) : 0.)); if (abs(s[ks]) <= eps*t) { s[ks] = 0.0; break; @@ -305,12 +305,12 @@ // Deflate negligible s(p). case 1: { - double f = e[p-2]; + Real f(e[p-2]); e[p-2] = 0.0; for (j = p-2; j >= k; j--) { - double t = hypot(s[j],f); - double cs = s[j]/t; - double sn = f/t; + Real t( hypot(s[j],f)); + Real cs(s[j]/t); + Real sn(f/t); s[j] = t; if (j != k) { f = -sn*e[j-1]; @@ -330,12 +330,12 @@ // Split at negligible s(k). case 2: { - double f = e[k-1]; + Real f(e[k-1]); e[k-1] = 0.0; for (j = k; j < p; j++) { - double t = hypot(s[j],f); - double cs = s[j]/t; - double sn = f/t; + Real t(hypot(s[j],f)); + Real cs( s[j]/t); + Real sn(f/t); s[j] = t; f = -sn*e[j]; e[j] = cs*e[j]; @@ -356,17 +356,17 @@ // Calculate the shift. - double scale = max(max(max(max( + Real scale = max(max(max(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; - double spm1 = s[p-2]/scale; - double epm1 = e[p-2]/scale; - double sk = s[k]/scale; - double ek = e[k]/scale; - double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0; - double c = (sp*epm1)*(sp*epm1); - double shift = 0.0; + Real sp = s[p-1]/scale; + Real spm1 = s[p-2]/scale; + Real epm1 = e[p-2]/scale; + Real sk = s[k]/scale; + Real ek = e[k]/scale; + Real b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0; + Real c = (sp*epm1)*(sp*epm1); + Real shift = 0.0; if ((b != 0.0) || (c != 0.0)) { shift = sqrt(b*b + c); if (b < 0.0) { @@ -374,15 +374,15 @@ } shift = c/(b + shift); } - double f = (sk + sp)*(sk - sp) + shift; - double g = sk*ek; + Real f = (sk + sp)*(sk - sp) + shift; + Real g = sk*ek; // Chase zeros. for (j = k; j < p-1; j++) { - double t = hypot(f,g); - double cs = f/t; - double sn = g/t; + Real t = hypot(f,g); + Real cs = f/t; + Real sn = g/t; if (j != k) { e[j-1] = t; } @@ -439,7 +439,7 @@ if (s[k] >= s[k+1]) { break; } - double t = s[k]; + Real t = s[k]; s[k] = s[k+1]; s[k+1] = t; if (wantv && (k < n-1)) { @@ -505,13 +505,13 @@ /** Two norm (max(S)) */ - double norm2 () { + Real norm2 () { return s[0]; } /** Two norm of condition number (max(S)/min(S)) */ - double cond () { + Real cond () { return s[0]/s[min(m,n)-1]; } @@ -521,8 +521,8 @@ int rank () { - double eps = pow(2.0,-52.0); - double tol = max(m,n)*s[0]*eps; + Real eps = 1e-10; //pow(2.0,-52.0); + Real tol = max(m,n)*s[0]*eps; int r = 0; for (int i = 0; i < s.dim(); i++) { if (s[i] > tol) { Index: tnt_array1d.h =================================================================== --- tnt_array1d.h (revision 20) +++ tnt_array1d.h (revision 49) @@ -235,7 +235,7 @@ template inline Array1D Array1D::subarray(int i0, int i1) { - if ((i0 > 0) && (i1 < n_) || (i0 <= i1)) + if (((i0 > 0) && (i1 < n_)) || (i0 <= i1)) { Array1D X(*this); /* create a new instance of this array. */ X.n_ = i1-i0+1;