From a0c8925b76d8d3fbd24dc4b9b3c84a8f1b823d4b Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 17 Jul 2014 09:29:35 +0000 Subject: Remove tnt/jama This was used in the past by SCM/PDM, but hasn't been used for a long time now. --- lib/include/tnt/jama_cholesky.h | 258 ------- lib/include/tnt/jama_eig.h | 1034 --------------------------- lib/include/tnt/jama_lu.h | 323 --------- lib/include/tnt/jama_qr.h | 326 --------- lib/include/tnt/jama_svd.h | 543 -------------- lib/include/tnt/tnt.h | 64 -- lib/include/tnt/tnt_array1d.h | 278 ------- lib/include/tnt/tnt_array1d_utils.h | 230 ------ lib/include/tnt/tnt_array2d.h | 315 -------- lib/include/tnt/tnt_array2d_utils.h | 287 -------- lib/include/tnt/tnt_array3d.h | 296 -------- lib/include/tnt/tnt_array3d_utils.h | 236 ------ lib/include/tnt/tnt_cmat.h | 580 --------------- lib/include/tnt/tnt_fortran_array1d.h | 267 ------- lib/include/tnt/tnt_fortran_array1d_utils.h | 242 ------- lib/include/tnt/tnt_fortran_array2d.h | 225 ------ lib/include/tnt/tnt_fortran_array2d_utils.h | 236 ------ lib/include/tnt/tnt_fortran_array3d.h | 223 ------ lib/include/tnt/tnt_fortran_array3d_utils.h | 249 ------- lib/include/tnt/tnt_i_refvec.h | 243 ------- lib/include/tnt/tnt_math_utils.h | 34 - lib/include/tnt/tnt_sparse_matrix_csr.h | 103 --- lib/include/tnt/tnt_stopwatch.h | 95 --- lib/include/tnt/tnt_subscript.h | 54 -- lib/include/tnt/tnt_vec.h | 404 ----------- lib/include/tnt/tnt_version.h | 39 - lib/licenses/tnt.txt | 14 - 27 files changed, 7198 deletions(-) delete mode 100644 lib/include/tnt/jama_cholesky.h delete mode 100644 lib/include/tnt/jama_eig.h delete mode 100644 lib/include/tnt/jama_lu.h delete mode 100644 lib/include/tnt/jama_qr.h delete mode 100644 lib/include/tnt/jama_svd.h delete mode 100644 lib/include/tnt/tnt.h delete mode 100644 lib/include/tnt/tnt_array1d.h delete mode 100644 lib/include/tnt/tnt_array1d_utils.h delete mode 100644 lib/include/tnt/tnt_array2d.h delete mode 100644 lib/include/tnt/tnt_array2d_utils.h delete mode 100644 lib/include/tnt/tnt_array3d.h delete mode 100644 lib/include/tnt/tnt_array3d_utils.h delete mode 100644 lib/include/tnt/tnt_cmat.h delete mode 100644 lib/include/tnt/tnt_fortran_array1d.h delete mode 100644 lib/include/tnt/tnt_fortran_array1d_utils.h delete mode 100644 lib/include/tnt/tnt_fortran_array2d.h delete mode 100644 lib/include/tnt/tnt_fortran_array2d_utils.h delete mode 100644 lib/include/tnt/tnt_fortran_array3d.h delete mode 100644 lib/include/tnt/tnt_fortran_array3d_utils.h delete mode 100644 lib/include/tnt/tnt_i_refvec.h delete mode 100644 lib/include/tnt/tnt_math_utils.h delete mode 100644 lib/include/tnt/tnt_sparse_matrix_csr.h delete mode 100644 lib/include/tnt/tnt_stopwatch.h delete mode 100644 lib/include/tnt/tnt_subscript.h delete mode 100644 lib/include/tnt/tnt_vec.h delete mode 100644 lib/include/tnt/tnt_version.h delete mode 100644 lib/licenses/tnt.txt (limited to 'lib') diff --git a/lib/include/tnt/jama_cholesky.h b/lib/include/tnt/jama_cholesky.h deleted file mode 100644 index 371d2f7..0000000 --- a/lib/include/tnt/jama_cholesky.h +++ /dev/null @@ -1,258 +0,0 @@ -#ifndef JAMA_CHOLESKY_H -#define JAMA_CHOLESKY_H - -#include "math.h" - /* needed for sqrt() below. */ - - -namespace JAMA -{ - -using namespace TNT; - -/** -

- For a symmetric, positive definite matrix A, this function - computes the Cholesky factorization, i.e. it computes a lower - triangular matrix L such that A = L*L'. - If the matrix is not symmetric or positive definite, the function - computes only a partial decomposition. This can be tested with - the is_spd() flag. - -

Typical usage looks like: - 

-	Array2D A(n,n);
-	Array2D L;
-
-	 ... 
-
-	Cholesky chol(A);
-
-	if (chol.is_spd())
-		L = chol.getL();
-		
-  	else
-		cout << "factorization was not complete.\n";
-
-
- - -

- (Adapted from JAMA, a Java Matrix Library, developed by jointly - by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). - - */ - -template -class Cholesky -{ - 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); - int is_spd() const; - -}; - -template -Cholesky::Cholesky() : L_(0,0), isspd(0) {} - -/** - @return 1, if original matrix to be factored was symmetric - positive-definite (SPD). -*/ -template -int Cholesky::is_spd() const -{ - return isspd; -} - -/** - @return the lower triangular factor, L, such that L*L'=A. -*/ -template -Array2D Cholesky::getL() const -{ - return L_; -} - -/** - Constructs a lower triangular matrix L, such that L*L'= A. - If A is not symmetric positive-definite (SPD), only a - partial factorization is performed. If is_spd() - evalutate true (1) then the factorizaiton was successful. -*/ -template -Cholesky::Cholesky(const Array2D &A) -{ - - - int m = A.dim1(); - int n = A.dim2(); - - isspd = (m == n); - - if (m != n) - { - L_ = Array2D(0,0); - return; - } - - L_ = Array2D(n,n); - - - // Main loop. - for (int j = 0; j < n; j++) - { - Real d(0.0); - for (int k = 0; k < j; k++) - { - Real s(0.0); - for (int i = 0; i < k; i++) - { - s += L_[k][i]*L_[j][i]; - } - L_[j][k] = s = (A[j][k] - s)/L_[k][k]; - d = d + s*s; - isspd = isspd && (A[k][j] == A[j][k]); - } - d = A[j][j] - d; - isspd = isspd && (d > 0.0); - L_[j][j] = sqrt(d > 0.0 ? d : 0.0); - for (int k = j+1; k < n; k++) - { - L_[j][k] = 0.0; - } - } -} - -/** - - Solve a linear system A*x = b, using the previously computed - cholesky factorization of A: L*L'. - - @param B A Matrix with as many rows as A and any number of columns. - @return x so that L*L'*x = b. If b is nonconformat, or if A - was not symmetric posidtive definite, a null (0x0) - array is returned. -*/ -template -Array1D Cholesky::solve(const Array1D &b) -{ - int n = L_.dim1(); - if (b.dim1() != n) - return Array1D(); - - - Array1D x = b.copy(); - - - // Solve L*y = b; - for (int k = 0; k < n; k++) - { - for (int i = 0; i < k; i++) - x[k] -= x[i]*L_[k][i]; - x[k] /= L_[k][k]; - - } - - // Solve L'*X = Y; - for (int k = n-1; k >= 0; k--) - { - for (int i = k+1; i < n; i++) - x[k] -= x[i]*L_[i][k]; - x[k] /= L_[k][k]; - } - - return x; -} - - -/** - - Solve a linear system A*X = B, using the previously computed - cholesky factorization of A: L*L'. - - @param B A Matrix with as many rows as A and any number of columns. - @return X so that L*L'*X = B. If B is nonconformat, or if A - was not symmetric posidtive definite, a null (0x0) - array is returned. -*/ -template -Array2D Cholesky::solve(const Array2D &B) -{ - int n = L_.dim1(); - if (B.dim1() != n) - return Array2D(); - - - Array2D X = B.copy(); - int nx = B.dim2(); - -// Cleve's original code -#if 0 - // Solve L*Y = B; - for (int k = 0; k < n; k++) { - for (int i = k+1; i < n; i++) { - for (int j = 0; j < nx; j++) { - X[i][j] -= X[k][j]*L_[k][i]; - } - } - for (int j = 0; j < nx; j++) { - X[k][j] /= L_[k][k]; - } - } - - // Solve L'*X = Y; - for (int k = n-1; k >= 0; k--) { - for (int j = 0; j < nx; j++) { - X[k][j] /= L_[k][k]; - } - for (int i = 0; i < k; i++) { - for (int j = 0; j < nx; j++) { - X[i][j] -= X[k][j]*L_[k][i]; - } - } - } -#endif - - - // Solve L*y = b; - for (int j=0; j< nx; j++) - { - for (int k = 0; k < n; k++) - { - for (int i = 0; i < k; i++) - X[k][j] -= X[i][j]*L_[k][i]; - X[k][j] /= L_[k][k]; - } - } - - // Solve L'*X = Y; - for (int j=0; j= 0; k--) - { - for (int i = k+1; i < n; i++) - X[k][j] -= X[i][j]*L_[i][k]; - X[k][j] /= L_[k][k]; - } - } - - - - return X; -} - - -} -// namespace JAMA - -#endif -// JAMA_CHOLESKY_H diff --git a/lib/include/tnt/jama_eig.h b/lib/include/tnt/jama_eig.h deleted file mode 100644 index 8e3572f..0000000 --- a/lib/include/tnt/jama_eig.h +++ /dev/null @@ -1,1034 +0,0 @@ -#ifndef JAMA_EIG_H -#define JAMA_EIG_H - - -#include "tnt_array1d.h" -#include "tnt_array2d.h" -#include "tnt_math_utils.h" - -#include -// for min(), max() below - -#include -// for abs() below - -using namespace TNT; -using namespace std; - -// Modification by Willem Jan Palenstijn, 2010-03-11: -// Use std::min() instead of min(), std::max() instead of max() - - -namespace JAMA -{ - -/** - - Computes eigenvalues and eigenvectors of a real (non-complex) - matrix. -

- If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is - diagonal and the eigenvector matrix V is orthogonal. That is, - the diagonal values of D are the eigenvalues, and - V*V' = I, where I is the identity matrix. The columns of V - represent the eigenvectors in the sense that A*V = V*D. - -

- If A is not symmetric, then the eigenvalue matrix D is block diagonal - with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, - a + i*b, in 2-by-2 blocks, [a, b; -b, a]. That is, if the complex - eigenvalues look like -

-
-          u + iv     .        .          .      .    .
-            .      u - iv     .          .      .    .
-            .        .      a + ib       .      .    .
-            .        .        .        a - ib   .    .
-            .        .        .          .      x    .
-            .        .        .          .      .    y
-
- then D looks like -
-
-            u        v        .          .      .    .
-           -v        u        .          .      .    . 
-            .        .        a          b      .    .
-            .        .       -b          a      .    .
-            .        .        .          .      x    .
-            .        .        .          .      .    y
-
- This keeps V a real matrix in both symmetric and non-symmetric - cases, and A*V = V*D. - - - -

- The matrix V may be badly - conditioned, or even singular, so the validity of the equation - A = V*D*inverse(V) depends upon the condition number of V. - -

- (Adapted from JAMA, a Java Matrix Library, developed by jointly - by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). -**/ - -template -class Eigenvalue -{ - - - /** Row and column dimension (square matrix). */ - int n; - - int issymmetric; /* boolean*/ - - /** Arrays for internal storage of eigenvalues. */ - - TNT::Array1D d; /* real part */ - TNT::Array1D e; /* img part */ - - /** Array for internal storage of eigenvectors. */ - TNT::Array2D V; - - /** Array for internal storage of nonsymmetric Hessenberg form. - @serial internal storage of nonsymmetric Hessenberg form. - */ - TNT::Array2D H; - - - /** Working storage for nonsymmetric algorithm. - @serial working storage for nonsymmetric algorithm. - */ - TNT::Array1D ort; - - - // Symmetric Householder reduction to tridiagonal form. - - void tred2() { - - // This is derived from the Algol procedures tred2 by - // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for - // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutine in EISPACK. - - for (int j = 0; j < n; j++) { - d[j] = V[n-1][j]; - } - - // Householder reduction to tridiagonal form. - - for (int i = n-1; i > 0; i--) { - - // Scale to avoid under/overflow. - - Real scale = 0.0; - Real h = 0.0; - for (int k = 0; k < i; k++) { - scale = scale + abs(d[k]); - } - if (scale == 0.0) { - e[i] = d[i-1]; - for (int j = 0; j < i; j++) { - d[j] = V[i-1][j]; - V[i][j] = 0.0; - V[j][i] = 0.0; - } - } else { - - // Generate Householder vector. - - for (int k = 0; k < i; k++) { - d[k] /= scale; - h += d[k] * d[k]; - } - Real f = d[i-1]; - Real g = sqrt(h); - if (f > 0) { - g = -g; - } - e[i] = scale * g; - h = h - f * g; - d[i-1] = f - g; - for (int j = 0; j < i; j++) { - e[j] = 0.0; - } - - // Apply similarity transformation to remaining columns. - - for (int j = 0; j < i; j++) { - f = d[j]; - V[j][i] = f; - g = e[j] + V[j][j] * f; - for (int k = j+1; k <= i-1; k++) { - g += V[k][j] * d[k]; - e[k] += V[k][j] * f; - } - e[j] = g; - } - f = 0.0; - for (int j = 0; j < i; j++) { - e[j] /= h; - f += e[j] * d[j]; - } - Real hh = f / (h + h); - for (int j = 0; j < i; j++) { - e[j] -= hh * d[j]; - } - for (int j = 0; j < i; j++) { - f = d[j]; - g = e[j]; - for (int k = j; k <= i-1; k++) { - V[k][j] -= (f * e[k] + g * d[k]); - } - d[j] = V[i-1][j]; - V[i][j] = 0.0; - } - } - d[i] = h; - } - - // Accumulate transformations. - - for (int i = 0; i < n-1; i++) { - V[n-1][i] = V[i][i]; - V[i][i] = 1.0; - Real h = d[i+1]; - if (h != 0.0) { - for (int k = 0; k <= i; k++) { - d[k] = V[k][i+1] / h; - } - for (int j = 0; j <= i; j++) { - Real g = 0.0; - for (int k = 0; k <= i; k++) { - g += V[k][i+1] * V[k][j]; - } - for (int k = 0; k <= i; k++) { - V[k][j] -= g * d[k]; - } - } - } - for (int k = 0; k <= i; k++) { - V[k][i+1] = 0.0; - } - } - for (int j = 0; j < n; j++) { - d[j] = V[n-1][j]; - V[n-1][j] = 0.0; - } - V[n-1][n-1] = 1.0; - e[0] = 0.0; - } - - // Symmetric tridiagonal QL algorithm. - - void tql2 () { - - // This is derived from the Algol procedures tql2, by - // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for - // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutine in EISPACK. - - for (int i = 1; i < n; i++) { - e[i-1] = e[i]; - } - e[n-1] = 0.0; - - Real f = 0.0; - Real tst1 = 0.0; - Real eps = pow(2.0,-52.0); - for (int l = 0; l < n; l++) { - - // Find small subdiagonal element - - tst1 = std::max(tst1,abs(d[l]) + abs(e[l])); - int m = l; - - // Original while-loop from Java code - while (m < n) { - if (abs(e[m]) <= eps*tst1) { - break; - } - m++; - } - - - // If m == l, d[l] is an eigenvalue, - // otherwise, iterate. - - if (m > l) { - int iter = 0; - do { - iter = iter + 1; // (Could check iteration count here.) - - // Compute implicit shift - - Real g = d[l]; - Real p = (d[l+1] - g) / (2.0 * e[l]); - Real r = hypot(p,1.0); - if (p < 0) { - r = -r; - } - d[l] = e[l] / (p + r); - d[l+1] = e[l] * (p + r); - Real dl1 = d[l+1]; - Real h = g - d[l]; - for (int i = l+2; i < n; i++) { - d[i] -= h; - } - f = f + h; - - // Implicit QL transformation. - - p = d[m]; - Real c = 1.0; - Real c2 = c; - Real c3 = c; - Real el1 = e[l+1]; - Real s = 0.0; - Real s2 = 0.0; - for (int i = m-1; i >= l; i--) { - c3 = c2; - c2 = c; - s2 = s; - g = c * e[i]; - h = c * p; - r = hypot(p,e[i]); - e[i+1] = s * r; - s = e[i] / r; - c = p / r; - p = c * d[i] - s * g; - d[i+1] = h + s * (c * g + s * d[i]); - - // Accumulate transformation. - - for (int k = 0; k < n; k++) { - h = V[k][i+1]; - V[k][i+1] = s * V[k][i] + c * h; - V[k][i] = c * V[k][i] - s * h; - } - } - p = -s * s2 * c3 * el1 * e[l] / dl1; - e[l] = s * p; - d[l] = c * p; - - // Check for convergence. - - } while (abs(e[l]) > eps*tst1); - } - d[l] = d[l] + f; - e[l] = 0.0; - } - - // Sort eigenvalues and corresponding vectors. - - for (int i = 0; i < n-1; i++) { - int k = i; - Real p = d[i]; - for (int j = i+1; j < n; j++) { - if (d[j] < p) { - k = j; - p = d[j]; - } - } - if (k != i) { - d[k] = d[i]; - d[i] = p; - for (int j = 0; j < n; j++) { - p = V[j][i]; - V[j][i] = V[j][k]; - V[j][k] = p; - } - } - } - } - - // Nonsymmetric reduction to Hessenberg form. - - void orthes () { - - // This is derived from the Algol procedures orthes and ortran, - // by Martin and Wilkinson, Handbook for Auto. Comp., - // Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutines in EISPACK. - - int low = 0; - int high = n-1; - - for (int m = low+1; m <= high-1; m++) { - - // Scale column. - - Real scale = 0.0; - for (int i = m; i <= high; i++) { - scale = scale + abs(H[i][m-1]); - } - if (scale != 0.0) { - - // Compute Householder transformation. - - Real h = 0.0; - for (int i = high; i >= m; i--) { - ort[i] = H[i][m-1]/scale; - h += ort[i] * ort[i]; - } - Real g = sqrt(h); - if (ort[m] > 0) { - g = -g; - } - h = h - ort[m] * g; - ort[m] = ort[m] - g; - - // Apply Householder similarity transformation - // H = (I-u*u'/h)*H*(I-u*u')/h) - - for (int j = m; j < n; j++) { - Real f = 0.0; - for (int i = high; i >= m; i--) { - f += ort[i]*H[i][j]; - } - f = f/h; - for (int i = m; i <= high; i++) { - H[i][j] -= f*ort[i]; - } - } - - for (int i = 0; i <= high; i++) { - Real f = 0.0; - for (int j = high; j >= m; j--) { - f += ort[j]*H[i][j]; - } - f = f/h; - for (int j = m; j <= high; j++) { - H[i][j] -= f*ort[j]; - } - } - ort[m] = scale*ort[m]; - H[m][m-1] = scale*g; - } - } - - // Accumulate transformations (Algol's ortran). - - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - V[i][j] = (i == j ? 1.0 : 0.0); - } - } - - for (int m = high-1; m >= low+1; m--) { - if (H[m][m-1] != 0.0) { - for (int i = m+1; i <= high; i++) { - ort[i] = H[i][m-1]; - } - for (int j = m; j <= high; j++) { - Real g = 0.0; - for (int i = m; i <= high; i++) { - g += ort[i] * V[i][j]; - } - // Double division avoids possible underflow - g = (g / ort[m]) / H[m][m-1]; - for (int i = m; i <= high; i++) { - V[i][j] += g * ort[i]; - } - } - } - } - } - - - // Complex scalar division. - - Real cdivr, cdivi; - void cdiv(Real xr, Real xi, Real yr, Real yi) { - Real r,d; - if (abs(yr) > abs(yi)) { - r = yi/yr; - d = yr + r*yi; - cdivr = (xr + r*xi)/d; - cdivi = (xi - r*xr)/d; - } else { - r = yr/yi; - d = yi + r*yr; - cdivr = (r*xr + xi)/d; - cdivi = (r*xi - xr)/d; - } - } - - - // Nonsymmetric reduction from Hessenberg to real Schur form. - - void hqr2 () { - - // This is derived from the Algol procedure hqr2, - // by Martin and Wilkinson, Handbook for Auto. Comp., - // Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutine in EISPACK. - - // Initialize - - int nn = this->n; - int n = nn-1; - int low = 0; - int high = nn-1; - Real eps = pow(2.0,-52.0); - Real exshift = 0.0; - Real p=0,q=0,r=0,s=0,z=0,t,w,x,y; - - // Store roots isolated by balanc and compute matrix norm - - Real norm = 0.0; - for (int i = 0; i < nn; i++) { - if ((i < low) || (i > high)) { - d[i] = H[i][i]; - e[i] = 0.0; - } - for (int j = std::max(i-1,0); j < nn; j++) { - norm = norm + abs(H[i][j]); - } - } - - // Outer loop over eigenvalue index - - int iter = 0; - while (n >= low) { - - // Look for single small sub-diagonal element - - int l = n; - while (l > low) { - s = abs(H[l-1][l-1]) + abs(H[l][l]); - if (s == 0.0) { - s = norm; - } - if (abs(H[l][l-1]) < eps * s) { - break; - } - l--; - } - - // Check for convergence - // One root found - - if (l == n) { - H[n][n] = H[n][n] + exshift; - d[n] = H[n][n]; - e[n] = 0.0; - n--; - iter = 0; - - // Two roots found - - } else if (l == n-1) { - 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)); - H[n][n] = H[n][n] + exshift; - H[n-1][n-1] = H[n-1][n-1] + exshift; - x = H[n][n]; - - // Real pair - - if (q >= 0) { - if (p >= 0) { - z = p + z; - } else { - z = p - z; - } - d[n-1] = x + z; - d[n] = d[n-1]; - if (z != 0.0) { - d[n] = x - w / z; - } - e[n-1] = 0.0; - e[n] = 0.0; - x = H[n][n-1]; - s = abs(x) + abs(z); - p = x / s; - q = z / s; - r = sqrt(p * p+q * q); - p = p / r; - q = q / r; - - // Row modification - - for (int j = n-1; j < nn; j++) { - z = H[n-1][j]; - H[n-1][j] = q * z + p * H[n][j]; - H[n][j] = q * H[n][j] - p * z; - } - - // Column modification - - for (int i = 0; i <= n; i++) { - z = H[i][n-1]; - H[i][n-1] = q * z + p * H[i][n]; - H[i][n] = q * H[i][n] - p * z; - } - - // Accumulate transformations - - for (int i = low; i <= high; i++) { - z = V[i][n-1]; - V[i][n-1] = q * z + p * V[i][n]; - V[i][n] = q * V[i][n] - p * z; - } - - // Complex pair - - } else { - d[n-1] = x + p; - d[n] = x + p; - e[n-1] = z; - e[n] = -z; - } - n = n - 2; - iter = 0; - - // No convergence yet - - } else { - - // Form shift - - x = H[n][n]; - y = 0.0; - w = 0.0; - if (l < n) { - y = H[n-1][n-1]; - w = H[n][n-1] * H[n-1][n]; - } - - // Wilkinson's original ad hoc shift - - if (iter == 10) { - exshift += x; - for (int i = low; i <= n; i++) { - H[i][i] -= x; - } - s = abs(H[n][n-1]) + abs(H[n-1][n-2]); - x = y = 0.75 * s; - w = -0.4375 * s * s; - } - - // MATLAB's new ad hoc shift - - if (iter == 30) { - s = (y - x) / 2.0; - s = s * s + w; - if (s > 0) { - s = sqrt(s); - if (y < x) { - s = -s; - } - s = x - w / ((y - x) / 2.0 + s); - for (int i = low; i <= n; i++) { - H[i][i] -= s; - } - exshift += s; - x = y = w = 0.964; - } - } - - iter = iter + 1; // (Could check iteration count here.) - - // Look for two consecutive small sub-diagonal elements - - int m = n-2; - while (m >= l) { - z = H[m][m]; - r = x - z; - s = y - z; - 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); - 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])))) { - break; - } - m--; - } - - for (int i = m+2; i <= n; i++) { - H[i][i-2] = 0.0; - if (i > m+2) { - H[i][i-3] = 0.0; - } - } - - // Double QR step involving rows l:n and columns m:n - - for (int k = m; k <= n-1; k++) { - int notlast = (k != n-1); - if (k != m) { - 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); - if (x != 0.0) { - p = p / x; - q = q / x; - r = r / x; - } - } - if (x == 0.0) { - break; - } - s = sqrt(p * p + q * q + r * r); - if (p < 0) { - s = -s; - } - if (s != 0) { - if (k != m) { - H[k][k-1] = -s * x; - } else if (l != m) { - H[k][k-1] = -H[k][k-1]; - } - p = p + s; - x = p / s; - y = q / s; - z = r / s; - q = q / p; - r = r / p; - - // Row modification - - for (int j = k; j < nn; j++) { - p = H[k][j] + q * H[k+1][j]; - if (notlast) { - p = p + r * H[k+2][j]; - H[k+2][j] = H[k+2][j] - p * z; - } - H[k][j] = H[k][j] - p * x; - H[k+1][j] = H[k+1][j] - p * y; - } - - // Column modification - - 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]; - H[i][k+2] = H[i][k+2] - p * r; - } - H[i][k] = H[i][k] - p; - H[i][k+1] = H[i][k+1] - p * q; - } - - // Accumulate transformations - - for (int i = low; i <= high; i++) { - p = x * V[i][k] + y * V[i][k+1]; - if (notlast) { - p = p + z * V[i][k+2]; - V[i][k+2] = V[i][k+2] - p * r; - } - V[i][k] = V[i][k] - p; - V[i][k+1] = V[i][k+1] - p * q; - } - } // (s != 0) - } // k loop - } // check convergence - } // while (n >= low) - - // Backsubstitute to find vectors of upper triangular form - - if (norm == 0.0) { - return; - } - - for (n = nn-1; n >= 0; n--) { - p = d[n]; - q = e[n]; - - // Real vector - - if (q == 0) { - int l = n; - H[n][n] = 1.0; - for (int i = n-1; i >= 0; i--) { - w = H[i][i] - p; - r = 0.0; - for (int j = l; j <= n; j++) { - r = r + H[i][j] * H[j][n]; - } - if (e[i] < 0.0) { - z = w; - s = r; - } else { - l = i; - if (e[i] == 0.0) { - if (w != 0.0) { - H[i][n] = -r / w; - } else { - H[i][n] = -r / (eps * norm); - } - - // Solve real equations - - } else { - x = H[i][i+1]; - y = H[i+1][i]; - 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)) { - H[i+1][n] = (-r - w * t) / x; - } else { - H[i+1][n] = (-s - y * t) / z; - } - } - - // Overflow control - - t = abs(H[i][n]); - if ((eps * t) * t > 1) { - for (int j = i; j <= n; j++) { - H[j][n] = H[j][n] / t; - } - } - } - } - - // Complex vector - - } else if (q < 0) { - int l = n-1; - - // Last vector component imaginary so matrix is triangular - - if (abs(H[n][n-1]) > 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 { - cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q); - H[n-1][n-1] = cdivr; - H[n-1][n] = cdivi; - } - H[n][n-1] = 0.0; - H[n][n] = 1.0; - for (int i = n-2; i >= 0; i--) { - Real ra,sa,vr,vi; - ra = 0.0; - sa = 0.0; - for (int j = l; j <= n; j++) { - ra = ra + H[i][j] * H[j][n-1]; - sa = sa + H[i][j] * H[j][n]; - } - w = H[i][i] - p; - - if (e[i] < 0.0) { - z = w; - r = ra; - s = sa; - } else { - l = i; - if (e[i] == 0) { - cdiv(-ra,-sa,w,q); - H[i][n-1] = cdivr; - H[i][n] = cdivi; - } else { - - // Solve complex equations - - x = H[i][i+1]; - y = H[i+1][i]; - 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)); - } - 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))) { - 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 { - cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q); - H[i+1][n-1] = cdivr; - H[i+1][n] = cdivi; - } - } - - // Overflow control - - 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; - H[j][n] = H[j][n] / t; - } - } - } - } - } - } - - // Vectors of isolated roots - - for (int i = 0; i < nn; i++) { - if (i < low || i > high) { - for (int j = i; j < nn; j++) { - V[i][j] = H[i][j]; - } - } - } - - // Back transformation to get eigenvectors of original matrix - - for (int j = nn-1; j >= low; j--) { - for (int i = low; i <= high; i++) { - z = 0.0; - for (int k = low; k <= std::min(j,high); k++) { - z = z + V[i][k] * H[k][j]; - } - V[i][j] = z; - } - } - } - -public: - - - /** Check for symmetry, then construct the eigenvalue decomposition - @param A Square real (non-complex) matrix - */ - - Eigenvalue(const TNT::Array2D &A) { - n = A.dim2(); - V = Array2D(n,n); - d = Array1D(n); - e = Array1D(n); - - issymmetric = 1; - for (int j = 0; (j < n) && issymmetric; j++) { - for (int i = 0; (i < n) && issymmetric; i++) { - issymmetric = (A[i][j] == A[j][i]); - } - } - - if (issymmetric) { - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - V[i][j] = A[i][j]; - } - } - - // Tridiagonalize. - tred2(); - - // Diagonalize. - tql2(); - - } else { - H = TNT::Array2D(n,n); - ort = TNT::Array1D(n); - - for (int j = 0; j < n; j++) { - for (int i = 0; i < n; i++) { - H[i][j] = A[i][j]; - } - } - - // Reduce to Hessenberg form. - orthes(); - - // Reduce Hessenberg to real Schur form. - hqr2(); - } - } - - - /** Return the eigenvector matrix - @return V - */ - - void getV (TNT::Array2D &V_) { - V_ = V; - return; - } - - /** Return the real parts of the eigenvalues - @return real(diag(D)) - */ - - void getRealEigenvalues (TNT::Array1D &d_) { - d_ = d; - return ; - } - - /** Return the imaginary parts of the eigenvalues - in parameter e_. - - @pararm e_: new matrix with imaginary parts of the eigenvalues. - */ - void getImagEigenvalues (TNT::Array1D &e_) { - e_ = e; - return; - } - - -/** - Computes the block diagonal eigenvalue matrix. - If the original matrix A is not symmetric, then the eigenvalue - matrix D is block diagonal with the real eigenvalues in 1-by-1 - blocks and any complex eigenvalues, - a + i*b, in 2-by-2 blocks, [a, b; -b, a]. That is, if the complex - eigenvalues look like -

-
-          u + iv     .        .          .      .    .
-            .      u - iv     .          .      .    .
-            .        .      a + ib       .      .    .
-            .        .        .        a - ib   .    .
-            .        .        .          .      x    .
-            .        .        .          .      .    y
-
- then D looks like -
-
-            u        v        .          .      .    .
-           -v        u        .          .      .    . 
-            .        .        a          b      .    .
-            .        .       -b          a      .    .
-            .        .        .          .      x    .
-            .        .        .          .      .    y
-
- This keeps V a real matrix in both symmetric and non-symmetric - cases, and A*V = V*D. - - @param D: upon return, the matrix is filled with the block diagonal - eigenvalue matrix. - -*/ - void getD (TNT::Array2D &D) { - D = Array2D(n,n); - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - D[i][j] = 0.0; - } - D[i][i] = d[i]; - if (e[i] > 0) { - D[i][i+1] = e[i]; - } else if (e[i] < 0) { - D[i][i-1] = e[i]; - } - } - } -}; - -} //namespace JAMA - - -#endif -// JAMA_EIG_H diff --git a/lib/include/tnt/jama_lu.h b/lib/include/tnt/jama_lu.h deleted file mode 100644 index e95b433..0000000 --- a/lib/include/tnt/jama_lu.h +++ /dev/null @@ -1,323 +0,0 @@ -#ifndef JAMA_LU_H -#define JAMA_LU_H - -#include "tnt.h" -#include -//for min(), max() below - -using namespace TNT; -using namespace std; - - -// Modification by Willem Jan Palenstijn, 2010-03-11: -// Use std::min() instead of min() - -namespace JAMA -{ - - /** LU Decomposition. -

- For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n - unit lower triangular matrix L, an n-by-n upper triangular matrix U, - and a permutation vector piv of length m so that A(piv,:) = L*U. - If m < n, then L is m-by-m and U is m-by-n. -

- The LU decompostion with pivoting always exists, even if the matrix is - singular, so the constructor will never fail. The primary use of the - LU decomposition is in the solution of square systems of simultaneous - linear equations. This will fail if isNonsingular() returns false. - */ -template -class LU -{ - - - - /* Array for internal storage of decomposition. */ - Array2D LU_; - int m, n, pivsign; - Array1D piv; - - - Array2D permute_copy(const Array2D &A, - const Array1D &piv, int j0, int j1) - { - int piv_length = piv.dim(); - - Array2D X(piv_length, j1-j0+1); - - - for (int i = 0; i < piv_length; i++) - for (int j = j0; j <= j1; j++) - X[i][j-j0] = A[piv[i]][j]; - - return X; - } - - Array1D permute_copy(const Array1D &A, - const Array1D &piv) - { - int piv_length = piv.dim(); - if (piv_length != A.dim()) - return Array1D(); - - Array1D x(piv_length); - - - for (int i = 0; i < piv_length; i++) - x[i] = A[piv[i]]; - - return x; - } - - - public : - - /** LU Decomposition - @param A Rectangular matrix - @return LU Decomposition object to access L, U and piv. - */ - - LU (const Array2D &A) : LU_(A.copy()), m(A.dim1()), n(A.dim2()), - piv(A.dim1()) - - { - - // Use a "left-looking", dot-product, Crout/Doolittle algorithm. - - - for (int i = 0; i < m; i++) { - piv[i] = i; - } - pivsign = 1; - Real *LUrowi = 0;; - Array1D LUcolj(m); - - // Outer loop. - - for (int j = 0; j < n; j++) { - - // Make a copy of the j-th column to localize references. - - for (int i = 0; i < m; i++) { - LUcolj[i] = LU_[i][j]; - } - - // Apply previous transformations. - - for (int i = 0; i < m; i++) { - LUrowi = LU_[i]; - - // Most of the time is spent in the following dot product. - - int kmax = std::min(i,j); - double s = 0.0; - for (int k = 0; k < kmax; k++) { - s += LUrowi[k]*LUcolj[k]; - } - - LUrowi[j] = LUcolj[i] -= s; - } - - // Find pivot and exchange if necessary. - - int p = j; - for (int i = j+1; i < m; i++) { - if (abs(LUcolj[i]) > abs(LUcolj[p])) { - p = i; - } - } - if (p != j) { - int k=0; - for (k = 0; k < n; k++) { - double t = LU_[p][k]; - LU_[p][k] = LU_[j][k]; - LU_[j][k] = t; - } - k = piv[p]; - piv[p] = piv[j]; - piv[j] = k; - pivsign = -pivsign; - } - - // Compute multipliers. - - if ((j < m) && (LU_[j][j] != 0.0)) { - for (int i = j+1; i < m; i++) { - LU_[i][j] /= LU_[j][j]; - } - } - } - } - - - /** Is the matrix nonsingular? - @return 1 (true) if upper triangular factor U (and hence A) - is nonsingular, 0 otherwise. - */ - - int isNonsingular () { - for (int j = 0; j < n; j++) { - if (LU_[j][j] == 0) - return 0; - } - return 1; - } - - /** Return lower triangular factor - @return L - */ - - Array2D getL () { - Array2D L_(m,n); - for (int i = 0; i < m; i++) { - for (int j = 0; j < n; j++) { - if (i > j) { - L_[i][j] = LU_[i][j]; - } else if (i == j) { - L_[i][j] = 1.0; - } else { - L_[i][j] = 0.0; - } - } - } - return L_; - } - - /** Return upper triangular factor - @return U portion of LU factorization. - */ - - Array2D getU () { - Array2D U_(n,n); - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - if (i <= j) { - U_[i][j] = LU_[i][j]; - } else { - U_[i][j] = 0.0; - } - } - } - return U_; - } - - /** Return pivot permutation vector - @return piv - */ - - Array1D getPivot () { - return piv; - } - - - /** Compute determinant using LU factors. - @return determinant of A, or 0 if A is not square. - */ - - Real det () { - if (m != n) { - return Real(0); - } - Real d = Real(pivsign); - for (int j = 0; j < n; j++) { - d *= LU_[j][j]; - } - return d; - } - - /** Solve A*X = B - @param B A Matrix with as many rows as A and any number of columns. - @return X so that L*U*X = B(piv,:), if B is nonconformant, returns - 0x0 (null) array. - */ - - Array2D solve (const Array2D &B) - { - - /* Dimensions: A is mxn, X is nxk, B is mxk */ - - if (B.dim1() != m) { - return Array2D(0,0); - } - if (!isNonsingular()) { - return Array2D(0,0); - } - - // Copy right hand side with pivoting - int nx = B.dim2(); - - - Array2D X = permute_copy(B, piv, 0, nx-1); - - // Solve L*Y = B(piv,:) - for (int k = 0; k < n; k++) { - for (int i = k+1; i < n; i++) { - for (int j = 0; j < nx; j++) { - X[i][j] -= X[k][j]*LU_[i][k]; - } - } - } - // Solve U*X = Y; - for (int k = n-1; k >= 0; k--) { - for (int j = 0; j < nx; j++) { - X[k][j] /= LU_[k][k]; - } - for (int i = 0; i < k; i++) { - for (int j = 0; j < nx; j++) { - X[i][j] -= X[k][j]*LU_[i][k]; - } - } - } - return X; - } - - - /** 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 - of A. - @return x a vector (Array1D> so that L*U*x = b(piv), if B is nonconformant, - returns 0x0 (null) array. - */ - - Array1D solve (const Array1D &b) - { - - /* Dimensions: A is mxn, X is nxk, B is mxk */ - - if (b.dim1() != m) { - return Array1D(); - } - if (!isNonsingular()) { - return Array1D(); - } - - - Array1D x = permute_copy(b, piv); - - // Solve L*Y = B(piv) - for (int k = 0; k < n; k++) { - for (int i = k+1; i < n; i++) { - x[i] -= x[k]*LU_[i][k]; - } - } - - // Solve U*X = Y; - for (int k = n-1; k >= 0; k--) { - x[k] /= LU_[k][k]; - for (int i = 0; i < k; i++) - x[i] -= x[k]*LU_[i][k]; - } - - - return x; - } - -}; /* class LU */ - -} /* namespace JAMA */ - -#endif -/* JAMA_LU_H */ diff --git a/lib/include/tnt/jama_qr.h b/lib/include/tnt/jama_qr.h deleted file mode 100644 index 98d37f5..0000000 --- a/lib/include/tnt/jama_qr.h +++ /dev/null @@ -1,326 +0,0 @@ -#ifndef JAMA_QR_H -#define JAMA_QR_H - -#include "tnt_array1d.h" -#include "tnt_array2d.h" -#include "tnt_math_utils.h" - -namespace JAMA -{ - -/** -

- Classical QR Decompisition: - for an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n - orthogonal matrix Q and an n-by-n upper triangular matrix R so that - A = Q*R. -

- The QR decompostion always exists, even if the matrix does not have - full rank, so the constructor will never fail. The primary use of the - QR decomposition is in the least squares solution of nonsquare systems - of simultaneous linear equations. This will fail if isFullRank() - returns 0 (false). - -

- The Q and R factors can be retrived via the getQ() and getR() - methods. Furthermore, a solve() method is provided to find the - least squares solution of Ax=b using the QR factors. - -

- (Adapted from JAMA, a Java Matrix Library, developed by jointly - by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). -*/ - -template -class QR { - - - /** Array for internal storage of decomposition. - @serial internal array storage. - */ - - TNT::Array2D QR_; - - /** Row and column dimensions. - @serial column dimension. - @serial row dimension. - */ - int m, n; - - /** Array for internal storage of diagonal of R. - @serial diagonal of R. - */ - TNT::Array1D Rdiag; - - -public: - -/** - Create a QR factorization object for A. - - @param A rectangular (m>=n) matrix. -*/ - QR(const TNT::Array2D &A) /* constructor */ - { - QR_ = A.copy(); - m = A.dim1(); - n = A.dim2(); - Rdiag = TNT::Array1D(n); - int i=0, j=0, k=0; - - // Main loop. - for (k = 0; k < n; k++) { - // Compute 2-norm of k-th column without under/overflow. - Real nrm = 0; - for (i = k; i < m; i++) { - nrm = TNT::hypot(nrm,QR_[i][k]); - } - - if (nrm != 0.0) { - // Form k-th Householder vector. - if (QR_[k][k] < 0) { - nrm = -nrm; - } - for (i = k; i < m; i++) { - QR_[i][k] /= nrm; - } - QR_[k][k] += 1.0; - - // Apply transformation to remaining columns. - for (j = k+1; j < n; j++) { - Real s = 0.0; - for (i = k; i < m; i++) { - s += QR_[i][k]*QR_[i][j]; - } - s = -s/QR_[k][k]; - for (i = k; i < m; i++) { - QR_[i][j] += s*QR_[i][k]; - } - } - } - Rdiag[k] = -nrm; - } - } - - -/** - Flag to denote the matrix is of full rank. - - @return 1 if matrix is full rank, 0 otherwise. -*/ - int isFullRank() const - { - for (int j = 0; j < n; j++) - { - if (Rdiag[j] == 0) - return 0; - } - return 1; - } - - - - - /** - - Retreive the Householder vectors from QR factorization - @returns lower trapezoidal matrix whose columns define the reflections - */ - - TNT::Array2D getHouseholder (void) const - { - TNT::Array2D H(m,n); - - /* note: H is completely filled in by algorithm, so - initializaiton of H is not necessary. - */ - for (int i = 0; i < m; i++) - { - for (int j = 0; j < n; j++) - { - if (i >= j) { - H[i][j] = QR_[i][j]; - } else { - H[i][j] = 0.0; - } - } - } - return H; - } - - - - /** Return the upper triangular factor, R, of the QR factorization - @return R - */ - - TNT::Array2D getR() const - { - TNT::Array2D R(n,n); - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - if (i < j) { - R[i][j] = QR_[i][j]; - } else if (i == j) { - R[i][j] = Rdiag[i]; - } else { - R[i][j] = 0.0; - } - } - } - return R; - } - - - - - - /** - Generate and return the (economy-sized) orthogonal factor - @param Q the (ecnomy-sized) orthogonal factor (Q*R=A). - */ - - TNT::Array2D getQ() const - { - int i=0, j=0, k=0; - - TNT::Array2D Q(m,n); - for (k = n-1; k >= 0; k--) { - for (i = 0; i < m; i++) { - Q[i][k] = 0.0; - } - Q[k][k] = 1.0; - for (j = k; j < n; j++) { - if (QR_[k][k] != 0) { - Real s = 0.0; - for (i = k; i < m; i++) { - s += QR_[i][k]*Q[i][j]; - } - s = -s/QR_[k][k]; - for (i = k; i < m; i++) { - Q[i][j] += s*QR_[i][k]; - } - } - } - } - return Q; - } - - - /** Least squares solution of A*x = b - @param B m-length array (vector). - @return x n-length array (vector) that minimizes the two norm of Q*R*X-B. - If B is non-conformant, or if QR.isFullRank() is false, - the routine returns a null (0-length) vector. - */ - - TNT::Array1D solve(const TNT::Array1D &b) const - { - if (b.dim1() != m) /* arrays must be conformant */ - return TNT::Array1D(); - - if ( !isFullRank() ) /* matrix is rank deficient */ - { - return TNT::Array1D(); - } - - TNT::Array1D x = b.copy(); - - // Compute Y = transpose(Q)*b - for (int k = 0; k < n; k++) - { - Real s = 0.0; - for (int i = k; i < m; i++) - { - s += QR_[i][k]*x[i]; - } - s = -s/QR_[k][k]; - for (int i = k; i < m; i++) - { - x[i] += s*QR_[i][k]; - } - } - // Solve R*X = Y; - for (int k = n-1; k >= 0; k--) - { - x[k] /= Rdiag[k]; - for (int i = 0; i < k; i++) { - x[i] -= x[k]*QR_[i][k]; - } - } - - - /* return n x nx portion of X */ - TNT::Array1D x_(n); - for (int i=0; i solve(const TNT::Array2D &B) const - { - if (B.dim1() != m) /* arrays must be conformant */ - return TNT::Array2D(0,0); - - if ( !isFullRank() ) /* matrix is rank deficient */ - { - return TNT::Array2D(0,0); - } - - int nx = B.dim2(); - TNT::Array2D X = B.copy(); - int i=0, j=0, k=0; - - // Compute Y = transpose(Q)*B - for (k = 0; k < n; k++) { - for (j = 0; j < nx; j++) { - Real s = 0.0; - for (i = k; i < m; i++) { - s += QR_[i][k]*X[i][j]; - } - s = -s/QR_[k][k]; - for (i = k; i < m; i++) { - X[i][j] += s*QR_[i][k]; - } - } - } - // Solve R*X = Y; - for (k = n-1; k >= 0; k--) { - for (j = 0; j < nx; j++) { - X[k][j] /= Rdiag[k]; - } - for (i = 0; i < k; i++) { - for (j = 0; j < nx; j++) { - X[i][j] -= X[k][j]*QR_[i][k]; - } - } - } - - - /* return n x nx portion of X */ - TNT::Array2D X_(n,nx); - for (i=0; i -// for min(), max() below -#include -// for abs() below - -using namespace TNT; -using namespace std; - - -// Modification by Willem Jan Palenstijn, 2010-03-11: -// Use std::min() instead of min(), std::max() instead of max() - - -namespace JAMA -{ - /** Singular Value Decomposition. -

- For an m-by-n matrix A with m >= n, the singular value decomposition is - an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and - an n-by-n orthogonal matrix V so that A = U*S*V'. -

- The singular values, sigma[k] = S[k][k], are ordered so that - sigma[0] >= sigma[1] >= ... >= sigma[n-1]. -

- The singular value decompostion always exists, so the constructor will - never fail. The matrix condition number and the effective numerical - rank can be computed from this decomposition. - -

- (Adapted from JAMA, a Java Matrix Library, developed by jointly - by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). - */ -template -class SVD -{ - - - Array2D U, V; - Array1D s; - int m, n; - - public: - - - SVD (const Array2D &Arg) { - - - m = Arg.dim1(); - n = Arg.dim2(); - int nu = std::min(m,n); - s = Array1D(std::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 wantu = 1; /* boolean */ - int wantv = 1; /* boolean */ - int i=0, j=0, k=0; - - // Reduce A to bidiagonal form, storing the diagonal elements - // in s and the super-diagonal elements in e. - - 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 - // place the k-th diagonal in s[k]. - // Compute 2-norm of k-th column without under/overflow. - s[k] = 0; - for (i = k; i < m; i++) { - s[k] = hypot(s[k],A[i][k]); - } - if (s[k] != 0.0) { - if (A[k][k] < 0.0) { - s[k] = -s[k]; - } - for (i = k; i < m; i++) { - A[i][k] /= s[k]; - } - A[k][k] += 1.0; - } - s[k] = -s[k]; - } - for (j = k+1; j < n; j++) { - if ((k < nct) && (s[k] != 0.0)) { - - // Apply the transformation. - - Real t(0.0); - for (i = k; i < m; i++) { - t += A[i][k]*A[i][j]; - } - t = -t/A[k][k]; - for (i = k; i < m; i++) { - A[i][j] += t*A[i][k]; - } - } - - // Place the k-th row of A into e for the - // subsequent calculation of the row transformation. - - e[j] = A[k][j]; - } - if (wantu & (k < nct)) { - - // Place the transformation in U for subsequent back - // multiplication. - - for (i = k; i < m; i++) { - U[i][k] = A[i][k]; - } - } - if (k < nrt) { - - // Compute the k-th row transformation and place the - // k-th super-diagonal in e[k]. - // Compute 2-norm without under/overflow. - e[k] = 0; - for (i = k+1; i < n; i++) { - e[k] = hypot(e[k],e[i]); - } - if (e[k] != 0.0) { - if (e[k+1] < 0.0) { - e[k] = -e[k]; - } - for (i = k+1; i < n; i++) { - e[i] /= e[k]; - } - e[k+1] += 1.0; - } - e[k] = -e[k]; - if ((k+1 < m) & (e[k] != 0.0)) { - - // Apply the transformation. - - for (i = k+1; i < m; i++) { - work[i] = 0.0; - } - for (j = k+1; j < n; j++) { - for (i = k+1; i < m; i++) { - work[i] += e[j]*A[i][j]; - } - } - for (j = k+1; j < n; j++) { - Real t(-e[j]/e[k+1]); - for (i = k+1; i < m; i++) { - A[i][j] += t*work[i]; - } - } - } - if (wantv) { - - // Place the transformation in V for subsequent - // back multiplication. - - for (i = k+1; i < n; i++) { - V[i][k] = e[i]; - } - } - } - } - - // Set up the final bidiagonal matrix or order p. - - int p = std::min(n,m+1); - if (nct < n) { - s[nct] = A[nct][nct]; - } - if (m < p) { - s[p-1] = 0.0; - } - if (nrt+1 < p) { - e[nrt] = A[nrt][p-1]; - } - e[p-1] = 0.0; - - // If required, generate U. - - if (wantu) { - for (j = nct; j < nu; j++) { - for (i = 0; i < m; i++) { - U[i][j] = 0.0; - } - U[j][j] = 1.0; - } - for (k = nct-1; k >= 0; k--) { - if (s[k] != 0.0) { - for (j = k+1; j < nu; j++) { - Real t(0.0); - for (i = k; i < m; i++) { - t += U[i][k]*U[i][j]; - } - t = -t/U[k][k]; - for (i = k; i < m; i++) { - U[i][j] += t*U[i][k]; - } - } - for (i = k; i < m; i++ ) { - U[i][k] = -U[i][k]; - } - U[k][k] = 1.0 + U[k][k]; - for (i = 0; i < k-1; i++) { - U[i][k] = 0.0; - } - } else { - for (i = 0; i < m; i++) { - U[i][k] = 0.0; - } - U[k][k] = 1.0; - } - } - } - - // If required, generate V. - - if (wantv) { - for (k = n-1; k >= 0; k--) { - if ((k < nrt) & (e[k] != 0.0)) { - for (j = k+1; j < nu; j++) { - Real t(0.0); - for (i = k+1; i < n; i++) { - t += V[i][k]*V[i][j]; - } - t = -t/V[k+1][k]; - for (i = k+1; i < n; i++) { - V[i][j] += t*V[i][k]; - } - } - } - for (i = 0; i < n; i++) { - V[i][k] = 0.0; - } - V[k][k] = 1.0; - } - } - - // Main iteration loop for the singular values. - - int pp = p-1; - int iter = 0; - Real eps(pow(2.0,-52.0)); - while (p > 0) { - int k=0; - int kase=0; - - // Here is where a test for too many iterations would go. - - // This section of the program inspects for - // negligible elements in the s and e arrays. On - // completion the variables kase and k are set as follows. - - // kase = 1 if s(p) and e[k-1] are negligible and k

= -1; k--) { - if (k == -1) { - break; - } - if (abs(e[k]) <= eps*(abs(s[k]) + abs(s[k+1]))) { - e[k] = 0.0; - break; - } - } - if (k == p-2) { - kase = 4; - } else { - int ks; - for (ks = p-1; ks >= k; ks--) { - 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) { - s[ks] = 0.0; - break; - } - } - if (ks == k) { - kase = 3; - } else if (ks == p-1) { - kase = 1; - } else { - kase = 2; - k = ks; - } - } - k++; - - // Perform the task indicated by kase. - - switch (kase) { - - // Deflate negligible s(p). - - case 1: { - Real f(e[p-2]); - e[p-2] = 0.0; - for (j = p-2; j >= k; j--) { - 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]; - e[j-1] = cs*e[j-1]; - } - if (wantv) { - for (i = 0; i < n; i++) { - t = cs*V[i][j] + sn*V[i][p-1]; - V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1]; - V[i][j] = t; - } - } - } - } - break; - - // Split at negligible s(k). - - case 2: { - Real f(e[k-1]); - e[k-1] = 0.0; - for (j = k; j < p; j++) { - 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]; - if (wantu) { - for (i = 0; i < m; i++) { - t = cs*U[i][j] + sn*U[i][k-1]; - U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1]; - U[i][j] = t; - } - } - } - } - break; - - // Perform one qr step. - - case 3: { - - // Calculate the shift. - - Real 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])); - 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) { - shift = -shift; - } - shift = c/(b + shift); - } - Real f = (sk + sp)*(sk - sp) + shift; - Real g = sk*ek; - - // Chase zeros. - - for (j = k; j < p-1; j++) { - Real t = hypot(f,g); - Real cs = f/t; - Real sn = g/t; - if (j != k) { - e[j-1] = t; - } - f = cs*s[j] + sn*e[j]; - e[j] = cs*e[j] - sn*s[j]; - g = sn*s[j+1]; - s[j+1] = cs*s[j+1]; - if (wantv) { - for (i = 0; i < n; i++) { - t = cs*V[i][j] + sn*V[i][j+1]; - V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1]; - V[i][j] = t; - } - } - t = hypot(f,g); - cs = f/t; - sn = g/t; - s[j] = t; - f = cs*e[j] + sn*s[j+1]; - s[j+1] = -sn*e[j] + cs*s[j+1]; - g = sn*e[j+1]; - e[j+1] = cs*e[j+1]; - if (wantu && (j < m-1)) { - for (i = 0; i < m; i++) { - t = cs*U[i][j] + sn*U[i][j+1]; - U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1]; - U[i][j] = t; - } - } - } - e[p-2] = f; - iter = iter + 1; - } - break; - - // Convergence. - - case 4: { - - // Make the singular values positive. - - if (s[k] <= 0.0) { - s[k] = (s[k] < 0.0 ? -s[k] : 0.0); - if (wantv) { - for (i = 0; i <= pp; i++) { - V[i][k] = -V[i][k]; - } - } - } - - // Order the singular values. - - while (k < pp) { - if (s[k] >= s[k+1]) { - break; - } - Real t = s[k]; - s[k] = s[k+1]; - s[k+1] = t; - if (wantv && (k < n-1)) { - for (i = 0; i < n; i++) { - t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t; - } - } - if (wantu && (k < m-1)) { - for (i = 0; i < m; i++) { - t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t; - } - } - k++; - } - iter = 0; - p--; - } - break; - } - } - } - - - void getU (Array2D &A) - { - int minm = std::min(m+1,n); - - A = Array2D(m, minm); - - for (int i=0; i &A) - { - A = V; - } - - /** Return the one-dimensional array of singular values */ - - void getSingularValues (Array1D &x) - { - x = s; - } - - /** Return the diagonal matrix of singular values - @return S - */ - - void getS (Array2D &A) { - A = Array2D(n,n); - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - A[i][j] = 0.0; - } - A[i][i] = s[i]; - } - } - - /** Two norm (max(S)) */ - - Real norm2 () { - return s[0]; - } - - /** Two norm of condition number (max(S)/min(S)) */ - - Real cond () { - return s[0]/s[std::min(m,n)-1]; - } - - /** Effective numerical matrix rank - @return Number of nonnegligible singular values. - */ - - int rank () - { - Real eps = 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) { - r++; - } - } - return r; - } -}; - -} -#endif -// JAMA_SVD_H diff --git a/lib/include/tnt/tnt.h b/lib/include/tnt/tnt.h deleted file mode 100644 index 92463e0..0000000 --- a/lib/include/tnt/tnt.h +++ /dev/null @@ -1,64 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - -#ifndef TNT_H -#define TNT_H - - - -//--------------------------------------------------------------------- -// Define this macro if you want TNT to track some of the out-of-bounds -// indexing. This can encur a small run-time overhead, but is recommended -// while developing code. It can be turned off for production runs. -// -// #define TNT_BOUNDS_CHECK -//--------------------------------------------------------------------- -// - -//#define TNT_BOUNDS_CHECK - - - -#include "tnt_version.h" -#include "tnt_math_utils.h" -#include "tnt_array1d.h" -#include "tnt_array2d.h" -#include "tnt_array3d.h" -#include "tnt_array1d_utils.h" -#include "tnt_array2d_utils.h" -#include "tnt_array3d_utils.h" - -#include "tnt_fortran_array1d.h" -#include "tnt_fortran_array2d.h" -#include "tnt_fortran_array3d.h" -#include "tnt_fortran_array1d_utils.h" -#include "tnt_fortran_array2d_utils.h" -#include "tnt_fortran_array3d_utils.h" - -#include "tnt_sparse_matrix_csr.h" - -#include "tnt_stopwatch.h" -#include "tnt_subscript.h" -#include "tnt_vec.h" -#include "tnt_cmat.h" - - -#endif -// TNT_H diff --git a/lib/include/tnt/tnt_array1d.h b/lib/include/tnt/tnt_array1d.h deleted file mode 100644 index 858df57..0000000 --- a/lib/include/tnt/tnt_array1d.h +++ /dev/null @@ -1,278 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT) -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - - -#ifndef TNT_ARRAY1D_H -#define TNT_ARRAY1D_H - -//#include -#include - -#ifdef TNT_BOUNDS_CHECK -#include -#endif - - -#include "tnt_i_refvec.h" - -namespace TNT -{ - -template -class Array1D -{ - - private: - - /* ... */ - i_refvec v_; - int n_; - T* data_; /* this normally points to v_.begin(), but - * could also point to a portion (subvector) - * of v_. - */ - - void copy_(T* p, const T* q, int len) const; - void set_(T* begin, T* end, const T& val); - - - public: - - typedef T value_type; - - - Array1D(); - explicit Array1D(int n); - Array1D(int n, const T &a); - Array1D(int n, T *a); - inline Array1D(const Array1D &A); - inline operator T*(); - inline operator const T*(); - inline Array1D & operator=(const T &a); - inline Array1D & operator=(const Array1D &A); - inline Array1D & ref(const Array1D &A); - Array1D copy() const; - Array1D & inject(const Array1D & A); - inline T& operator[](int i); - inline const T& operator[](int i) const; - inline int dim1() const; - inline int dim() const; - ~Array1D(); - - - /* ... extended interface ... */ - - inline int ref_count() const; - inline Array1D subarray(int i0, int i1); - -}; - - - - -template -Array1D::Array1D() : v_(), n_(0), data_(0) {} - -template -Array1D::Array1D(const Array1D &A) : v_(A.v_), n_(A.n_), - data_(A.data_) -{ -#ifdef TNT_DEBUG - std::cout << "Created Array1D(const Array1D &A) \n"; -#endif - -} - - -template -Array1D::Array1D(int n) : v_(n), n_(n), data_(v_.begin()) -{ -#ifdef TNT_DEBUG - std::cout << "Created Array1D(int n) \n"; -#endif -} - -template -Array1D::Array1D(int n, const T &val) : v_(n), n_(n), data_(v_.begin()) -{ -#ifdef TNT_DEBUG - std::cout << "Created Array1D(int n, const T& val) \n"; -#endif - set_(data_, data_+ n, val); - -} - -template -Array1D::Array1D(int n, T *a) : v_(a), n_(n) , data_(v_.begin()) -{ -#ifdef TNT_DEBUG - std::cout << "Created Array1D(int n, T* a) \n"; -#endif -} - -template -inline Array1D::operator T*() -{ - return &(v_[0]); -} - - -template -inline Array1D::operator const T*() -{ - return &(v_[0]); -} - - - -template -inline T& Array1D::operator[](int i) -{ -#ifdef TNT_BOUNDS_CHECK - assert(i>= 0); - assert(i < n_); -#endif - return data_[i]; -} - -template -inline const T& Array1D::operator[](int i) const -{ -#ifdef TNT_BOUNDS_CHECK - assert(i>= 0); - assert(i < n_); -#endif - return data_[i]; -} - - - - -template -Array1D & Array1D::operator=(const T &a) -{ - set_(data_, data_+n_, a); - return *this; -} - -template -Array1D Array1D::copy() const -{ - Array1D A( n_); - copy_(A.data_, data_, n_); - - return A; -} - - -template -Array1D & Array1D::inject(const Array1D &A) -{ - if (A.n_ == n_) - copy_(data_, A.data_, n_); - - return *this; -} - - - - - -template -Array1D & Array1D::ref(const Array1D &A) -{ - if (this != &A) - { - v_ = A.v_; /* operator= handles the reference counting. */ - n_ = A.n_; - data_ = A.data_; - - } - return *this; -} - -template -Array1D & Array1D::operator=(const Array1D &A) -{ - return ref(A); -} - -template -inline int Array1D::dim1() const { return n_; } - -template -inline int Array1D::dim() const { return n_; } - -template -Array1D::~Array1D() {} - - -/* ............................ exented interface ......................*/ - -template -inline int Array1D::ref_count() const -{ - return v_.ref_count(); -} - -template -inline Array1D Array1D::subarray(int i0, int i1) -{ - if ((i0 > 0) && (i1 < n_) || (i0 <= i1)) - { - Array1D X(*this); /* create a new instance of this array. */ - X.n_ = i1-i0+1; - X.data_ += i0; - - return X; - } - else - { - return Array1D(); - } -} - - -/* private internal functions */ - - -template -void Array1D::set_(T* begin, T* end, const T& a) -{ - for (T* p=begin; p -void Array1D::copy_(T* p, const T* q, int len) const -{ - T *end = p + len; - while (p -#include - -namespace TNT -{ - - -template -std::ostream& operator<<(std::ostream &s, const Array1D &A) -{ - int N=A.dim1(); - -#ifdef TNT_DEBUG - s << "addr: " << (void *) &A[0] << "\n"; -#endif - s << N << "\n"; - for (int j=0; j -std::istream& operator>>(std::istream &s, Array1D &A) -{ - int N; - s >> N; - - Array1D B(N); - for (int i=0; i> B[i]; - A = B; - return s; -} - - - -template -Array1D operator+(const Array1D &A, const Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() != n ) - return Array1D(); - - else - { - Array1D C(n); - - for (int i=0; i -Array1D operator-(const Array1D &A, const Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() != n ) - return Array1D(); - - else - { - Array1D C(n); - - for (int i=0; i -Array1D operator*(const Array1D &A, const Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() != n ) - return Array1D(); - - else - { - Array1D C(n); - - for (int i=0; i -Array1D operator/(const Array1D &A, const Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() != n ) - return Array1D(); - - else - { - Array1D C(n); - - for (int i=0; i -Array1D& operator+=(Array1D &A, const Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() == n) - { - for (int i=0; i -Array1D& operator-=(Array1D &A, const Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() == n) - { - for (int i=0; i -Array1D& operator*=(Array1D &A, const Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() == n) - { - for (int i=0; i -Array1D& operator/=(Array1D &A, const Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() == n) - { - for (int i=0; i -#include -#ifdef TNT_BOUNDS_CHECK -#include -#endif - -#include "tnt_array1d.h" - -namespace TNT -{ - -template -class Array2D -{ - - - private: - - - - Array1D data_; - Array1D v_; - int m_; - int n_; - - public: - - typedef T value_type; - Array2D(); - Array2D(int m, int n); - Array2D(int m, int n, T *a); - Array2D(int m, int n, const T &a); - inline Array2D(const Array2D &A); - inline operator T**(); - inline operator const T**(); - inline Array2D & operator=(const T &a); - inline Array2D & operator=(const Array2D &A); - inline Array2D & ref(const Array2D &A); - Array2D copy() const; - Array2D & inject(const Array2D & A); - inline T* operator[](int i); - inline const T* operator[](int i) const; - inline int dim1() const; - inline int dim2() const; - ~Array2D(); - - /* extended interface (not part of the standard) */ - - - inline int ref_count(); - inline int ref_count_data(); - inline int ref_count_dim1(); - Array2D subarray(int i0, int i1, int j0, int j1); - -}; - - -template -Array2D::Array2D() : data_(), v_(), m_(0), n_(0) {} - -template -Array2D::Array2D(const Array2D &A) : data_(A.data_), v_(A.v_), - m_(A.m_), n_(A.n_) {} - - - - -template -Array2D::Array2D(int m, int n) : data_(m*n), v_(m), m_(m), n_(n) -{ - if (m>0 && n>0) - { - T* p = &(data_[0]); - for (int i=0; i -Array2D::Array2D(int m, int n, const T &val) : data_(m*n), v_(m), - m_(m), n_(n) -{ - if (m>0 && n>0) - { - data_ = val; - T* p = &(data_[0]); - for (int i=0; i -Array2D::Array2D(int m, int n, T *a) : data_(m*n, a), v_(m), m_(m), n_(n) -{ - if (m>0 && n>0) - { - T* p = &(data_[0]); - - for (int i=0; i -inline T* Array2D::operator[](int i) -{ -#ifdef TNT_BOUNDS_CHECK - assert(i >= 0); - assert(i < m_); -#endif - -return v_[i]; - -} - - -template -inline const T* Array2D::operator[](int i) const -{ -#ifdef TNT_BOUNDS_CHECK - assert(i >= 0); - assert(i < m_); -#endif - -return v_[i]; - -} - -template -Array2D & Array2D::operator=(const T &a) -{ - /* non-optimzied, but will work with subarrays in future verions */ - - for (int i=0; i -Array2D Array2D::copy() const -{ - Array2D A(m_, n_); - - for (int i=0; i -Array2D & Array2D::inject(const Array2D &A) -{ - if (A.m_ == m_ && A.n_ == n_) - { - for (int i=0; i -Array2D & Array2D::ref(const Array2D &A) -{ - if (this != &A) - { - v_ = A.v_; - data_ = A.data_; - m_ = A.m_; - n_ = A.n_; - - } - return *this; -} - - - -template -Array2D & Array2D::operator=(const Array2D &A) -{ - return ref(A); -} - -template -inline int Array2D::dim1() const { return m_; } - -template -inline int Array2D::dim2() const { return n_; } - - -template -Array2D::~Array2D() {} - - - - -template -inline Array2D::operator T**() -{ - return &(v_[0]); -} -template -inline Array2D::operator const T**() -{ - return &(v_[0]); -} - -/* ............... extended interface ............... */ -/** - Create a new view to a subarray defined by the boundaries - [i0][i0] and [i1][j1]. The size of the subarray is - (i1-i0) by (j1-j0). If either of these lengths are zero - or negative, the subarray view is null. - -*/ -template -Array2D Array2D::subarray(int i0, int i1, int j0, int j1) -{ - Array2D A; - int m = i1-i0+1; - int n = j1-j0+1; - - /* if either length is zero or negative, this is an invalide - subarray. return a null view. - */ - if (m<1 || n<1) - return A; - - A.data_ = data_; - A.m_ = m; - A.n_ = n; - A.v_ = Array1D(m); - T* p = &(data_[0]) + i0 * n_ + j0; - for (int i=0; i -inline int Array2D::ref_count() -{ - return ref_count_data(); -} - - - -template -inline int Array2D::ref_count_data() -{ - return data_.ref_count(); -} - -template -inline int Array2D::ref_count_dim1() -{ - return v_.ref_count(); -} - - - - -} /* namespace TNT */ - -#endif -/* TNT_ARRAY2D_H */ - diff --git a/lib/include/tnt/tnt_array2d_utils.h b/lib/include/tnt/tnt_array2d_utils.h deleted file mode 100644 index 7041ed3..0000000 --- a/lib/include/tnt/tnt_array2d_utils.h +++ /dev/null @@ -1,287 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT) -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - -#ifndef TNT_ARRAY2D_UTILS_H -#define TNT_ARRAY2D_UTILS_H - -#include -#include - -namespace TNT -{ - - -template -std::ostream& operator<<(std::ostream &s, const Array2D &A) -{ - int M=A.dim1(); - int N=A.dim2(); - - s << M << " " << N << "\n"; - - for (int i=0; i -std::istream& operator>>(std::istream &s, Array2D &A) -{ - - int M, N; - - s >> M >> N; - - Array2D B(M,N); - - for (int i=0; i> B[i][j]; - } - - A = B; - return s; -} - - -template -Array2D operator+(const Array2D &A, const Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() != m || B.dim2() != n ) - return Array2D(); - - else - { - Array2D C(m,n); - - for (int i=0; i -Array2D operator-(const Array2D &A, const Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() != m || B.dim2() != n ) - return Array2D(); - - else - { - Array2D C(m,n); - - for (int i=0; i -Array2D operator*(const Array2D &A, const Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() != m || B.dim2() != n ) - return Array2D(); - - else - { - Array2D C(m,n); - - for (int i=0; i -Array2D operator/(const Array2D &A, const Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() != m || B.dim2() != n ) - return Array2D(); - - else - { - Array2D C(m,n); - - for (int i=0; i -Array2D& operator+=(Array2D &A, const Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() == m || B.dim2() == n ) - { - for (int i=0; i -Array2D& operator-=(Array2D &A, const Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() == m || B.dim2() == n ) - { - for (int i=0; i -Array2D& operator*=(Array2D &A, const Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() == m || B.dim2() == n ) - { - for (int i=0; i -Array2D& operator/=(Array2D &A, const Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() == m || B.dim2() == n ) - { - for (int i=0; i -Array2D matmult(const Array2D &A, const Array2D &B) -{ - if (A.dim2() != B.dim1()) - return Array2D(); - - int M = A.dim1(); - int N = A.dim2(); - int K = B.dim2(); - - Array2D C(M,K); - - for (int i=0; i -#include -#ifdef TNT_BOUNDS_CHECK -#include -#endif - -#include "tnt_array1d.h" -#include "tnt_array2d.h" - -namespace TNT -{ - -template -class Array3D -{ - - - private: - Array1D data_; - Array2D v_; - int m_; - int n_; - int g_; - - - public: - - typedef T value_type; - - Array3D(); - Array3D(int m, int n, int g); - Array3D(int m, int n, int g, T val); - Array3D(int m, int n, int g, T *a); - - inline operator T***(); - inline operator const T***(); - inline Array3D(const Array3D &A); - inline Array3D & operator=(const T &a); - inline Array3D & operator=(const Array3D &A); - inline Array3D & ref(const Array3D &A); - Array3D copy() const; - Array3D & inject(const Array3D & A); - - inline T** operator[](int i); - inline const T* const * operator[](int i) const; - inline int dim1() const; - inline int dim2() const; - inline int dim3() const; - ~Array3D(); - - /* extended interface */ - - inline int ref_count(){ return data_.ref_count(); } - Array3D subarray(int i0, int i1, int j0, int j1, - int k0, int k1); -}; - -template -Array3D::Array3D() : data_(), v_(), m_(0), n_(0) {} - -template -Array3D::Array3D(const Array3D &A) : data_(A.data_), - v_(A.v_), m_(A.m_), n_(A.n_), g_(A.g_) -{ -} - - - -template -Array3D::Array3D(int m, int n, int g) : data_(m*n*g), v_(m,n), - m_(m), n_(n), g_(g) -{ - - if (m>0 && n>0 && g>0) - { - T* p = & (data_[0]); - int ng = n_*g_; - - for (int i=0; i -Array3D::Array3D(int m, int n, int g, T val) : data_(m*n*g, val), - v_(m,n), m_(m), n_(n), g_(g) -{ - if (m>0 && n>0 && g>0) - { - - T* p = & (data_[0]); - int ng = n_*g_; - - for (int i=0; i -Array3D::Array3D(int m, int n, int g, T* a) : - data_(m*n*g, a), v_(m,n), m_(m), n_(n), g_(g) -{ - - if (m>0 && n>0 && g>0) - { - T* p = & (data_[0]); - int ng = n_*g_; - - for (int i=0; i -inline T** Array3D::operator[](int i) -{ -#ifdef TNT_BOUNDS_CHECK - assert(i >= 0); - assert(i < m_); -#endif - -return v_[i]; - -} - -template -inline const T* const * Array3D::operator[](int i) const -{ return v_[i]; } - -template -Array3D & Array3D::operator=(const T &a) -{ - for (int i=0; i -Array3D Array3D::copy() const -{ - Array3D A(m_, n_, g_); - for (int i=0; i -Array3D & Array3D::inject(const Array3D &A) -{ - if (A.m_ == m_ && A.n_ == n_ && A.g_ == g_) - - for (int i=0; i -Array3D & Array3D::ref(const Array3D &A) -{ - if (this != &A) - { - m_ = A.m_; - n_ = A.n_; - g_ = A.g_; - v_ = A.v_; - data_ = A.data_; - } - return *this; -} - -template -Array3D & Array3D::operator=(const Array3D &A) -{ - return ref(A); -} - - -template -inline int Array3D::dim1() const { return m_; } - -template -inline int Array3D::dim2() const { return n_; } - -template -inline int Array3D::dim3() const { return g_; } - - - -template -Array3D::~Array3D() {} - -template -inline Array3D::operator T***() -{ - return v_; -} - - -template -inline Array3D::operator const T***() -{ - return v_; -} - -/* extended interface */ -template -Array3D Array3D::subarray(int i0, int i1, int j0, - int j1, int k0, int k1) -{ - - /* check that ranges are valid. */ - if (!( 0 <= i0 && i0 <= i1 && i1 < m_ && - 0 <= j0 && j0 <= j1 && j1 < n_ && - 0 <= k0 && k0 <= k1 && k1 < g_)) - return Array3D(); /* null array */ - - - Array3D A; - A.data_ = data_; - A.m_ = i1-i0+1; - A.n_ = j1-j0+1; - A.g_ = k1-k0+1; - A.v_ = Array2D(A.m_,A.n_); - T* p = &(data_[0]) + i0*n_*g_ + j0*g_ + k0; - - for (int i=0; i -#include - -namespace TNT -{ - - -template -std::ostream& operator<<(std::ostream &s, const Array3D &A) -{ - int M=A.dim1(); - int N=A.dim2(); - int K=A.dim3(); - - s << M << " " << N << " " << K << "\n"; - - for (int i=0; i -std::istream& operator>>(std::istream &s, Array3D &A) -{ - - int M, N, K; - - s >> M >> N >> K; - - Array3D B(M,N,K); - - for (int i=0; i> B[i][j][k]; - - A = B; - return s; -} - - - -template -Array3D operator+(const Array3D &A, const Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) - return Array3D(); - - else - { - Array3D C(m,n,p); - - for (int i=0; i -Array3D operator-(const Array3D &A, const Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) - return Array3D(); - - else - { - Array3D C(m,n,p); - - for (int i=0; i -Array3D operator*(const Array3D &A, const Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) - return Array3D(); - - else - { - Array3D C(m,n,p); - - for (int i=0; i -Array3D operator/(const Array3D &A, const Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) - return Array3D(); - - else - { - Array3D C(m,n,p); - - for (int i=0; i -Array3D& operator+=(Array3D &A, const Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) - { - for (int i=0; i -Array3D& operator-=(Array3D &A, const Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) - { - for (int i=0; i -Array3D& operator*=(Array3D &A, const Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) - { - for (int i=0; i -Array3D& operator/=(Array3D &A, const Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) - { - for (int i=0; i -#include -#include -#include - -namespace TNT -{ - - -template -class Matrix -{ - - - public: - - typedef Subscript size_type; - typedef T value_type; - typedef T element_type; - typedef T* pointer; - typedef T* iterator; - typedef T& reference; - typedef const T* const_iterator; - typedef const T& const_reference; - - Subscript lbound() const { return 1;} - - protected: - Subscript m_; - Subscript n_; - Subscript mn_; // total size - T* v_; - T** row_; - T* vm1_ ; // these point to the same data, but are 1-based - T** rowm1_; - - // internal helper function to create the array - // of row pointers - - void initialize(Subscript M, Subscript N) - { - mn_ = M*N; - m_ = M; - n_ = N; - - v_ = new T[mn_]; - row_ = new T*[M]; - rowm1_ = new T*[M]; - - assert(v_ != NULL); - assert(row_ != NULL); - assert(rowm1_ != NULL); - - T* p = v_; - vm1_ = v_ - 1; - for (Subscript i=0; i &A) - { - initialize(A.m_, A.n_); - copy(A.v_); - } - - Matrix(Subscript M, Subscript N, const T& value = T()) - { - initialize(M,N); - set(value); - } - - Matrix(Subscript M, Subscript N, const T* v) - { - initialize(M,N); - copy(v); - } - - Matrix(Subscript M, Subscript N, const char *s) - { - initialize(M,N); - //std::istrstream ins(s); - std::istringstream ins(s); - - Subscript i, j; - - for (i=0; i> row_[i][j]; - } - - // destructor - // - ~Matrix() - { - destroy(); - } - - - // reallocating - // - Matrix& newsize(Subscript M, Subscript N) - { - if (num_rows() == M && num_cols() == N) - return *this; - - destroy(); - initialize(M,N); - - return *this; - } - - - - - // assignments - // - Matrix& operator=(const Matrix &A) - { - if (v_ == A.v_) - return *this; - - if (m_ == A.m_ && n_ == A.n_) // no need to re-alloc - copy(A.v_); - - else - { - destroy(); - initialize(A.m_, A.n_); - copy(A.v_); - } - - return *this; - } - - Matrix& operator=(const T& scalar) - { - set(scalar); - return *this; - } - - - Subscript dim(Subscript d) const - { -#ifdef TNT_BOUNDS_CHECK - assert( d >= 1); - assert( d <= 2); -#endif - return (d==1) ? m_ : ((d==2) ? n_ : 0); - } - - Subscript num_rows() const { return m_; } - Subscript num_cols() const { return n_; } - - - - - inline T* operator[](Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i); - assert(i < m_) ; -#endif - return row_[i]; - } - - inline const T* operator[](Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i); - assert(i < m_) ; -#endif - return row_[i]; - } - - inline reference operator()(Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= mn_) ; -#endif - return vm1_[i]; - } - - inline const_reference operator()(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= mn_) ; -#endif - return vm1_[i]; - } - - - - inline reference operator()(Subscript i, Subscript j) - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= m_) ; - assert(1<=j); - assert(j <= n_); -#endif - return rowm1_[i][j]; - } - - - - inline const_reference operator() (Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= m_) ; - assert(1<=j); - assert(j <= n_); -#endif - return rowm1_[i][j]; - } - - - - -}; - - -/* *************************** I/O ********************************/ - -template -std::ostream& operator<<(std::ostream &s, const Matrix &A) -{ - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - s << M << " " << N << "\n"; - - for (Subscript i=0; i -std::istream& operator>>(std::istream &s, Matrix &A) -{ - - Subscript M, N; - - s >> M >> N; - - if ( !(M == A.num_rows() && N == A.num_cols() )) - { - A.newsize(M,N); - } - - - for (Subscript i=0; i> A[i][j]; - } - - - return s; -} - -// *******************[ basic matrix algorithms ]*************************** - - -template -Matrix operator+(const Matrix &A, - const Matrix &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==B.num_rows()); - assert(N==B.num_cols()); - - Matrix tmp(M,N); - Subscript i,j; - - for (i=0; i -Matrix operator-(const Matrix &A, - const Matrix &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==B.num_rows()); - assert(N==B.num_cols()); - - Matrix tmp(M,N); - Subscript i,j; - - for (i=0; i -Matrix mult_element(const Matrix &A, - const Matrix &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==B.num_rows()); - assert(N==B.num_cols()); - - Matrix tmp(M,N); - Subscript i,j; - - for (i=0; i -Matrix transpose(const Matrix &A) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - Matrix S(N,M); - Subscript i, j; - - for (i=0; i -inline Matrix matmult(const Matrix &A, - const Matrix &B) -{ - -#ifdef TNT_BOUNDS_CHECK - assert(A.num_cols() == B.num_rows()); -#endif - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - Subscript K = B.num_cols(); - - Matrix tmp(M,K); - T sum; - - for (Subscript i=0; i -inline Matrix operator*(const Matrix &A, - const Matrix &B) -{ - return matmult(A,B); -} - -template -inline int matmult(Matrix& C, const Matrix &A, - const Matrix &B) -{ - - assert(A.num_cols() == B.num_rows()); - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - Subscript K = B.num_cols(); - - C.newsize(M,K); - - T sum; - - const T* row_i; - const T* col_k; - - for (Subscript i=0; i -Vector matmult(const Matrix &A, const Vector &x) -{ - -#ifdef TNT_BOUNDS_CHECK - assert(A.num_cols() == x.dim()); -#endif - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - Vector tmp(M); - T sum; - - for (Subscript i=0; i -inline Vector operator*(const Matrix &A, const Vector &x) -{ - return matmult(A,x); -} - -} // namespace TNT - -#endif -// CMAT_H diff --git a/lib/include/tnt/tnt_fortran_array1d.h b/lib/include/tnt/tnt_fortran_array1d.h deleted file mode 100644 index ad3bba0..0000000 --- a/lib/include/tnt/tnt_fortran_array1d.h +++ /dev/null @@ -1,267 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT) -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - - -#ifndef TNT_FORTRAN_ARRAY1D_H -#define TNT_FORTRAN_ARRAY1D_H - -#include -#include - -#ifdef TNT_BOUNDS_CHECK -#include -#endif - - -#include "tnt_i_refvec.h" - -namespace TNT -{ - -template -class Fortran_Array1D -{ - - private: - - i_refvec v_; - int n_; - T* data_; /* this normally points to v_.begin(), but - * could also point to a portion (subvector) - * of v_. - */ - - void initialize_(int n); - void copy_(T* p, const T* q, int len) const; - void set_(T* begin, T* end, const T& val); - - - public: - - typedef T value_type; - - - Fortran_Array1D(); - explicit Fortran_Array1D(int n); - Fortran_Array1D(int n, const T &a); - Fortran_Array1D(int n, T *a); - inline Fortran_Array1D(const Fortran_Array1D &A); - inline Fortran_Array1D & operator=(const T &a); - inline Fortran_Array1D & operator=(const Fortran_Array1D &A); - inline Fortran_Array1D & ref(const Fortran_Array1D &A); - Fortran_Array1D copy() const; - Fortran_Array1D & inject(const Fortran_Array1D & A); - inline T& operator()(int i); - inline const T& operator()(int i) const; - inline int dim1() const; - inline int dim() const; - ~Fortran_Array1D(); - - - /* ... extended interface ... */ - - inline int ref_count() const; - inline Fortran_Array1D subarray(int i0, int i1); - -}; - - - - -template -Fortran_Array1D::Fortran_Array1D() : v_(), n_(0), data_(0) {} - -template -Fortran_Array1D::Fortran_Array1D(const Fortran_Array1D &A) : v_(A.v_), n_(A.n_), - data_(A.data_) -{ -#ifdef TNT_DEBUG - std::cout << "Created Fortran_Array1D(const Fortran_Array1D &A) \n"; -#endif - -} - - -template -Fortran_Array1D::Fortran_Array1D(int n) : v_(n), n_(n), data_(v_.begin()) -{ -#ifdef TNT_DEBUG - std::cout << "Created Fortran_Array1D(int n) \n"; -#endif -} - -template -Fortran_Array1D::Fortran_Array1D(int n, const T &val) : v_(n), n_(n), data_(v_.begin()) -{ -#ifdef TNT_DEBUG - std::cout << "Created Fortran_Array1D(int n, const T& val) \n"; -#endif - set_(data_, data_+ n, val); - -} - -template -Fortran_Array1D::Fortran_Array1D(int n, T *a) : v_(a), n_(n) , data_(v_.begin()) -{ -#ifdef TNT_DEBUG - std::cout << "Created Fortran_Array1D(int n, T* a) \n"; -#endif -} - -template -inline T& Fortran_Array1D::operator()(int i) -{ -#ifdef TNT_BOUNDS_CHECK - assert(i>= 1); - assert(i <= n_); -#endif - return data_[i-1]; -} - -template -inline const T& Fortran_Array1D::operator()(int i) const -{ -#ifdef TNT_BOUNDS_CHECK - assert(i>= 1); - assert(i <= n_); -#endif - return data_[i-1]; -} - - - - -template -Fortran_Array1D & Fortran_Array1D::operator=(const T &a) -{ - set_(data_, data_+n_, a); - return *this; -} - -template -Fortran_Array1D Fortran_Array1D::copy() const -{ - Fortran_Array1D A( n_); - copy_(A.data_, data_, n_); - - return A; -} - - -template -Fortran_Array1D & Fortran_Array1D::inject(const Fortran_Array1D &A) -{ - if (A.n_ == n_) - copy_(data_, A.data_, n_); - - return *this; -} - - - - - -template -Fortran_Array1D & Fortran_Array1D::ref(const Fortran_Array1D &A) -{ - if (this != &A) - { - v_ = A.v_; /* operator= handles the reference counting. */ - n_ = A.n_; - data_ = A.data_; - - } - return *this; -} - -template -Fortran_Array1D & Fortran_Array1D::operator=(const Fortran_Array1D &A) -{ - return ref(A); -} - -template -inline int Fortran_Array1D::dim1() const { return n_; } - -template -inline int Fortran_Array1D::dim() const { return n_; } - -template -Fortran_Array1D::~Fortran_Array1D() {} - - -/* ............................ exented interface ......................*/ - -template -inline int Fortran_Array1D::ref_count() const -{ - return v_.ref_count(); -} - -template -inline Fortran_Array1D Fortran_Array1D::subarray(int i0, int i1) -{ -#ifdef TNT_DEBUG - std::cout << "entered subarray. \n"; -#endif - if ((i0 > 0) && (i1 < n_) || (i0 <= i1)) - { - Fortran_Array1D X(*this); /* create a new instance of this array. */ - X.n_ = i1-i0+1; - X.data_ += i0; - - return X; - } - else - { -#ifdef TNT_DEBUG - std::cout << "subarray: null return.\n"; -#endif - return Fortran_Array1D(); - } -} - - -/* private internal functions */ - - -template -void Fortran_Array1D::set_(T* begin, T* end, const T& a) -{ - for (T* p=begin; p -void Fortran_Array1D::copy_(T* p, const T* q, int len) const -{ - T *end = p + len; - while (p - -namespace TNT -{ - - -/** - Write an array to a character outstream. Output format is one that can - be read back in via the in-stream operator: one integer - denoting the array dimension (n), followed by n elements, - one per line. - -*/ -template -std::ostream& operator<<(std::ostream &s, const Fortran_Array1D &A) -{ - int N=A.dim1(); - - s << N << "\n"; - for (int j=1; j<=N; j++) - { - s << A(j) << "\n"; - } - s << "\n"; - - return s; -} - -/** - Read an array from a character stream. Input format - is one integer, denoting the dimension (n), followed - by n whitespace-separated elments. Newlines are ignored - -

- Note: the array being read into references new memory - storage. If the intent is to fill an existing conformant - array, use cin >> B; A.inject(B) ); - instead or read the elements in one-a-time by hand. - - @param s the charater to read from (typically std::in) - @param A the array to read into. -*/ -template -std::istream& operator>>(std::istream &s, Fortran_Array1D &A) -{ - int N; - s >> N; - - Fortran_Array1D B(N); - for (int i=1; i<=N; i++) - s >> B(i); - A = B; - return s; -} - - -template -Fortran_Array1D operator+(const Fortran_Array1D &A, const Fortran_Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() != n ) - return Fortran_Array1D(); - - else - { - Fortran_Array1D C(n); - - for (int i=1; i<=n; i++) - { - C(i) = A(i) + B(i); - } - return C; - } -} - - - -template -Fortran_Array1D operator-(const Fortran_Array1D &A, const Fortran_Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() != n ) - return Fortran_Array1D(); - - else - { - Fortran_Array1D C(n); - - for (int i=1; i<=n; i++) - { - C(i) = A(i) - B(i); - } - return C; - } -} - - -template -Fortran_Array1D operator*(const Fortran_Array1D &A, const Fortran_Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() != n ) - return Fortran_Array1D(); - - else - { - Fortran_Array1D C(n); - - for (int i=1; i<=n; i++) - { - C(i) = A(i) * B(i); - } - return C; - } -} - - -template -Fortran_Array1D operator/(const Fortran_Array1D &A, const Fortran_Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() != n ) - return Fortran_Array1D(); - - else - { - Fortran_Array1D C(n); - - for (int i=1; i<=n; i++) - { - C(i) = A(i) / B(i); - } - return C; - } -} - - - - - - - - - -template -Fortran_Array1D& operator+=(Fortran_Array1D &A, const Fortran_Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() == n) - { - for (int i=1; i<=n; i++) - { - A(i) += B(i); - } - } - return A; -} - - - - -template -Fortran_Array1D& operator-=(Fortran_Array1D &A, const Fortran_Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() == n) - { - for (int i=1; i<=n; i++) - { - A(i) -= B(i); - } - } - return A; -} - - - -template -Fortran_Array1D& operator*=(Fortran_Array1D &A, const Fortran_Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() == n) - { - for (int i=1; i<=n; i++) - { - A(i) *= B(i); - } - } - return A; -} - - - - -template -Fortran_Array1D& operator/=(Fortran_Array1D &A, const Fortran_Array1D &B) -{ - int n = A.dim1(); - - if (B.dim1() == n) - { - for (int i=1; i<=n; i++) - { - A(i) /= B(i); - } - } - return A; -} - - -} // namespace TNT - -#endif diff --git a/lib/include/tnt/tnt_fortran_array2d.h b/lib/include/tnt/tnt_fortran_array2d.h deleted file mode 100644 index f307536..0000000 --- a/lib/include/tnt/tnt_fortran_array2d.h +++ /dev/null @@ -1,225 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT): Two-dimensional Fortran numerical array -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - - -#ifndef TNT_FORTRAN_ARRAY2D_H -#define TNT_FORTRAN_ARRAY2D_H - -#include -#include - -#ifdef TNT_BOUNDS_CHECK -#include -#endif - -#include "tnt_i_refvec.h" - -namespace TNT -{ - -template -class Fortran_Array2D -{ - - - private: - i_refvec v_; - int m_; - int n_; - T* data_; - - - void initialize_(int n); - void copy_(T* p, const T* q, int len); - void set_(T* begin, T* end, const T& val); - - public: - - typedef T value_type; - - Fortran_Array2D(); - Fortran_Array2D(int m, int n); - Fortran_Array2D(int m, int n, T *a); - Fortran_Array2D(int m, int n, const T &a); - inline Fortran_Array2D(const Fortran_Array2D &A); - inline Fortran_Array2D & operator=(const T &a); - inline Fortran_Array2D & operator=(const Fortran_Array2D &A); - inline Fortran_Array2D & ref(const Fortran_Array2D &A); - Fortran_Array2D copy() const; - Fortran_Array2D & inject(const Fortran_Array2D & A); - inline T& operator()(int i, int j); - inline const T& operator()(int i, int j) const ; - inline int dim1() const; - inline int dim2() const; - ~Fortran_Array2D(); - - /* extended interface */ - - inline int ref_count() const; - -}; - -template -Fortran_Array2D::Fortran_Array2D() : v_(), m_(0), n_(0), data_(0) {} - - -template -Fortran_Array2D::Fortran_Array2D(const Fortran_Array2D &A) : v_(A.v_), - m_(A.m_), n_(A.n_), data_(A.data_) {} - - - -template -Fortran_Array2D::Fortran_Array2D(int m, int n) : v_(m*n), m_(m), n_(n), - data_(v_.begin()) {} - -template -Fortran_Array2D::Fortran_Array2D(int m, int n, const T &val) : - v_(m*n), m_(m), n_(n), data_(v_.begin()) -{ - set_(data_, data_+m*n, val); -} - - -template -Fortran_Array2D::Fortran_Array2D(int m, int n, T *a) : v_(a), - m_(m), n_(n), data_(v_.begin()) {} - - - - -template -inline T& Fortran_Array2D::operator()(int i, int j) -{ -#ifdef TNT_BOUNDS_CHECK - assert(i >= 1); - assert(i <= m_); - assert(j >= 1); - assert(j <= n_); -#endif - - return v_[ (j-1)*m_ + (i-1) ]; - -} - -template -inline const T& Fortran_Array2D::operator()(int i, int j) const -{ -#ifdef TNT_BOUNDS_CHECK - assert(i >= 1); - assert(i <= m_); - assert(j >= 1); - assert(j <= n_); -#endif - - return v_[ (j-1)*m_ + (i-1) ]; - -} - - -template -Fortran_Array2D & Fortran_Array2D::operator=(const T &a) -{ - set_(data_, data_+m_*n_, a); - return *this; -} - -template -Fortran_Array2D Fortran_Array2D::copy() const -{ - - Fortran_Array2D B(m_,n_); - - B.inject(*this); - return B; -} - - -template -Fortran_Array2D & Fortran_Array2D::inject(const Fortran_Array2D &A) -{ - if (m_ == A.m_ && n_ == A.n_) - copy_(data_, A.data_, m_*n_); - - return *this; -} - - - -template -Fortran_Array2D & Fortran_Array2D::ref(const Fortran_Array2D &A) -{ - if (this != &A) - { - v_ = A.v_; - m_ = A.m_; - n_ = A.n_; - data_ = A.data_; - } - return *this; -} - -template -Fortran_Array2D & Fortran_Array2D::operator=(const Fortran_Array2D &A) -{ - return ref(A); -} - -template -inline int Fortran_Array2D::dim1() const { return m_; } - -template -inline int Fortran_Array2D::dim2() const { return n_; } - - -template -Fortran_Array2D::~Fortran_Array2D() -{ -} - -template -inline int Fortran_Array2D::ref_count() const { return v_.ref_count(); } - - - - -template -void Fortran_Array2D::set_(T* begin, T* end, const T& a) -{ - for (T* p=begin; p -void Fortran_Array2D::copy_(T* p, const T* q, int len) -{ - T *end = p + len; - while (p - -namespace TNT -{ - - -template -std::ostream& operator<<(std::ostream &s, const Fortran_Array2D &A) -{ - int M=A.dim1(); - int N=A.dim2(); - - s << M << " " << N << "\n"; - - for (int i=1; i<=M; i++) - { - for (int j=1; j<=N; j++) - { - s << A(i,j) << " "; - } - s << "\n"; - } - - - return s; -} - -template -std::istream& operator>>(std::istream &s, Fortran_Array2D &A) -{ - - int M, N; - - s >> M >> N; - - Fortran_Array2D B(M,N); - - for (int i=1; i<=M; i++) - for (int j=1; j<=N; j++) - { - s >> B(i,j); - } - - A = B; - return s; -} - - - - -template -Fortran_Array2D operator+(const Fortran_Array2D &A, const Fortran_Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() != m || B.dim2() != n ) - return Fortran_Array2D(); - - else - { - Fortran_Array2D C(m,n); - - for (int i=1; i<=m; i++) - { - for (int j=1; j<=n; j++) - C(i,j) = A(i,j) + B(i,j); - } - return C; - } -} - -template -Fortran_Array2D operator-(const Fortran_Array2D &A, const Fortran_Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() != m || B.dim2() != n ) - return Fortran_Array2D(); - - else - { - Fortran_Array2D C(m,n); - - for (int i=1; i<=m; i++) - { - for (int j=1; j<=n; j++) - C(i,j) = A(i,j) - B(i,j); - } - return C; - } -} - - -template -Fortran_Array2D operator*(const Fortran_Array2D &A, const Fortran_Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() != m || B.dim2() != n ) - return Fortran_Array2D(); - - else - { - Fortran_Array2D C(m,n); - - for (int i=1; i<=m; i++) - { - for (int j=1; j<=n; j++) - C(i,j) = A(i,j) * B(i,j); - } - return C; - } -} - - -template -Fortran_Array2D operator/(const Fortran_Array2D &A, const Fortran_Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() != m || B.dim2() != n ) - return Fortran_Array2D(); - - else - { - Fortran_Array2D C(m,n); - - for (int i=1; i<=m; i++) - { - for (int j=1; j<=n; j++) - C(i,j) = A(i,j) / B(i,j); - } - return C; - } -} - - - -template -Fortran_Array2D& operator+=(Fortran_Array2D &A, const Fortran_Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() == m || B.dim2() == n ) - { - for (int i=1; i<=m; i++) - { - for (int j=1; j<=n; j++) - A(i,j) += B(i,j); - } - } - return A; -} - -template -Fortran_Array2D& operator-=(Fortran_Array2D &A, const Fortran_Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() == m || B.dim2() == n ) - { - for (int i=1; i<=m; i++) - { - for (int j=1; j<=n; j++) - A(i,j) -= B(i,j); - } - } - return A; -} - -template -Fortran_Array2D& operator*=(Fortran_Array2D &A, const Fortran_Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() == m || B.dim2() == n ) - { - for (int i=1; i<=m; i++) - { - for (int j=1; j<=n; j++) - A(i,j) *= B(i,j); - } - } - return A; -} - -template -Fortran_Array2D& operator/=(Fortran_Array2D &A, const Fortran_Array2D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - - if (B.dim1() == m || B.dim2() == n ) - { - for (int i=1; i<=m; i++) - { - for (int j=1; j<=n; j++) - A(i,j) /= B(i,j); - } - } - return A; -} - -} // namespace TNT - -#endif diff --git a/lib/include/tnt/tnt_fortran_array3d.h b/lib/include/tnt/tnt_fortran_array3d.h deleted file mode 100644 index e51affb..0000000 --- a/lib/include/tnt/tnt_fortran_array3d.h +++ /dev/null @@ -1,223 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT): Three-dimensional Fortran numerical array -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - - -#ifndef TNT_FORTRAN_ARRAY3D_H -#define TNT_FORTRAN_ARRAY3D_H - -#include -#include -#ifdef TNT_BOUNDS_CHECK -#include -#endif -#include "tnt_i_refvec.h" - -namespace TNT -{ - -template -class Fortran_Array3D -{ - - - private: - - - i_refvec v_; - int m_; - int n_; - int k_; - T* data_; - - public: - - typedef T value_type; - - Fortran_Array3D(); - Fortran_Array3D(int m, int n, int k); - Fortran_Array3D(int m, int n, int k, T *a); - Fortran_Array3D(int m, int n, int k, const T &a); - inline Fortran_Array3D(const Fortran_Array3D &A); - inline Fortran_Array3D & operator=(const T &a); - inline Fortran_Array3D & operator=(const Fortran_Array3D &A); - inline Fortran_Array3D & ref(const Fortran_Array3D &A); - Fortran_Array3D copy() const; - Fortran_Array3D & inject(const Fortran_Array3D & A); - inline T& operator()(int i, int j, int k); - inline const T& operator()(int i, int j, int k) const ; - inline int dim1() const; - inline int dim2() const; - inline int dim3() const; - inline int ref_count() const; - ~Fortran_Array3D(); - - -}; - -template -Fortran_Array3D::Fortran_Array3D() : v_(), m_(0), n_(0), k_(0), data_(0) {} - - -template -Fortran_Array3D::Fortran_Array3D(const Fortran_Array3D &A) : - v_(A.v_), m_(A.m_), n_(A.n_), k_(A.k_), data_(A.data_) {} - - - -template -Fortran_Array3D::Fortran_Array3D(int m, int n, int k) : - v_(m*n*k), m_(m), n_(n), k_(k), data_(v_.begin()) {} - - - -template -Fortran_Array3D::Fortran_Array3D(int m, int n, int k, const T &val) : - v_(m*n*k), m_(m), n_(n), k_(k), data_(v_.begin()) -{ - for (T* p = data_; p < data_ + m*n*k; p++) - *p = val; -} - -template -Fortran_Array3D::Fortran_Array3D(int m, int n, int k, T *a) : - v_(a), m_(m), n_(n), k_(k), data_(v_.begin()) {} - - - - -template -inline T& Fortran_Array3D::operator()(int i, int j, int k) -{ -#ifdef TNT_BOUNDS_CHECK - assert(i >= 1); - assert(i <= m_); - assert(j >= 1); - assert(j <= n_); - assert(k >= 1); - assert(k <= k_); -#endif - - return data_[(k-1)*m_*n_ + (j-1) * m_ + i-1]; - -} - -template -inline const T& Fortran_Array3D::operator()(int i, int j, int k) const -{ -#ifdef TNT_BOUNDS_CHECK - assert(i >= 1); - assert(i <= m_); - assert(j >= 1); - assert(j <= n_); - assert(k >= 1); - assert(k <= k_); -#endif - - return data_[(k-1)*m_*n_ + (j-1) * m_ + i-1]; -} - - -template -Fortran_Array3D & Fortran_Array3D::operator=(const T &a) -{ - - T *end = data_ + m_*n_*k_; - - for (T *p=data_; p != end; *p++ = a); - - return *this; -} - -template -Fortran_Array3D Fortran_Array3D::copy() const -{ - - Fortran_Array3D B(m_, n_, k_); - B.inject(*this); - return B; - -} - - -template -Fortran_Array3D & Fortran_Array3D::inject(const Fortran_Array3D &A) -{ - - if (m_ == A.m_ && n_ == A.n_ && k_ == A.k_) - { - T *p = data_; - T *end = data_ + m_*n_*k_; - const T* q = A.data_; - for (; p < end; *p++ = *q++); - } - return *this; -} - - - - -template -Fortran_Array3D & Fortran_Array3D::ref(const Fortran_Array3D &A) -{ - - if (this != &A) - { - v_ = A.v_; - m_ = A.m_; - n_ = A.n_; - k_ = A.k_; - data_ = A.data_; - } - return *this; -} - -template -Fortran_Array3D & Fortran_Array3D::operator=(const Fortran_Array3D &A) -{ - return ref(A); -} - -template -inline int Fortran_Array3D::dim1() const { return m_; } - -template -inline int Fortran_Array3D::dim2() const { return n_; } - -template -inline int Fortran_Array3D::dim3() const { return k_; } - - -template -inline int Fortran_Array3D::ref_count() const -{ - return v_.ref_count(); -} - -template -Fortran_Array3D::~Fortran_Array3D() -{ -} - - -} /* namespace TNT */ - -#endif -/* TNT_FORTRAN_ARRAY3D_H */ - diff --git a/lib/include/tnt/tnt_fortran_array3d_utils.h b/lib/include/tnt/tnt_fortran_array3d_utils.h deleted file mode 100644 index a13a275..0000000 --- a/lib/include/tnt/tnt_fortran_array3d_utils.h +++ /dev/null @@ -1,249 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT) -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - -#ifndef TNT_FORTRAN_ARRAY3D_UTILS_H -#define TNT_FORTRAN_ARRAY3D_UTILS_H - -#include -#include - -namespace TNT -{ - - -template -std::ostream& operator<<(std::ostream &s, const Fortran_Array3D &A) -{ - int M=A.dim1(); - int N=A.dim2(); - int K=A.dim3(); - - s << M << " " << N << " " << K << "\n"; - - for (int i=1; i<=M; i++) - { - for (int j=1; j<=N; j++) - { - for (int k=1; k<=K; k++) - s << A(i,j,k) << " "; - s << "\n"; - } - s << "\n"; - } - - - return s; -} - -template -std::istream& operator>>(std::istream &s, Fortran_Array3D &A) -{ - - int M, N, K; - - s >> M >> N >> K; - - Fortran_Array3D B(M,N,K); - - for (int i=1; i<=M; i++) - for (int j=1; j<=N; j++) - for (int k=1; k<=K; k++) - s >> B(i,j,k); - - A = B; - return s; -} - - -template -Fortran_Array3D operator+(const Fortran_Array3D &A, const Fortran_Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) - return Fortran_Array3D(); - - else - { - Fortran_Array3D C(m,n,p); - - for (int i=1; i<=m; i++) - for (int j=1; j<=n; j++) - for (int k=1; k<=p; k++) - C(i,j,k) = A(i,j,k)+ B(i,j,k); - - return C; - } -} - - -template -Fortran_Array3D operator-(const Fortran_Array3D &A, const Fortran_Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) - return Fortran_Array3D(); - - else - { - Fortran_Array3D C(m,n,p); - - for (int i=1; i<=m; i++) - for (int j=1; j<=n; j++) - for (int k=1; k<=p; k++) - C(i,j,k) = A(i,j,k)- B(i,j,k); - - return C; - } -} - - -template -Fortran_Array3D operator*(const Fortran_Array3D &A, const Fortran_Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) - return Fortran_Array3D(); - - else - { - Fortran_Array3D C(m,n,p); - - for (int i=1; i<=m; i++) - for (int j=1; j<=n; j++) - for (int k=1; k<=p; k++) - C(i,j,k) = A(i,j,k)* B(i,j,k); - - return C; - } -} - - -template -Fortran_Array3D operator/(const Fortran_Array3D &A, const Fortran_Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) - return Fortran_Array3D(); - - else - { - Fortran_Array3D C(m,n,p); - - for (int i=1; i<=m; i++) - for (int j=1; j<=n; j++) - for (int k=1; k<=p; k++) - C(i,j,k) = A(i,j,k)/ B(i,j,k); - - return C; - } -} - - -template -Fortran_Array3D& operator+=(Fortran_Array3D &A, const Fortran_Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) - { - for (int i=1; i<=m; i++) - for (int j=1; j<=n; j++) - for (int k=1; k<=p; k++) - A(i,j,k) += B(i,j,k); - } - - return A; -} - - -template -Fortran_Array3D& operator-=(Fortran_Array3D &A, const Fortran_Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) - { - for (int i=1; i<=m; i++) - for (int j=1; j<=n; j++) - for (int k=1; k<=p; k++) - A(i,j,k) -= B(i,j,k); - } - - return A; -} - - -template -Fortran_Array3D& operator*=(Fortran_Array3D &A, const Fortran_Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) - { - for (int i=1; i<=m; i++) - for (int j=1; j<=n; j++) - for (int k=1; k<=p; k++) - A(i,j,k) *= B(i,j,k); - } - - return A; -} - - -template -Fortran_Array3D& operator/=(Fortran_Array3D &A, const Fortran_Array3D &B) -{ - int m = A.dim1(); - int n = A.dim2(); - int p = A.dim3(); - - if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) - { - for (int i=1; i<=m; i++) - for (int j=1; j<=n; j++) - for (int k=1; k<=p; k++) - A(i,j,k) /= B(i,j,k); - } - - return A; -} - - -} // namespace TNT - -#endif diff --git a/lib/include/tnt/tnt_i_refvec.h b/lib/include/tnt/tnt_i_refvec.h deleted file mode 100644 index 5a67eb5..0000000 --- a/lib/include/tnt/tnt_i_refvec.h +++ /dev/null @@ -1,243 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT) -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - - -#ifndef TNT_I_REFVEC_H -#define TNT_I_REFVEC_H - -#include -#include - -#ifdef TNT_BOUNDS_CHECK -#include -#endif - -#ifndef NULL -#define NULL 0 -#endif - -namespace TNT -{ -/* - Internal representation of ref-counted array. The TNT - arrays all use this building block. - -

- If an array block is created by TNT, then every time - an assignment is made, the left-hand-side reference - is decreased by one, and the right-hand-side refernce - count is increased by one. If the array block was - external to TNT, the refernce count is a NULL pointer - regardless of how many references are made, since the - memory is not freed by TNT. - - - -*/ -template -class i_refvec -{ - - - private: - T* data_; - int *ref_count_; - - - public: - - i_refvec(); - explicit i_refvec(int n); - inline i_refvec(T* data); - inline i_refvec(const i_refvec &v); - inline T* begin(); - inline const T* begin() const; - inline T& operator[](int i); - inline const T& operator[](int i) const; - inline i_refvec & operator=(const i_refvec &V); - void copy_(T* p, const T* q, const T* e); - void set_(T* p, const T* b, const T* e); - inline int ref_count() const; - inline int is_null() const; - inline void destroy(); - ~i_refvec(); - -}; - -template -void i_refvec::copy_(T* p, const T* q, const T* e) -{ - for (T* t=p; q -i_refvec::i_refvec() : data_(NULL), ref_count_(NULL) {} - -/** - In case n is 0 or negative, it does NOT call new. -*/ -template -i_refvec::i_refvec(int n) : data_(NULL), ref_count_(NULL) -{ - if (n >= 1) - { -#ifdef TNT_DEBUG - std::cout << "new data storage.\n"; -#endif - data_ = new T[n]; - ref_count_ = new int; - *ref_count_ = 1; - } -} - -template -inline i_refvec::i_refvec(const i_refvec &V): data_(V.data_), - ref_count_(V.ref_count_) -{ - if (V.ref_count_ != NULL) - (*(V.ref_count_))++; -} - - -template -i_refvec::i_refvec(T* data) : data_(data), ref_count_(NULL) {} - -template -inline T* i_refvec::begin() -{ - return data_; -} - -template -inline const T& i_refvec::operator[](int i) const -{ - return data_[i]; -} - -template -inline T& i_refvec::operator[](int i) -{ - return data_[i]; -} - - -template -inline const T* i_refvec::begin() const -{ - return data_; -} - - - -template -i_refvec & i_refvec::operator=(const i_refvec &V) -{ - if (this == &V) - return *this; - - - if (ref_count_ != NULL) - { - (*ref_count_) --; - if ((*ref_count_) == 0) - destroy(); - } - - data_ = V.data_; - ref_count_ = V.ref_count_; - - if (V.ref_count_ != NULL) - (*(V.ref_count_))++; - - return *this; -} - -template -void i_refvec::destroy() -{ - if (ref_count_ != NULL) - { -#ifdef TNT_DEBUG - std::cout << "destorying data... \n"; -#endif - delete ref_count_; - -#ifdef TNT_DEBUG - std::cout << "deleted ref_count_ ...\n"; -#endif - if (data_ != NULL) - delete []data_; -#ifdef TNT_DEBUG - std::cout << "deleted data_[] ...\n"; -#endif - data_ = NULL; - } -} - -/* -* return 1 is vector is empty, 0 otherwise -* -* if is_null() is false and ref_count() is 0, then -* -*/ -template -int i_refvec::is_null() const -{ - return (data_ == NULL ? 1 : 0); -} - -/* -* returns -1 if data is external, -* returns 0 if a is NULL array, -* otherwise returns the positive number of vectors sharing -* this data space. -*/ -template -int i_refvec::ref_count() const -{ - if (data_ == NULL) - return 0; - else - return (ref_count_ != NULL ? *ref_count_ : -1) ; -} - -template -i_refvec::~i_refvec() -{ - if (ref_count_ != NULL) - { - (*ref_count_)--; - - if (*ref_count_ == 0) - destroy(); - } -} - - -} /* namespace TNT */ - - - - - -#endif -/* TNT_I_REFVEC_H */ - diff --git a/lib/include/tnt/tnt_math_utils.h b/lib/include/tnt/tnt_math_utils.h deleted file mode 100644 index f9c1c91..0000000 --- a/lib/include/tnt/tnt_math_utils.h +++ /dev/null @@ -1,34 +0,0 @@ -#ifndef MATH_UTILS_H -#define MATH_UTILS_H - -/* needed for fabs, sqrt() below */ -#include - - - -namespace TNT -{ -/** - @returns hypotenuse of real (non-complex) scalars a and b by - avoiding underflow/overflow - using (a * sqrt( 1 + (b/a) * (b/a))), rather than - sqrt(a*a + b*b). -*/ -template -Real hypot(const Real &a, const Real &b) -{ - - if (a== 0) - return abs(b); - else - { - Real c = b/a; - return fabs(a) * sqrt(1 + c*c); - } -} -} /* TNT namespace */ - - - -#endif -/* MATH_UTILS_H */ diff --git a/lib/include/tnt/tnt_sparse_matrix_csr.h b/lib/include/tnt/tnt_sparse_matrix_csr.h deleted file mode 100644 index 0d4fde1..0000000 --- a/lib/include/tnt/tnt_sparse_matrix_csr.h +++ /dev/null @@ -1,103 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT) -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - -#ifndef TNT_SPARSE_MATRIX_CSR_H -#define TNT_SPARSE_MATRIX_CSR_H - -#include "tnt_array1d.h" - -namespace TNT -{ - - -/** - Read-only view of a sparse matrix in compressed-row storage - format. Neither array elements (nonzeros) nor sparsity - structure can be modified. If modifications are required, - create a new view. - -

- Index values begin at 0. - -

- Storage requirements: An (m x n) matrix with - nz nonzeros requires no more than ((T+I)*nz + M*I) - bytes, where T is the size of data elements and - I is the size of integers. - - -*/ -template -class Sparse_Matrix_CompRow { - -private: - Array1D val_; // data values (nz_ elements) - Array1D rowptr_; // row_ptr (dim_[0]+1 elements) - Array1D colind_; // col_ind (nz_ elements) - - int dim1_; // number of rows - int dim2_; // number of cols - -public: - - 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); - - - - inline const T& val(int i) const { return val_[i]; } - inline const int& row_ptr(int i) const { return rowptr_[i]; } - inline const int& col_ind(int i) const { return colind_[i];} - - inline int dim1() const {return dim1_;} - inline int dim2() const {return dim2_;} - int NumNonzeros() const {return val_.dim1();} - - - Sparse_Matrix_CompRow& operator=( - const Sparse_Matrix_CompRow &R); - - - -}; - -/** - Construct a read-only view of existing sparse matrix in - compressed-row storage format. - - @param M the number of rows of sparse matrix - @param N the number of columns of sparse matrix - @param nz the number of nonzeros - @param val a contiguous list of nonzero values - @param r row-pointers: r[i] denotes the begining position of row i - (i.e. the ith row begins at val[row[i]]). - @param c column-indices: c[i] denotes the column location of val[i] -*/ -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), - rowptr_(M, r), colind_(nz, c), dim1_(M), dim2_(N) {} - - -} -// namespace TNT - -#endif diff --git a/lib/include/tnt/tnt_stopwatch.h b/lib/include/tnt/tnt_stopwatch.h deleted file mode 100644 index 8dc5d23..0000000 --- a/lib/include/tnt/tnt_stopwatch.h +++ /dev/null @@ -1,95 +0,0 @@ -/* -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - - -#ifndef STOPWATCH_H -#define STOPWATCH_H - -// for clock() and CLOCKS_PER_SEC -#include - - -namespace TNT -{ - -inline static double seconds(void) -{ - const double secs_per_tick = 1.0 / CLOCKS_PER_SEC; - return ( (double) clock() ) * secs_per_tick; -} - -class Stopwatch { - private: - int running_; - double start_time_; - double total_; - - public: - inline Stopwatch(); - inline void start(); - inline double stop(); - inline double read(); - inline void resume(); - inline int running(); -}; - -inline Stopwatch::Stopwatch() : running_(0), start_time_(0.0), total_(0.0) {} - -void Stopwatch::start() -{ - running_ = 1; - total_ = 0.0; - start_time_ = seconds(); -} - -double Stopwatch::stop() -{ - if (running_) - { - total_ += (seconds() - start_time_); - running_ = 0; - } - return total_; -} - -inline void Stopwatch::resume() -{ - if (!running_) - { - start_time_ = seconds(); - running_ = 1; - } -} - - -inline double Stopwatch::read() -{ - if (running_) - { - stop(); - resume(); - } - return total_; -} - - -} /* TNT namespace */ -#endif - - - diff --git a/lib/include/tnt/tnt_subscript.h b/lib/include/tnt/tnt_subscript.h deleted file mode 100644 index d8fe120..0000000 --- a/lib/include/tnt/tnt_subscript.h +++ /dev/null @@ -1,54 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT) -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - -#ifndef TNT_SUBSCRPT_H -#define TNT_SUBSCRPT_H - - -//--------------------------------------------------------------------- -// This definition describes the default TNT data type used for -// indexing into TNT matrices and vectors. The data type should -// be wide enough to index into large arrays. It defaults to an -// "int", but can be overriden at compile time redefining TNT_SUBSCRIPT_TYPE, -// e.g. -// -// c++ -DTNT_SUBSCRIPT_TYPE='unsigned int' ... -// -//--------------------------------------------------------------------- -// - -#ifndef TNT_SUBSCRIPT_TYPE -#define TNT_SUBSCRIPT_TYPE int -#endif - -namespace TNT -{ - typedef TNT_SUBSCRIPT_TYPE Subscript; -} /* namespace TNT */ - - -// () indexing in TNT means 1-offset, i.e. x(1) and A(1,1) are the -// first elements. This offset is left as a macro for future -// purposes, but should not be changed in the current release. -// -// -#define TNT_BASE_OFFSET (1) - -#endif diff --git a/lib/include/tnt/tnt_vec.h b/lib/include/tnt/tnt_vec.h deleted file mode 100644 index 3455d79..0000000 --- a/lib/include/tnt/tnt_vec.h +++ /dev/null @@ -1,404 +0,0 @@ -/* -* -* Template Numerical Toolkit (TNT) -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -*/ - - - -#ifndef TNT_VEC_H -#define TNT_VEC_H - -#include "tnt_subscript.h" -#include -#include -#include -#include - -namespace TNT -{ - -/** - [Deprecatred] Value-based vector class from pre-1.0 - TNT version. Kept here for backward compatiblity, but should - use the newer TNT::Array1D classes instead. - -*/ - -template -class Vector -{ - - - public: - - typedef Subscript size_type; - typedef T value_type; - typedef T element_type; - typedef T* pointer; - typedef T* iterator; - typedef T& reference; - typedef const T* const_iterator; - typedef const T& const_reference; - - Subscript lbound() const { return 1;} - - protected: - T* v_; - T* vm1_; // pointer adjustment for optimzied 1-offset indexing - Subscript n_; - - // internal helper function to create the array - // of row pointers - - void initialize(Subscript N) - { - // adjust pointers so that they are 1-offset: - // v_[] is the internal contiguous array, it is still 0-offset - // - assert(v_ == NULL); - v_ = new T[N]; - assert(v_ != NULL); - vm1_ = v_-1; - n_ = N; - } - - void copy(const T* v) - { - Subscript N = n_; - Subscript i; - -#ifdef TNT_UNROLL_LOOPS - Subscript Nmod4 = N & 3; - Subscript N4 = N - Nmod4; - - for (i=0; i &A) : v_(0), vm1_(0), n_(0) - { - initialize(A.n_); - copy(A.v_); - } - - Vector(Subscript N, const T& value = T()) : v_(0), vm1_(0), n_(0) - { - initialize(N); - set(value); - } - - Vector(Subscript N, const T* v) : v_(0), vm1_(0), n_(0) - { - initialize(N); - copy(v); - } - - Vector(Subscript N, char *s) : v_(0), vm1_(0), n_(0) - { - initialize(N); - std::istringstream ins(s); - - Subscript i; - - for (i=0; i> v_[i]; - } - - - // methods - // - Vector& newsize(Subscript N) - { - if (n_ == N) return *this; - - destroy(); - initialize(N); - - return *this; - } - - - // assignments - // - Vector& operator=(const Vector &A) - { - if (v_ == A.v_) - return *this; - - if (n_ == A.n_) // no need to re-alloc - copy(A.v_); - - else - { - destroy(); - initialize(A.n_); - copy(A.v_); - } - - return *this; - } - - Vector& operator=(const T& scalar) - { - set(scalar); - return *this; - } - - inline Subscript dim() const - { - return n_; - } - - inline Subscript size() const - { - return n_; - } - - - inline reference operator()(Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= n_) ; -#endif - return vm1_[i]; - } - - inline const_reference operator() (Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= n_) ; -#endif - return vm1_[i]; - } - - inline reference operator[](Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i); - assert(i < n_) ; -#endif - return v_[i]; - } - - inline const_reference operator[](Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i); - - - - - - - assert(i < n_) ; -#endif - return v_[i]; - } - - - -}; - - -/* *************************** I/O ********************************/ - -template -std::ostream& operator<<(std::ostream &s, const Vector &A) -{ - Subscript N=A.dim(); - - s << N << "\n"; - - for (Subscript i=0; i -std::istream & operator>>(std::istream &s, Vector &A) -{ - - Subscript N; - - s >> N; - - if ( !(N == A.size() )) - { - A.newsize(N); - } - - - for (Subscript i=0; i> A[i]; - - - return s; -} - -// *******************[ basic matrix algorithms ]*************************** - - -template -Vector operator+(const Vector &A, - const Vector &B) -{ - Subscript N = A.dim(); - - assert(N==B.dim()); - - Vector tmp(N); - Subscript i; - - for (i=0; i -Vector operator-(const Vector &A, - const Vector &B) -{ - Subscript N = A.dim(); - - assert(N==B.dim()); - - Vector tmp(N); - Subscript i; - - for (i=0; i -Vector operator*(const Vector &A, - const Vector &B) -{ - Subscript N = A.dim(); - - assert(N==B.dim()); - - Vector tmp(N); - Subscript i; - - for (i=0; i -T dot_prod(const Vector &A, const Vector &B) -{ - Subscript N = A.dim(); - assert(N == B.dim()); - - Subscript i; - T sum = 0; - - for (i=0; i