From a0c8925b76d8d3fbd24dc4b9b3c84a8f1b823d4b Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
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 -
26 files changed, 7184 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
(limited to 'lib/include')
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;
-
-/**
- <P>
- 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.
-
- <p>Typical usage looks like:
- <pre>
- Array2D<double> A(n,n);
- Array2D<double> L;
-
- ...
-
- Cholesky<double> chol(A);
-
- if (chol.is_spd())
- L = chol.getL();
-
- else
- cout << "factorization was not complete.\n";
-
- </pre>
-
-
- <p>
- (Adapted from JAMA, a Java Matrix Library, developed by jointly
- by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
-
- */
-
-template <class Real>
-class Cholesky
-{
- Array2D<Real> L_; // lower triangular factor
- int isspd; // 1 if matrix to be factored was SPD
-
-public:
-
- Cholesky();
- Cholesky(const Array2D<Real> &A);
- Array2D<Real> getL() const;
- Array1D<Real> solve(const Array1D<Real> &B);
- Array2D<Real> solve(const Array2D<Real> &B);
- int is_spd() const;
-
-};
-
-template <class Real>
-Cholesky<Real>::Cholesky() : L_(0,0), isspd(0) {}
-
-/**
- @return 1, if original matrix to be factored was symmetric
- positive-definite (SPD).
-*/
-template <class Real>
-int Cholesky<Real>::is_spd() const
-{
- return isspd;
-}
-
-/**
- @return the lower triangular factor, L, such that L*L'=A.
-*/
-template <class Real>
-Array2D<Real> Cholesky<Real>::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 <class Real>
-Cholesky<Real>::Cholesky(const Array2D<Real> &A)
-{
-
-
- int m = A.dim1();
- int n = A.dim2();
-
- isspd = (m == n);
-
- if (m != n)
- {
- L_ = Array2D<Real>(0,0);
- return;
- }
-
- L_ = Array2D<Real>(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 <class Real>
-Array1D<Real> Cholesky<Real>::solve(const Array1D<Real> &b)
-{
- int n = L_.dim1();
- if (b.dim1() != n)
- return Array1D<Real>();
-
-
- Array1D<Real> 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 <class Real>
-Array2D<Real> Cholesky<Real>::solve(const Array2D<Real> &B)
-{
- int n = L_.dim1();
- if (B.dim1() != n)
- return Array2D<Real>();
-
-
- Array2D<Real> 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<nx; j++)
- {
- for (int k = n-1; k >= 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 <algorithm>
-// for min(), max() below
-
-#include <cmath>
-// 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.
-<P>
- 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.
-
-<P>
- 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
-<pre>
-
- u + iv . . . . .
- . u - iv . . . .
- . . a + ib . . .
- . . . a - ib . .
- . . . . x .
- . . . . . y
-</pre>
- then D looks like
-<pre>
-
- u v . . . .
- -v u . . . .
- . . a b . .
- . . -b a . .
- . . . . x .
- . . . . . y
-</pre>
- This keeps V a real matrix in both symmetric and non-symmetric
- cases, and A*V = V*D.
-
-
-
- <p>
- 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.
-
- <p>
- (Adapted from JAMA, a Java Matrix Library, developed by jointly
- by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
-**/
-
-template <class Real>
-class Eigenvalue
-{
-
-
- /** Row and column dimension (square matrix). */
- int n;
-
- int issymmetric; /* boolean*/
-
- /** Arrays for internal storage of eigenvalues. */
-
- TNT::Array1D<Real> d; /* real part */
- TNT::Array1D<Real> e; /* img part */
-
- /** Array for internal storage of eigenvectors. */
- TNT::Array2D<Real> V;
-
- /** Array for internal storage of nonsymmetric Hessenberg form.
- @serial internal storage of nonsymmetric Hessenberg form.
- */
- TNT::Array2D<Real> H;
-
-
- /** Working storage for nonsymmetric algorithm.
- @serial working storage for nonsymmetric algorithm.
- */
- TNT::Array1D<Real> 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<Real> &A) {
- n = A.dim2();
- V = Array2D<Real>(n,n);
- d = Array1D<Real>(n);
- e = Array1D<Real>(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<Real>(n,n);
- ort = TNT::Array1D<Real>(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<Real> &V_) {
- V_ = V;
- return;
- }
-
- /** Return the real parts of the eigenvalues
- @return real(diag(D))
- */
-
- void getRealEigenvalues (TNT::Array1D<Real> &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<Real> &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
-<pre>
-
- u + iv . . . . .
- . u - iv . . . .
- . . a + ib . . .
- . . . a - ib . .
- . . . . x .
- . . . . . y
-</pre>
- then D looks like
-<pre>
-
- u v . . . .
- -v u . . . .
- . . a b . .
- . . -b a . .
- . . . . x .
- . . . . . y
-</pre>
- 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<Real> &D) {
- D = Array2D<Real>(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 <algorithm>
-//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.
- <P>
- 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.
- <P>
- 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 Real>
-class LU
-{
-
-
-
- /* Array for internal storage of decomposition. */
- Array2D<Real> LU_;
- int m, n, pivsign;
- Array1D<int> piv;
-
-
- Array2D<Real> permute_copy(const Array2D<Real> &A,
- const Array1D<int> &piv, int j0, int j1)
- {
- int piv_length = piv.dim();
-
- Array2D<Real> 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<Real> permute_copy(const Array1D<Real> &A,
- const Array1D<int> &piv)
- {
- int piv_length = piv.dim();
- if (piv_length != A.dim())
- return Array1D<Real>();
-
- Array1D<Real> 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<Real> &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<Real> 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<Real> getL () {
- Array2D<Real> 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<Real> getU () {
- Array2D<Real> 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<int> 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<Real> solve (const Array2D<Real> &B)
- {
-
- /* Dimensions: A is mxn, X is nxk, B is mxk */
-
- if (B.dim1() != m) {
- return Array2D<Real>(0,0);
- }
- if (!isNonsingular()) {
- return Array2D<Real>(0,0);
- }
-
- // Copy right hand side with pivoting
- int nx = B.dim2();
-
-
- Array2D<Real> 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<Real> solve (const Array1D<Real> &b)
- {
-
- /* Dimensions: A is mxn, X is nxk, B is mxk */
-
- if (b.dim1() != m) {
- return Array1D<Real>();
- }
- if (!isNonsingular()) {
- return Array1D<Real>();
- }
-
-
- Array1D<Real> 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
-{
-
-/**
-<p>
- 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.
-<P>
- 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).
-
-<p>
- 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.
-
- <p>
- (Adapted from JAMA, a Java Matrix Library, developed by jointly
- by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
-*/
-
-template <class Real>
-class QR {
-
-
- /** Array for internal storage of decomposition.
- @serial internal array storage.
- */
-
- TNT::Array2D<Real> 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<Real> Rdiag;
-
-
-public:
-
-/**
- Create a QR factorization object for A.
-
- @param A rectangular (m>=n) matrix.
-*/
- QR(const TNT::Array2D<Real> &A) /* constructor */
- {
- QR_ = A.copy();
- m = A.dim1();
- n = A.dim2();
- Rdiag = TNT::Array1D<Real>(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<Real> getHouseholder (void) const
- {
- TNT::Array2D<Real> 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<Real> getR() const
- {
- TNT::Array2D<Real> 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<Real> getQ() const
- {
- int i=0, j=0, k=0;
-
- TNT::Array2D<Real> 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<Real> solve(const TNT::Array1D<Real> &b) const
- {
- if (b.dim1() != m) /* arrays must be conformant */
- return TNT::Array1D<Real>();
-
- if ( !isFullRank() ) /* matrix is rank deficient */
- {
- return TNT::Array1D<Real>();
- }
-
- TNT::Array1D<Real> 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<Real> x_(n);
- for (int i=0; i<n; i++)
- x_[i] = x[i];
-
- return x_;
- }
-
- /** Least squares solution of A*X = B
- @param B m x k Array (must conform).
- @return X n x k Array 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 (0x0) array.
- */
-
- TNT::Array2D<Real> solve(const TNT::Array2D<Real> &B) const
- {
- if (B.dim1() != m) /* arrays must be conformant */
- return TNT::Array2D<Real>(0,0);
-
- if ( !isFullRank() ) /* matrix is rank deficient */
- {
- return TNT::Array2D<Real>(0,0);
- }
-
- int nx = B.dim2();
- TNT::Array2D<Real> 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<Real> X_(n,nx);
- for (i=0; i<n; i++)
- for (j=0; j<nx; j++)
- X_[i][j] = X[i][j];
-
- return X_;
- }
-
-
-};
-
-
-}
-// namespace JAMA
-
-#endif
-// JAMA_QR__H
-
diff --git a/lib/include/tnt/jama_svd.h b/lib/include/tnt/jama_svd.h
deleted file mode 100644
index 72ce3a7..0000000
--- a/lib/include/tnt/jama_svd.h
+++ /dev/null
@@ -1,543 +0,0 @@
-#ifndef JAMA_SVD_H
-#define JAMA_SVD_H
-
-
-#include "tnt_array1d.h"
-#include "tnt_array1d_utils.h"
-#include "tnt_array2d.h"
-#include "tnt_array2d_utils.h"
-#include "tnt_math_utils.h"
-
-#include <algorithm>
-// for min(), max() below
-#include <cmath>
-// 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.
- <P>
- 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'.
- <P>
- The singular values, sigma[k] = S[k][k], are ordered so that
- sigma[0] >= sigma[1] >= ... >= sigma[n-1].
- <P>
- 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.
-
- <p>
- (Adapted from JAMA, a Java Matrix Library, developed by jointly
- by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
- */
-template <class Real>
-class SVD
-{
-
-
- Array2D<Real> U, V;
- Array1D<Real> s;
- int m, n;
-
- public:
-
-
- SVD (const Array2D<Real> &Arg) {
-
-
- m = Arg.dim1();
- n = Arg.dim2();
- int nu = std::min(m,n);
- s = Array1D<Real>(std::min(m+1,n));
- U = Array2D<Real>(m, nu, Real(0));
- V = Array2D<Real>(n,n);
- Array1D<Real> e(n);
- Array1D<Real> work(m);
- Array2D<Real> A(Arg.copy());
- int 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<p
- // kase = 2 if s(k) is negligible and k<p
- // kase = 3 if e[k-1] is negligible, k<p, and
- // s(k), ..., s(p) are not negligible (qr step).
- // kase = 4 if e(p-1) is negligible (convergence).
-
- for (k = p-2; 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<Real> &A)
- {
- int minm = std::min(m+1,n);
-
- A = Array2D<Real>(m, minm);
-
- for (int i=0; i<m; i++)
- for (int j=0; j<minm; j++)
- A[i][j] = U[i][j];
-
- }
-
- /* Return the right singular vectors */
-
- void getV (Array2D<Real> &A)
- {
- A = V;
- }
-
- /** Return the one-dimensional array of singular values */
-
- void getSingularValues (Array1D<Real> &x)
- {
- x = s;
- }
-
- /** Return the diagonal matrix of singular values
- @return S
- */
-
- void getS (Array2D<Real> &A) {
- A = Array2D<Real>(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 <cstdlib>
-#include <iostream>
-
-#ifdef TNT_BOUNDS_CHECK
-#include <assert.h>
-#endif
-
-
-#include "tnt_i_refvec.h"
-
-namespace TNT
-{
-
-template <class T>
-class Array1D
-{
-
- private:
-
- /* ... */
- i_refvec<T> 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<T> subarray(int i0, int i1);
-
-};
-
-
-
-
-template <class T>
-Array1D<T>::Array1D() : v_(), n_(0), data_(0) {}
-
-template <class T>
-Array1D<T>::Array1D(const Array1D<T> &A) : v_(A.v_), n_(A.n_),
- data_(A.data_)
-{
-#ifdef TNT_DEBUG
- std::cout << "Created Array1D(const Array1D<T> &A) \n";
-#endif
-
-}
-
-
-template <class T>
-Array1D<T>::Array1D(int n) : v_(n), n_(n), data_(v_.begin())
-{
-#ifdef TNT_DEBUG
- std::cout << "Created Array1D(int n) \n";
-#endif
-}
-
-template <class T>
-Array1D<T>::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 <class T>
-Array1D<T>::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 <class T>
-inline Array1D<T>::operator T*()
-{
- return &(v_[0]);
-}
-
-
-template <class T>
-inline Array1D<T>::operator const T*()
-{
- return &(v_[0]);
-}
-
-
-
-template <class T>
-inline T& Array1D<T>::operator[](int i)
-{
-#ifdef TNT_BOUNDS_CHECK
- assert(i>= 0);
- assert(i < n_);
-#endif
- return data_[i];
-}
-
-template <class T>
-inline const T& Array1D<T>::operator[](int i) const
-{
-#ifdef TNT_BOUNDS_CHECK
- assert(i>= 0);
- assert(i < n_);
-#endif
- return data_[i];
-}
-
-
-
-
-template <class T>
-Array1D<T> & Array1D<T>::operator=(const T &a)
-{
- set_(data_, data_+n_, a);
- return *this;
-}
-
-template <class T>
-Array1D<T> Array1D<T>::copy() const
-{
- Array1D A( n_);
- copy_(A.data_, data_, n_);
-
- return A;
-}
-
-
-template <class T>
-Array1D<T> & Array1D<T>::inject(const Array1D &A)
-{
- if (A.n_ == n_)
- copy_(data_, A.data_, n_);
-
- return *this;
-}
-
-
-
-
-
-template <class T>
-Array1D<T> & Array1D<T>::ref(const Array1D<T> &A)
-{
- if (this != &A)
- {
- v_ = A.v_; /* operator= handles the reference counting. */
- n_ = A.n_;
- data_ = A.data_;
-
- }
- return *this;
-}
-
-template <class T>
-Array1D<T> & Array1D<T>::operator=(const Array1D<T> &A)
-{
- return ref(A);
-}
-
-template <class T>
-inline int Array1D<T>::dim1() const { return n_; }
-
-template <class T>
-inline int Array1D<T>::dim() const { return n_; }
-
-template <class T>
-Array1D<T>::~Array1D() {}
-
-
-/* ............................ exented interface ......................*/
-
-template <class T>
-inline int Array1D<T>::ref_count() const
-{
- return v_.ref_count();
-}
-
-template <class T>
-inline Array1D<T> Array1D<T>::subarray(int i0, int i1)
-{
- if ((i0 > 0) && (i1 < n_) || (i0 <= i1))
- {
- Array1D<T> X(*this); /* create a new instance of this array. */
- X.n_ = i1-i0+1;
- X.data_ += i0;
-
- return X;
- }
- else
- {
- return Array1D<T>();
- }
-}
-
-
-/* private internal functions */
-
-
-template <class T>
-void Array1D<T>::set_(T* begin, T* end, const T& a)
-{
- for (T* p=begin; p<end; p++)
- *p = a;
-
-}
-
-template <class T>
-void Array1D<T>::copy_(T* p, const T* q, int len) const
-{
- T *end = p + len;
- while (p<end )
- *p++ = *q++;
-
-}
-
-
-} /* namespace TNT */
-
-#endif
-/* TNT_ARRAY1D_H */
-
diff --git a/lib/include/tnt/tnt_array1d_utils.h b/lib/include/tnt/tnt_array1d_utils.h
deleted file mode 100644
index 683e0e2..0000000
--- a/lib/include/tnt/tnt_array1d_utils.h
+++ /dev/null
@@ -1,230 +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_UTILS_H
-#define TNT_ARRAY1D_UTILS_H
-
-#include <cstdlib>
-#include <cassert>
-
-namespace TNT
-{
-
-
-template <class T>
-std::ostream& operator<<(std::ostream &s, const Array1D<T> &A)
-{
- int N=A.dim1();
-
-#ifdef TNT_DEBUG
- s << "addr: " << (void *) &A[0] << "\n";
-#endif
- s << N << "\n";
- for (int j=0; j<N; j++)
- {
- s << A[j] << "\n";
- }
- s << "\n";
-
- return s;
-}
-
-template <class T>
-std::istream& operator>>(std::istream &s, Array1D<T> &A)
-{
- int N;
- s >> N;
-
- Array1D<T> B(N);
- for (int i=0; i<N; i++)
- s >> B[i];
- A = B;
- return s;
-}
-
-
-
-template <class T>
-Array1D<T> operator+(const Array1D<T> &A, const Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() != n )
- return Array1D<T>();
-
- else
- {
- Array1D<T> C(n);
-
- for (int i=0; i<n; i++)
- {
- C[i] = A[i] + B[i];
- }
- return C;
- }
-}
-
-
-
-template <class T>
-Array1D<T> operator-(const Array1D<T> &A, const Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() != n )
- return Array1D<T>();
-
- else
- {
- Array1D<T> C(n);
-
- for (int i=0; i<n; i++)
- {
- C[i] = A[i] - B[i];
- }
- return C;
- }
-}
-
-
-template <class T>
-Array1D<T> operator*(const Array1D<T> &A, const Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() != n )
- return Array1D<T>();
-
- else
- {
- Array1D<T> C(n);
-
- for (int i=0; i<n; i++)
- {
- C[i] = A[i] * B[i];
- }
- return C;
- }
-}
-
-
-template <class T>
-Array1D<T> operator/(const Array1D<T> &A, const Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() != n )
- return Array1D<T>();
-
- else
- {
- Array1D<T> C(n);
-
- for (int i=0; i<n; i++)
- {
- C[i] = A[i] / B[i];
- }
- return C;
- }
-}
-
-
-
-
-
-
-
-
-
-template <class T>
-Array1D<T>& operator+=(Array1D<T> &A, const Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() == n)
- {
- for (int i=0; i<n; i++)
- {
- A[i] += B[i];
- }
- }
- return A;
-}
-
-
-
-
-template <class T>
-Array1D<T>& operator-=(Array1D<T> &A, const Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() == n)
- {
- for (int i=0; i<n; i++)
- {
- A[i] -= B[i];
- }
- }
- return A;
-}
-
-
-
-template <class T>
-Array1D<T>& operator*=(Array1D<T> &A, const Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() == n)
- {
- for (int i=0; i<n; i++)
- {
- A[i] *= B[i];
- }
- }
- return A;
-}
-
-
-
-
-template <class T>
-Array1D<T>& operator/=(Array1D<T> &A, const Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() == n)
- {
- for (int i=0; i<n; i++)
- {
- A[i] /= B[i];
- }
- }
- return A;
-}
-
-
-
-
-
-
-} // namespace TNT
-
-#endif
diff --git a/lib/include/tnt/tnt_array2d.h b/lib/include/tnt/tnt_array2d.h
deleted file mode 100644
index c791575..0000000
--- a/lib/include/tnt/tnt_array2d.h
+++ /dev/null
@@ -1,315 +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_H
-#define TNT_ARRAY2D_H
-
-#include <cstdlib>
-#include <iostream>
-#ifdef TNT_BOUNDS_CHECK
-#include <assert.h>
-#endif
-
-#include "tnt_array1d.h"
-
-namespace TNT
-{
-
-template <class T>
-class Array2D
-{
-
-
- private:
-
-
-
- Array1D<T> data_;
- Array1D<T*> 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 <class T>
-Array2D<T>::Array2D() : data_(), v_(), m_(0), n_(0) {}
-
-template <class T>
-Array2D<T>::Array2D(const Array2D<T> &A) : data_(A.data_), v_(A.v_),
- m_(A.m_), n_(A.n_) {}
-
-
-
-
-template <class T>
-Array2D<T>::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<m; i++)
- {
- v_[i] = p;
- p += n;
- }
- }
-}
-
-
-
-template <class T>
-Array2D<T>::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<m; i++)
- {
- v_[i] = p;
- p += n;
- }
- }
-}
-
-template <class T>
-Array2D<T>::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<m; i++)
- {
- v_[i] = p;
- p += n;
- }
- }
-}
-
-
-template <class T>
-inline T* Array2D<T>::operator[](int i)
-{
-#ifdef TNT_BOUNDS_CHECK
- assert(i >= 0);
- assert(i < m_);
-#endif
-
-return v_[i];
-
-}
-
-
-template <class T>
-inline const T* Array2D<T>::operator[](int i) const
-{
-#ifdef TNT_BOUNDS_CHECK
- assert(i >= 0);
- assert(i < m_);
-#endif
-
-return v_[i];
-
-}
-
-template <class T>
-Array2D<T> & Array2D<T>::operator=(const T &a)
-{
- /* non-optimzied, but will work with subarrays in future verions */
-
- for (int i=0; i<m_; i++)
- for (int j=0; j<n_; j++)
- v_[i][j] = a;
- return *this;
-}
-
-
-
-
-template <class T>
-Array2D<T> Array2D<T>::copy() const
-{
- Array2D A(m_, n_);
-
- for (int i=0; i<m_; i++)
- for (int j=0; j<n_; j++)
- A[i][j] = v_[i][j];
-
-
- return A;
-}
-
-
-template <class T>
-Array2D<T> & Array2D<T>::inject(const Array2D &A)
-{
- if (A.m_ == m_ && A.n_ == n_)
- {
- for (int i=0; i<m_; i++)
- for (int j=0; j<n_; j++)
- v_[i][j] = A[i][j];
- }
- return *this;
-}
-
-
-
-
-template <class T>
-Array2D<T> & Array2D<T>::ref(const Array2D<T> &A)
-{
- if (this != &A)
- {
- v_ = A.v_;
- data_ = A.data_;
- m_ = A.m_;
- n_ = A.n_;
-
- }
- return *this;
-}
-
-
-
-template <class T>
-Array2D<T> & Array2D<T>::operator=(const Array2D<T> &A)
-{
- return ref(A);
-}
-
-template <class T>
-inline int Array2D<T>::dim1() const { return m_; }
-
-template <class T>
-inline int Array2D<T>::dim2() const { return n_; }
-
-
-template <class T>
-Array2D<T>::~Array2D() {}
-
-
-
-
-template <class T>
-inline Array2D<T>::operator T**()
-{
- return &(v_[0]);
-}
-template <class T>
-inline Array2D<T>::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 <class T>
-Array2D<T> Array2D<T>::subarray(int i0, int i1, int j0, int j1)
-{
- Array2D<T> 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<T*>(m);
- T* p = &(data_[0]) + i0 * n_ + j0;
- for (int i=0; i<m; i++)
- {
- A.v_[i] = p + i*n_;
-
- }
- return A;
-}
-
-template <class T>
-inline int Array2D<T>::ref_count()
-{
- return ref_count_data();
-}
-
-
-
-template <class T>
-inline int Array2D<T>::ref_count_data()
-{
- return data_.ref_count();
-}
-
-template <class T>
-inline int Array2D<T>::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 <cstdlib>
-#include <cassert>
-
-namespace TNT
-{
-
-
-template <class T>
-std::ostream& operator<<(std::ostream &s, const Array2D<T> &A)
-{
- int M=A.dim1();
- int N=A.dim2();
-
- s << M << " " << N << "\n";
-
- for (int i=0; i<M; i++)
- {
- for (int j=0; j<N; j++)
- {
- s << A[i][j] << " ";
- }
- s << "\n";
- }
-
-
- return s;
-}
-
-template <class T>
-std::istream& operator>>(std::istream &s, Array2D<T> &A)
-{
-
- int M, N;
-
- s >> M >> N;
-
- Array2D<T> B(M,N);
-
- for (int i=0; i<M; i++)
- for (int j=0; j<N; j++)
- {
- s >> B[i][j];
- }
-
- A = B;
- return s;
-}
-
-
-template <class T>
-Array2D<T> operator+(const Array2D<T> &A, const Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() != m || B.dim2() != n )
- return Array2D<T>();
-
- else
- {
- Array2D<T> C(m,n);
-
- for (int i=0; i<m; i++)
- {
- for (int j=0; j<n; j++)
- C[i][j] = A[i][j] + B[i][j];
- }
- return C;
- }
-}
-
-template <class T>
-Array2D<T> operator-(const Array2D<T> &A, const Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() != m || B.dim2() != n )
- return Array2D<T>();
-
- else
- {
- Array2D<T> C(m,n);
-
- for (int i=0; i<m; i++)
- {
- for (int j=0; j<n; j++)
- C[i][j] = A[i][j] - B[i][j];
- }
- return C;
- }
-}
-
-
-template <class T>
-Array2D<T> operator*(const Array2D<T> &A, const Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() != m || B.dim2() != n )
- return Array2D<T>();
-
- else
- {
- Array2D<T> C(m,n);
-
- for (int i=0; i<m; i++)
- {
- for (int j=0; j<n; j++)
- C[i][j] = A[i][j] * B[i][j];
- }
- return C;
- }
-}
-
-
-
-
-template <class T>
-Array2D<T> operator/(const Array2D<T> &A, const Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() != m || B.dim2() != n )
- return Array2D<T>();
-
- else
- {
- Array2D<T> C(m,n);
-
- for (int i=0; i<m; i++)
- {
- for (int j=0; j<n; j++)
- C[i][j] = A[i][j] / B[i][j];
- }
- return C;
- }
-}
-
-
-
-
-
-template <class T>
-Array2D<T>& operator+=(Array2D<T> &A, const Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() == m || B.dim2() == n )
- {
- for (int i=0; i<m; i++)
- {
- for (int j=0; j<n; j++)
- A[i][j] += B[i][j];
- }
- }
- return A;
-}
-
-
-
-template <class T>
-Array2D<T>& operator-=(Array2D<T> &A, const Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() == m || B.dim2() == n )
- {
- for (int i=0; i<m; i++)
- {
- for (int j=0; j<n; j++)
- A[i][j] -= B[i][j];
- }
- }
- return A;
-}
-
-
-
-template <class T>
-Array2D<T>& operator*=(Array2D<T> &A, const Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() == m || B.dim2() == n )
- {
- for (int i=0; i<m; i++)
- {
- for (int j=0; j<n; j++)
- A[i][j] *= B[i][j];
- }
- }
- return A;
-}
-
-
-
-
-
-template <class T>
-Array2D<T>& operator/=(Array2D<T> &A, const Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() == m || B.dim2() == n )
- {
- for (int i=0; i<m; i++)
- {
- for (int j=0; j<n; j++)
- A[i][j] /= B[i][j];
- }
- }
- return A;
-}
-
-/**
- Matrix Multiply: compute C = A*B, where C[i][j]
- is the dot-product of row i of A and column j of B.
-
-
- @param A an (m x n) array
- @param B an (n x k) array
- @return the (m x k) array A*B, or a null array (0x0)
- if the matrices are non-conformant (i.e. the number
- of columns of A are different than the number of rows of B.)
-
-
-*/
-template <class T>
-Array2D<T> matmult(const Array2D<T> &A, const Array2D<T> &B)
-{
- if (A.dim2() != B.dim1())
- return Array2D<T>();
-
- int M = A.dim1();
- int N = A.dim2();
- int K = B.dim2();
-
- Array2D<T> C(M,K);
-
- for (int i=0; i<M; i++)
- for (int j=0; j<K; j++)
- {
- T sum = 0;
-
- for (int k=0; k<N; k++)
- sum += A[i][k] * B [k][j];
-
- C[i][j] = sum;
- }
-
- return C;
-
-}
-
-} // namespace TNT
-
-#endif
diff --git a/lib/include/tnt/tnt_array3d.h b/lib/include/tnt/tnt_array3d.h
deleted file mode 100644
index c210d2e..0000000
--- a/lib/include/tnt/tnt_array3d.h
+++ /dev/null
@@ -1,296 +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_ARRAY3D_H
-#define TNT_ARRAY3D_H
-
-#include <cstdlib>
-#include <iostream>
-#ifdef TNT_BOUNDS_CHECK
-#include <assert.h>
-#endif
-
-#include "tnt_array1d.h"
-#include "tnt_array2d.h"
-
-namespace TNT
-{
-
-template <class T>
-class Array3D
-{
-
-
- private:
- Array1D<T> data_;
- Array2D<T*> 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 <class T>
-Array3D<T>::Array3D() : data_(), v_(), m_(0), n_(0) {}
-
-template <class T>
-Array3D<T>::Array3D(const Array3D<T> &A) : data_(A.data_),
- v_(A.v_), m_(A.m_), n_(A.n_), g_(A.g_)
-{
-}
-
-
-
-template <class T>
-Array3D<T>::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<m_; i++)
- {
- T* ping = p+ i*ng;
- for (int j=0; j<n; j++)
- v_[i][j] = ping + j*g_;
- }
- }
-}
-
-
-
-template <class T>
-Array3D<T>::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<m_; i++)
- {
- T* ping = p+ i*ng;
- for (int j=0; j<n; j++)
- v_[i][j] = ping + j*g_;
- }
- }
-}
-
-
-
-template <class T>
-Array3D<T>::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<m_; i++)
- {
- T* ping = p+ i*ng;
- for (int j=0; j<n; j++)
- v_[i][j] = ping + j*g_;
- }
- }
-}
-
-
-
-template <class T>
-inline T** Array3D<T>::operator[](int i)
-{
-#ifdef TNT_BOUNDS_CHECK
- assert(i >= 0);
- assert(i < m_);
-#endif
-
-return v_[i];
-
-}
-
-template <class T>
-inline const T* const * Array3D<T>::operator[](int i) const
-{ return v_[i]; }
-
-template <class T>
-Array3D<T> & Array3D<T>::operator=(const T &a)
-{
- for (int i=0; i<m_; i++)
- for (int j=0; j<n_; j++)
- for (int k=0; k<g_; k++)
- v_[i][j][k] = a;
-
- return *this;
-}
-
-template <class T>
-Array3D<T> Array3D<T>::copy() const
-{
- Array3D A(m_, n_, g_);
- for (int i=0; i<m_; i++)
- for (int j=0; j<n_; j++)
- for (int k=0; k<g_; k++)
- A.v_[i][j][k] = v_[i][j][k];
-
- return A;
-}
-
-
-template <class T>
-Array3D<T> & Array3D<T>::inject(const Array3D &A)
-{
- if (A.m_ == m_ && A.n_ == n_ && A.g_ == g_)
-
- for (int i=0; i<m_; i++)
- for (int j=0; j<n_; j++)
- for (int k=0; k<g_; k++)
- v_[i][j][k] = A.v_[i][j][k];
-
- return *this;
-}
-
-
-
-template <class T>
-Array3D<T> & Array3D<T>::ref(const Array3D<T> &A)
-{
- if (this != &A)
- {
- m_ = A.m_;
- n_ = A.n_;
- g_ = A.g_;
- v_ = A.v_;
- data_ = A.data_;
- }
- return *this;
-}
-
-template <class T>
-Array3D<T> & Array3D<T>::operator=(const Array3D<T> &A)
-{
- return ref(A);
-}
-
-
-template <class T>
-inline int Array3D<T>::dim1() const { return m_; }
-
-template <class T>
-inline int Array3D<T>::dim2() const { return n_; }
-
-template <class T>
-inline int Array3D<T>::dim3() const { return g_; }
-
-
-
-template <class T>
-Array3D<T>::~Array3D() {}
-
-template <class T>
-inline Array3D<T>::operator T***()
-{
- return v_;
-}
-
-
-template <class T>
-inline Array3D<T>::operator const T***()
-{
- return v_;
-}
-
-/* extended interface */
-template <class T>
-Array3D<T> Array3D<T>::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<T>(); /* null array */
-
-
- Array3D<T> A;
- A.data_ = data_;
- A.m_ = i1-i0+1;
- A.n_ = j1-j0+1;
- A.g_ = k1-k0+1;
- A.v_ = Array2D<T*>(A.m_,A.n_);
- T* p = &(data_[0]) + i0*n_*g_ + j0*g_ + k0;
-
- for (int i=0; i<A.m_; i++)
- {
- T* ping = p + i*n_*g_;
- for (int j=0; j<A.n_; j++)
- A.v_[i][j] = ping + j*g_ ;
- }
-
- return A;
-}
-
-
-
-} /* namespace TNT */
-
-#endif
-/* TNT_ARRAY3D_H */
-
diff --git a/lib/include/tnt/tnt_array3d_utils.h b/lib/include/tnt/tnt_array3d_utils.h
deleted file mode 100644
index 5acdc1d..0000000
--- a/lib/include/tnt/tnt_array3d_utils.h
+++ /dev/null
@@ -1,236 +0,0 @@
-
-
-#ifndef TNT_ARRAY3D_UTILS_H
-#define TNT_ARRAY3D_UTILS_H
-
-#include <cstdlib>
-#include <cassert>
-
-namespace TNT
-{
-
-
-template <class T>
-std::ostream& operator<<(std::ostream &s, const Array3D<T> &A)
-{
- int M=A.dim1();
- int N=A.dim2();
- int K=A.dim3();
-
- s << M << " " << N << " " << K << "\n";
-
- for (int i=0; i<M; i++)
- {
- for (int j=0; j<N; j++)
- {
- for (int k=0; k<K; k++)
- s << A[i][j][k] << " ";
- s << "\n";
- }
- s << "\n";
- }
-
-
- return s;
-}
-
-template <class T>
-std::istream& operator>>(std::istream &s, Array3D<T> &A)
-{
-
- int M, N, K;
-
- s >> M >> N >> K;
-
- Array3D<T> B(M,N,K);
-
- for (int i=0; i<M; i++)
- for (int j=0; j<N; j++)
- for (int k=0; k<K; k++)
- s >> B[i][j][k];
-
- A = B;
- return s;
-}
-
-
-
-template <class T>
-Array3D<T> operator+(const Array3D<T> &A, const Array3D<T> &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<T>();
-
- else
- {
- Array3D<T> C(m,n,p);
-
- for (int i=0; i<m; i++)
- for (int j=0; j<n; j++)
- for (int k=0; k<p; k++)
- C[i][j][k] = A[i][j][k] + B[i][j][k];
-
- return C;
- }
-}
-
-
-template <class T>
-Array3D<T> operator-(const Array3D<T> &A, const Array3D<T> &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<T>();
-
- else
- {
- Array3D<T> C(m,n,p);
-
- for (int i=0; i<m; i++)
- for (int j=0; j<n; j++)
- for (int k=0; k<p; k++)
- C[i][j][k] = A[i][j][k] - B[i][j][k];
-
- return C;
- }
-}
-
-
-
-
-template <class T>
-Array3D<T> operator*(const Array3D<T> &A, const Array3D<T> &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<T>();
-
- else
- {
- Array3D<T> C(m,n,p);
-
- for (int i=0; i<m; i++)
- for (int j=0; j<n; j++)
- for (int k=0; k<p; k++)
- C[i][j][k] = A[i][j][k] * B[i][j][k];
-
- return C;
- }
-}
-
-
-template <class T>
-Array3D<T> operator/(const Array3D<T> &A, const Array3D<T> &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<T>();
-
- else
- {
- Array3D<T> C(m,n,p);
-
- for (int i=0; i<m; i++)
- for (int j=0; j<n; j++)
- for (int k=0; k<p; k++)
- C[i][j][k] = A[i][j][k] / B[i][j][k];
-
- return C;
- }
-}
-
-
-
-template <class T>
-Array3D<T>& operator+=(Array3D<T> &A, const Array3D<T> &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<m; i++)
- for (int j=0; j<n; j++)
- for (int k=0; k<p; k++)
- A[i][j][k] += B[i][j][k];
- }
-
- return A;
-}
-
-template <class T>
-Array3D<T>& operator-=(Array3D<T> &A, const Array3D<T> &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<m; i++)
- for (int j=0; j<n; j++)
- for (int k=0; k<p; k++)
- A[i][j][k] -= B[i][j][k];
- }
-
- return A;
-}
-
-template <class T>
-Array3D<T>& operator*=(Array3D<T> &A, const Array3D<T> &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<m; i++)
- for (int j=0; j<n; j++)
- for (int k=0; k<p; k++)
- A[i][j][k] *= B[i][j][k];
- }
-
- return A;
-}
-
-
-template <class T>
-Array3D<T>& operator/=(Array3D<T> &A, const Array3D<T> &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<m; i++)
- for (int j=0; j<n; j++)
- for (int k=0; k<p; k++)
- A[i][j][k] /= B[i][j][k];
- }
-
- return A;
-}
-
-
-
-
-
-} // namespace TNT
-
-#endif
diff --git a/lib/include/tnt/tnt_cmat.h b/lib/include/tnt/tnt_cmat.h
deleted file mode 100644
index 5ff4c48..0000000
--- a/lib/include/tnt/tnt_cmat.h
+++ /dev/null
@@ -1,580 +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.
-*
-*/
-
-
-// C compatible matrix: row-oriented, 0-based [i][j] and 1-based (i,j) indexing
-//
-
-#ifndef TNT_CMAT_H
-#define TNT_CMAT_H
-
-#include "tnt_subscript.h"
-#include "tnt_vec.h"
-#include <cstdlib>
-#include <cassert>
-#include <iostream>
-#include <sstream>
-
-namespace TNT
-{
-
-
-template <class T>
-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<M; i++)
- {
- row_[i] = p;
- rowm1_[i] = p-1;
- p += N ;
-
- }
-
- rowm1_ -- ; // compensate for 1-based offset
- }
-
- void copy(const T* v)
- {
- Subscript N = m_ * n_;
- Subscript i;
-
-#ifdef TNT_UNROLL_LOOPS
- Subscript Nmod4 = N & 3;
- Subscript N4 = N - Nmod4;
-
- for (i=0; i<N4; i+=4)
- {
- v_[i] = v[i];
- v_[i+1] = v[i+1];
- v_[i+2] = v[i+2];
- v_[i+3] = v[i+3];
- }
-
- for (i=N4; i< N; i++)
- v_[i] = v[i];
-#else
-
- for (i=0; i< N; i++)
- v_[i] = v[i];
-#endif
- }
-
- void set(const T& val)
- {
- Subscript N = m_ * n_;
- Subscript i;
-
-#ifdef TNT_UNROLL_LOOPS
- Subscript Nmod4 = N & 3;
- Subscript N4 = N - Nmod4;
-
- for (i=0; i<N4; i+=4)
- {
- v_[i] = val;
- v_[i+1] = val;
- v_[i+2] = val;
- v_[i+3] = val;
- }
-
- for (i=N4; i< N; i++)
- v_[i] = val;
-#else
-
- for (i=0; i< N; i++)
- v_[i] = val;
-
-#endif
- }
-
-
-
- void destroy()
- {
- /* do nothing, if no memory has been previously allocated */
- if (v_ == NULL) return ;
-
- /* if we are here, then matrix was previously allocated */
- if (v_ != NULL) delete [] (v_);
- if (row_ != NULL) delete [] (row_);
-
- /* return rowm1_ back to original value */
- rowm1_ ++;
- if (rowm1_ != NULL ) delete [] (rowm1_);
- }
-
-
- public:
-
- operator T**(){ return row_; }
- operator T**() const { return row_; }
-
-
- Subscript size() const { return mn_; }
-
- // constructors
-
- Matrix() : m_(0), n_(0), mn_(0), v_(0), row_(0), vm1_(0), rowm1_(0) {};
-
- Matrix(const Matrix<T> &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<M; i++)
- for (j=0; j<N; j++)
- ins >> row_[i][j];
- }
-
- // destructor
- //
- ~Matrix()
- {
- destroy();
- }
-
-
- // reallocating
- //
- Matrix<T>& newsize(Subscript M, Subscript N)
- {
- if (num_rows() == M && num_cols() == N)
- return *this;
-
- destroy();
- initialize(M,N);
-
- return *this;
- }
-
-
-
-
- // assignments
- //
- Matrix<T>& operator=(const Matrix<T> &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<T>& 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 <class T>
-std::ostream& operator<<(std::ostream &s, const Matrix<T> &A)
-{
- Subscript M=A.num_rows();
- Subscript N=A.num_cols();
-
- s << M << " " << N << "\n";
-
- for (Subscript i=0; i<M; i++)
- {
- for (Subscript j=0; j<N; j++)
- {
- s << A[i][j] << " ";
- }
- s << "\n";
- }
-
-
- return s;
-}
-
-template <class T>
-std::istream& operator>>(std::istream &s, Matrix<T> &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<M; i++)
- for (Subscript j=0; j<N; j++)
- {
- s >> A[i][j];
- }
-
-
- return s;
-}
-
-// *******************[ basic matrix algorithms ]***************************
-
-
-template <class T>
-Matrix<T> operator+(const Matrix<T> &A,
- const Matrix<T> &B)
-{
- Subscript M = A.num_rows();
- Subscript N = A.num_cols();
-
- assert(M==B.num_rows());
- assert(N==B.num_cols());
-
- Matrix<T> tmp(M,N);
- Subscript i,j;
-
- for (i=0; i<M; i++)
- for (j=0; j<N; j++)
- tmp[i][j] = A[i][j] + B[i][j];
-
- return tmp;
-}
-
-template <class T>
-Matrix<T> operator-(const Matrix<T> &A,
- const Matrix<T> &B)
-{
- Subscript M = A.num_rows();
- Subscript N = A.num_cols();
-
- assert(M==B.num_rows());
- assert(N==B.num_cols());
-
- Matrix<T> tmp(M,N);
- Subscript i,j;
-
- for (i=0; i<M; i++)
- for (j=0; j<N; j++)
- tmp[i][j] = A[i][j] - B[i][j];
-
- return tmp;
-}
-
-template <class T>
-Matrix<T> mult_element(const Matrix<T> &A,
- const Matrix<T> &B)
-{
- Subscript M = A.num_rows();
- Subscript N = A.num_cols();
-
- assert(M==B.num_rows());
- assert(N==B.num_cols());
-
- Matrix<T> tmp(M,N);
- Subscript i,j;
-
- for (i=0; i<M; i++)
- for (j=0; j<N; j++)
- tmp[i][j] = A[i][j] * B[i][j];
-
- return tmp;
-}
-
-
-template <class T>
-Matrix<T> transpose(const Matrix<T> &A)
-{
- Subscript M = A.num_rows();
- Subscript N = A.num_cols();
-
- Matrix<T> S(N,M);
- Subscript i, j;
-
- for (i=0; i<M; i++)
- for (j=0; j<N; j++)
- S[j][i] = A[i][j];
-
- return S;
-}
-
-
-
-template <class T>
-inline Matrix<T> matmult(const Matrix<T> &A,
- const Matrix<T> &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<T> tmp(M,K);
- T sum;
-
- for (Subscript i=0; i<M; i++)
- for (Subscript k=0; k<K; k++)
- {
- sum = 0;
- for (Subscript j=0; j<N; j++)
- sum = sum + A[i][j] * B[j][k];
-
- tmp[i][k] = sum;
- }
-
- return tmp;
-}
-
-template <class T>
-inline Matrix<T> operator*(const Matrix<T> &A,
- const Matrix<T> &B)
-{
- return matmult(A,B);
-}
-
-template <class T>
-inline int matmult(Matrix<T>& C, const Matrix<T> &A,
- const Matrix<T> &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<M; i++)
- for (Subscript k=0; k<K; k++)
- {
- row_i = &(A[i][0]);
- col_k = &(B[0][k]);
- sum = 0;
- for (Subscript j=0; j<N; j++)
- {
- sum += *row_i * *col_k;
- row_i++;
- col_k += K;
- }
- C[i][k] = sum;
- }
-
- return 0;
-}
-
-
-template <class T>
-Vector<T> matmult(const Matrix<T> &A, const Vector<T> &x)
-{
-
-#ifdef TNT_BOUNDS_CHECK
- assert(A.num_cols() == x.dim());
-#endif
-
- Subscript M = A.num_rows();
- Subscript N = A.num_cols();
-
- Vector<T> tmp(M);
- T sum;
-
- for (Subscript i=0; i<M; i++)
- {
- sum = 0;
- const T* rowi = A[i];
- for (Subscript j=0; j<N; j++)
- sum = sum + rowi[j] * x[j];
-
- tmp[i] = sum;
- }
-
- return tmp;
-}
-
-template <class T>
-inline Vector<T> operator*(const Matrix<T> &A, const Vector<T> &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 <cstdlib>
-#include <iostream>
-
-#ifdef TNT_BOUNDS_CHECK
-#include <assert.h>
-#endif
-
-
-#include "tnt_i_refvec.h"
-
-namespace TNT
-{
-
-template <class T>
-class Fortran_Array1D
-{
-
- private:
-
- i_refvec<T> 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<T> subarray(int i0, int i1);
-
-};
-
-
-
-
-template <class T>
-Fortran_Array1D<T>::Fortran_Array1D() : v_(), n_(0), data_(0) {}
-
-template <class T>
-Fortran_Array1D<T>::Fortran_Array1D(const Fortran_Array1D<T> &A) : v_(A.v_), n_(A.n_),
- data_(A.data_)
-{
-#ifdef TNT_DEBUG
- std::cout << "Created Fortran_Array1D(const Fortran_Array1D<T> &A) \n";
-#endif
-
-}
-
-
-template <class T>
-Fortran_Array1D<T>::Fortran_Array1D(int n) : v_(n), n_(n), data_(v_.begin())
-{
-#ifdef TNT_DEBUG
- std::cout << "Created Fortran_Array1D(int n) \n";
-#endif
-}
-
-template <class T>
-Fortran_Array1D<T>::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 <class T>
-Fortran_Array1D<T>::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 <class T>
-inline T& Fortran_Array1D<T>::operator()(int i)
-{
-#ifdef TNT_BOUNDS_CHECK
- assert(i>= 1);
- assert(i <= n_);
-#endif
- return data_[i-1];
-}
-
-template <class T>
-inline const T& Fortran_Array1D<T>::operator()(int i) const
-{
-#ifdef TNT_BOUNDS_CHECK
- assert(i>= 1);
- assert(i <= n_);
-#endif
- return data_[i-1];
-}
-
-
-
-
-template <class T>
-Fortran_Array1D<T> & Fortran_Array1D<T>::operator=(const T &a)
-{
- set_(data_, data_+n_, a);
- return *this;
-}
-
-template <class T>
-Fortran_Array1D<T> Fortran_Array1D<T>::copy() const
-{
- Fortran_Array1D A( n_);
- copy_(A.data_, data_, n_);
-
- return A;
-}
-
-
-template <class T>
-Fortran_Array1D<T> & Fortran_Array1D<T>::inject(const Fortran_Array1D &A)
-{
- if (A.n_ == n_)
- copy_(data_, A.data_, n_);
-
- return *this;
-}
-
-
-
-
-
-template <class T>
-Fortran_Array1D<T> & Fortran_Array1D<T>::ref(const Fortran_Array1D<T> &A)
-{
- if (this != &A)
- {
- v_ = A.v_; /* operator= handles the reference counting. */
- n_ = A.n_;
- data_ = A.data_;
-
- }
- return *this;
-}
-
-template <class T>
-Fortran_Array1D<T> & Fortran_Array1D<T>::operator=(const Fortran_Array1D<T> &A)
-{
- return ref(A);
-}
-
-template <class T>
-inline int Fortran_Array1D<T>::dim1() const { return n_; }
-
-template <class T>
-inline int Fortran_Array1D<T>::dim() const { return n_; }
-
-template <class T>
-Fortran_Array1D<T>::~Fortran_Array1D() {}
-
-
-/* ............................ exented interface ......................*/
-
-template <class T>
-inline int Fortran_Array1D<T>::ref_count() const
-{
- return v_.ref_count();
-}
-
-template <class T>
-inline Fortran_Array1D<T> Fortran_Array1D<T>::subarray(int i0, int i1)
-{
-#ifdef TNT_DEBUG
- std::cout << "entered subarray. \n";
-#endif
- if ((i0 > 0) && (i1 < n_) || (i0 <= i1))
- {
- Fortran_Array1D<T> 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<T>();
- }
-}
-
-
-/* private internal functions */
-
-
-template <class T>
-void Fortran_Array1D<T>::set_(T* begin, T* end, const T& a)
-{
- for (T* p=begin; p<end; p++)
- *p = a;
-
-}
-
-template <class T>
-void Fortran_Array1D<T>::copy_(T* p, const T* q, int len) const
-{
- T *end = p + len;
- while (p<end )
- *p++ = *q++;
-
-}
-
-
-} /* namespace TNT */
-
-#endif
-/* TNT_FORTRAN_ARRAY1D_H */
-
diff --git a/lib/include/tnt/tnt_fortran_array1d_utils.h b/lib/include/tnt/tnt_fortran_array1d_utils.h
deleted file mode 100644
index b037b17..0000000
--- a/lib/include/tnt/tnt_fortran_array1d_utils.h
+++ /dev/null
@@ -1,242 +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_UTILS_H
-#define TNT_FORTRAN_ARRAY1D_UTILS_H
-
-#include <iostream>
-
-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 <class T>
-std::ostream& operator<<(std::ostream &s, const Fortran_Array1D<T> &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
-
- <p>
- Note: the array being read into references new memory
- storage. If the intent is to fill an existing conformant
- array, use <code> cin >> B; A.inject(B) ); </code>
- instead or read the elements in one-a-time by hand.
-
- @param s the charater to read from (typically <code>std::in</code>)
- @param A the array to read into.
-*/
-template <class T>
-std::istream& operator>>(std::istream &s, Fortran_Array1D<T> &A)
-{
- int N;
- s >> N;
-
- Fortran_Array1D<T> B(N);
- for (int i=1; i<=N; i++)
- s >> B(i);
- A = B;
- return s;
-}
-
-
-template <class T>
-Fortran_Array1D<T> operator+(const Fortran_Array1D<T> &A, const Fortran_Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() != n )
- return Fortran_Array1D<T>();
-
- else
- {
- Fortran_Array1D<T> C(n);
-
- for (int i=1; i<=n; i++)
- {
- C(i) = A(i) + B(i);
- }
- return C;
- }
-}
-
-
-
-template <class T>
-Fortran_Array1D<T> operator-(const Fortran_Array1D<T> &A, const Fortran_Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() != n )
- return Fortran_Array1D<T>();
-
- else
- {
- Fortran_Array1D<T> C(n);
-
- for (int i=1; i<=n; i++)
- {
- C(i) = A(i) - B(i);
- }
- return C;
- }
-}
-
-
-template <class T>
-Fortran_Array1D<T> operator*(const Fortran_Array1D<T> &A, const Fortran_Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() != n )
- return Fortran_Array1D<T>();
-
- else
- {
- Fortran_Array1D<T> C(n);
-
- for (int i=1; i<=n; i++)
- {
- C(i) = A(i) * B(i);
- }
- return C;
- }
-}
-
-
-template <class T>
-Fortran_Array1D<T> operator/(const Fortran_Array1D<T> &A, const Fortran_Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() != n )
- return Fortran_Array1D<T>();
-
- else
- {
- Fortran_Array1D<T> C(n);
-
- for (int i=1; i<=n; i++)
- {
- C(i) = A(i) / B(i);
- }
- return C;
- }
-}
-
-
-
-
-
-
-
-
-
-template <class T>
-Fortran_Array1D<T>& operator+=(Fortran_Array1D<T> &A, const Fortran_Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() == n)
- {
- for (int i=1; i<=n; i++)
- {
- A(i) += B(i);
- }
- }
- return A;
-}
-
-
-
-
-template <class T>
-Fortran_Array1D<T>& operator-=(Fortran_Array1D<T> &A, const Fortran_Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() == n)
- {
- for (int i=1; i<=n; i++)
- {
- A(i) -= B(i);
- }
- }
- return A;
-}
-
-
-
-template <class T>
-Fortran_Array1D<T>& operator*=(Fortran_Array1D<T> &A, const Fortran_Array1D<T> &B)
-{
- int n = A.dim1();
-
- if (B.dim1() == n)
- {
- for (int i=1; i<=n; i++)
- {
- A(i) *= B(i);
- }
- }
- return A;
-}
-
-
-
-
-template <class T>
-Fortran_Array1D<T>& operator/=(Fortran_Array1D<T> &A, const Fortran_Array1D<T> &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 <cstdlib>
-#include <iostream>
-
-#ifdef TNT_BOUNDS_CHECK
-#include <assert.h>
-#endif
-
-#include "tnt_i_refvec.h"
-
-namespace TNT
-{
-
-template <class T>
-class Fortran_Array2D
-{
-
-
- private:
- i_refvec<T> 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 <class T>
-Fortran_Array2D<T>::Fortran_Array2D() : v_(), m_(0), n_(0), data_(0) {}
-
-
-template <class T>
-Fortran_Array2D<T>::Fortran_Array2D(const Fortran_Array2D<T> &A) : v_(A.v_),
- m_(A.m_), n_(A.n_), data_(A.data_) {}
-
-
-
-template <class T>
-Fortran_Array2D<T>::Fortran_Array2D(int m, int n) : v_(m*n), m_(m), n_(n),
- data_(v_.begin()) {}
-
-template <class T>
-Fortran_Array2D<T>::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 <class T>
-Fortran_Array2D<T>::Fortran_Array2D(int m, int n, T *a) : v_(a),
- m_(m), n_(n), data_(v_.begin()) {}
-
-
-
-
-template <class T>
-inline T& Fortran_Array2D<T>::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 <class T>
-inline const T& Fortran_Array2D<T>::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 <class T>
-Fortran_Array2D<T> & Fortran_Array2D<T>::operator=(const T &a)
-{
- set_(data_, data_+m_*n_, a);
- return *this;
-}
-
-template <class T>
-Fortran_Array2D<T> Fortran_Array2D<T>::copy() const
-{
-
- Fortran_Array2D B(m_,n_);
-
- B.inject(*this);
- return B;
-}
-
-
-template <class T>
-Fortran_Array2D<T> & Fortran_Array2D<T>::inject(const Fortran_Array2D &A)
-{
- if (m_ == A.m_ && n_ == A.n_)
- copy_(data_, A.data_, m_*n_);
-
- return *this;
-}
-
-
-
-template <class T>
-Fortran_Array2D<T> & Fortran_Array2D<T>::ref(const Fortran_Array2D<T> &A)
-{
- if (this != &A)
- {
- v_ = A.v_;
- m_ = A.m_;
- n_ = A.n_;
- data_ = A.data_;
- }
- return *this;
-}
-
-template <class T>
-Fortran_Array2D<T> & Fortran_Array2D<T>::operator=(const Fortran_Array2D<T> &A)
-{
- return ref(A);
-}
-
-template <class T>
-inline int Fortran_Array2D<T>::dim1() const { return m_; }
-
-template <class T>
-inline int Fortran_Array2D<T>::dim2() const { return n_; }
-
-
-template <class T>
-Fortran_Array2D<T>::~Fortran_Array2D()
-{
-}
-
-template <class T>
-inline int Fortran_Array2D<T>::ref_count() const { return v_.ref_count(); }
-
-
-
-
-template <class T>
-void Fortran_Array2D<T>::set_(T* begin, T* end, const T& a)
-{
- for (T* p=begin; p<end; p++)
- *p = a;
-
-}
-
-template <class T>
-void Fortran_Array2D<T>::copy_(T* p, const T* q, int len)
-{
- T *end = p + len;
- while (p<end )
- *p++ = *q++;
-
-}
-
-
-} /* namespace TNT */
-
-#endif
-/* TNT_FORTRAN_ARRAY2D_H */
-
diff --git a/lib/include/tnt/tnt_fortran_array2d_utils.h b/lib/include/tnt/tnt_fortran_array2d_utils.h
deleted file mode 100644
index bb68673..0000000
--- a/lib/include/tnt/tnt_fortran_array2d_utils.h
+++ /dev/null
@@ -1,236 +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_ARRAY2D_UTILS_H
-#define TNT_FORTRAN_ARRAY2D_UTILS_H
-
-#include <iostream>
-
-namespace TNT
-{
-
-
-template <class T>
-std::ostream& operator<<(std::ostream &s, const Fortran_Array2D<T> &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 <class T>
-std::istream& operator>>(std::istream &s, Fortran_Array2D<T> &A)
-{
-
- int M, N;
-
- s >> M >> N;
-
- Fortran_Array2D<T> 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 <class T>
-Fortran_Array2D<T> operator+(const Fortran_Array2D<T> &A, const Fortran_Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() != m || B.dim2() != n )
- return Fortran_Array2D<T>();
-
- else
- {
- Fortran_Array2D<T> 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 <class T>
-Fortran_Array2D<T> operator-(const Fortran_Array2D<T> &A, const Fortran_Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() != m || B.dim2() != n )
- return Fortran_Array2D<T>();
-
- else
- {
- Fortran_Array2D<T> 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 <class T>
-Fortran_Array2D<T> operator*(const Fortran_Array2D<T> &A, const Fortran_Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() != m || B.dim2() != n )
- return Fortran_Array2D<T>();
-
- else
- {
- Fortran_Array2D<T> 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 <class T>
-Fortran_Array2D<T> operator/(const Fortran_Array2D<T> &A, const Fortran_Array2D<T> &B)
-{
- int m = A.dim1();
- int n = A.dim2();
-
- if (B.dim1() != m || B.dim2() != n )
- return Fortran_Array2D<T>();
-
- else
- {
- Fortran_Array2D<T> 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 <class T>
-Fortran_Array2D<T>& operator+=(Fortran_Array2D<T> &A, const Fortran_Array2D<T> &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 <class T>
-Fortran_Array2D<T>& operator-=(Fortran_Array2D<T> &A, const Fortran_Array2D<T> &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 <class T>
-Fortran_Array2D<T>& operator*=(Fortran_Array2D<T> &A, const Fortran_Array2D<T> &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 <class T>
-Fortran_Array2D<T>& operator/=(Fortran_Array2D<T> &A, const Fortran_Array2D<T> &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 <cstdlib>
-#include <iostream>
-#ifdef TNT_BOUNDS_CHECK
-#include <assert.h>
-#endif
-#include "tnt_i_refvec.h"
-
-namespace TNT
-{
-
-template <class T>
-class Fortran_Array3D
-{
-
-
- private:
-
-
- i_refvec<T> 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 <class T>
-Fortran_Array3D<T>::Fortran_Array3D() : v_(), m_(0), n_(0), k_(0), data_(0) {}
-
-
-template <class T>
-Fortran_Array3D<T>::Fortran_Array3D(const Fortran_Array3D<T> &A) :
- v_(A.v_), m_(A.m_), n_(A.n_), k_(A.k_), data_(A.data_) {}
-
-
-
-template <class T>
-Fortran_Array3D<T>::Fortran_Array3D(int m, int n, int k) :
- v_(m*n*k), m_(m), n_(n), k_(k), data_(v_.begin()) {}
-
-
-
-template <class T>
-Fortran_Array3D<T>::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 <class T>
-Fortran_Array3D<T>::Fortran_Array3D(int m, int n, int k, T *a) :
- v_(a), m_(m), n_(n), k_(k), data_(v_.begin()) {}
-
-
-
-
-template <class T>
-inline T& Fortran_Array3D<T>::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 <class T>
-inline const T& Fortran_Array3D<T>::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 <class T>
-Fortran_Array3D<T> & Fortran_Array3D<T>::operator=(const T &a)
-{
-
- T *end = data_ + m_*n_*k_;
-
- for (T *p=data_; p != end; *p++ = a);
-
- return *this;
-}
-
-template <class T>
-Fortran_Array3D<T> Fortran_Array3D<T>::copy() const
-{
-
- Fortran_Array3D B(m_, n_, k_);
- B.inject(*this);
- return B;
-
-}
-
-
-template <class T>
-Fortran_Array3D<T> & Fortran_Array3D<T>::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 <class T>
-Fortran_Array3D<T> & Fortran_Array3D<T>::ref(const Fortran_Array3D<T> &A)
-{
-
- if (this != &A)
- {
- v_ = A.v_;
- m_ = A.m_;
- n_ = A.n_;
- k_ = A.k_;
- data_ = A.data_;
- }
- return *this;
-}
-
-template <class T>
-Fortran_Array3D<T> & Fortran_Array3D<T>::operator=(const Fortran_Array3D<T> &A)
-{
- return ref(A);
-}
-
-template <class T>
-inline int Fortran_Array3D<T>::dim1() const { return m_; }
-
-template <class T>
-inline int Fortran_Array3D<T>::dim2() const { return n_; }
-
-template <class T>
-inline int Fortran_Array3D<T>::dim3() const { return k_; }
-
-
-template <class T>
-inline int Fortran_Array3D<T>::ref_count() const
-{
- return v_.ref_count();
-}
-
-template <class T>
-Fortran_Array3D<T>::~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 <cstdlib>
-#include <cassert>
-
-namespace TNT
-{
-
-
-template <class T>
-std::ostream& operator<<(std::ostream &s, const Fortran_Array3D<T> &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 <class T>
-std::istream& operator>>(std::istream &s, Fortran_Array3D<T> &A)
-{
-
- int M, N, K;
-
- s >> M >> N >> K;
-
- Fortran_Array3D<T> 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 <class T>
-Fortran_Array3D<T> operator+(const Fortran_Array3D<T> &A, const Fortran_Array3D<T> &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<T>();
-
- else
- {
- Fortran_Array3D<T> 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 <class T>
-Fortran_Array3D<T> operator-(const Fortran_Array3D<T> &A, const Fortran_Array3D<T> &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<T>();
-
- else
- {
- Fortran_Array3D<T> 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 <class T>
-Fortran_Array3D<T> operator*(const Fortran_Array3D<T> &A, const Fortran_Array3D<T> &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<T>();
-
- else
- {
- Fortran_Array3D<T> 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 <class T>
-Fortran_Array3D<T> operator/(const Fortran_Array3D<T> &A, const Fortran_Array3D<T> &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<T>();
-
- else
- {
- Fortran_Array3D<T> 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 <class T>
-Fortran_Array3D<T>& operator+=(Fortran_Array3D<T> &A, const Fortran_Array3D<T> &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 <class T>
-Fortran_Array3D<T>& operator-=(Fortran_Array3D<T> &A, const Fortran_Array3D<T> &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 <class T>
-Fortran_Array3D<T>& operator*=(Fortran_Array3D<T> &A, const Fortran_Array3D<T> &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 <class T>
-Fortran_Array3D<T>& operator/=(Fortran_Array3D<T> &A, const Fortran_Array3D<T> &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 <cstdlib>
-#include <iostream>
-
-#ifdef TNT_BOUNDS_CHECK
-#include <assert.h>
-#endif
-
-#ifndef NULL
-#define NULL 0
-#endif
-
-namespace TNT
-{
-/*
- Internal representation of ref-counted array. The TNT
- arrays all use this building block.
-
- <p>
- 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 T>
-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<T> & operator=(const i_refvec<T> &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 <class T>
-void i_refvec<T>::copy_(T* p, const T* q, const T* e)
-{
- for (T* t=p; q<e; t++, q++)
- *t= *q;
-}
-
-template <class T>
-i_refvec<T>::i_refvec() : data_(NULL), ref_count_(NULL) {}
-
-/**
- In case n is 0 or negative, it does NOT call new.
-*/
-template <class T>
-i_refvec<T>::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 <class T>
-inline i_refvec<T>::i_refvec(const i_refvec<T> &V): data_(V.data_),
- ref_count_(V.ref_count_)
-{
- if (V.ref_count_ != NULL)
- (*(V.ref_count_))++;
-}
-
-
-template <class T>
-i_refvec<T>::i_refvec(T* data) : data_(data), ref_count_(NULL) {}
-
-template <class T>
-inline T* i_refvec<T>::begin()
-{
- return data_;
-}
-
-template <class T>
-inline const T& i_refvec<T>::operator[](int i) const
-{
- return data_[i];
-}
-
-template <class T>
-inline T& i_refvec<T>::operator[](int i)
-{
- return data_[i];
-}
-
-
-template <class T>
-inline const T* i_refvec<T>::begin() const
-{
- return data_;
-}
-
-
-
-template <class T>
-i_refvec<T> & i_refvec<T>::operator=(const i_refvec<T> &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 <class T>
-void i_refvec<T>::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<class T>
-int i_refvec<T>::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 <class T>
-int i_refvec<T>::ref_count() const
-{
- if (data_ == NULL)
- return 0;
- else
- return (ref_count_ != NULL ? *ref_count_ : -1) ;
-}
-
-template <class T>
-i_refvec<T>::~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 <cmath>
-
-
-
-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 <class Real>
-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.
-
- <p>
- Index values begin at 0.
-
- <p>
- <b>Storage requirements:</b> 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 T>
-class Sparse_Matrix_CompRow {
-
-private:
- Array1D<T> val_; // data values (nz_ elements)
- Array1D<int> rowptr_; // row_ptr (dim_[0]+1 elements)
- Array1D<int> 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 <class T>
-Sparse_Matrix_CompRow<T>::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 <time.h>
-
-
-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 <cstdlib>
-#include <cassert>
-#include <iostream>
-#include <sstream>
-
-namespace TNT
-{
-
-/**
- <b>[Deprecatred]</b> 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 T>
-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<N4; i+=4)
- {
- v_[i] = v[i];
- v_[i+1] = v[i+1];
- v_[i+2] = v[i+2];
- v_[i+3] = v[i+3];
- }
-
- for (i=N4; i< N; i++)
- v_[i] = v[i];
-#else
-
- for (i=0; i< N; i++)
- v_[i] = v[i];
-#endif
- }
-
- void set(const T& val)
- {
- Subscript N = n_;
- Subscript i;
-
-#ifdef TNT_UNROLL_LOOPS
- Subscript Nmod4 = N & 3;
- Subscript N4 = N - Nmod4;
-
- for (i=0; i<N4; i+=4)
- {
- v_[i] = val;
- v_[i+1] = val;
- v_[i+2] = val;
- v_[i+3] = val;
- }
-
- for (i=N4; i< N; i++)
- v_[i] = val;
-#else
-
- for (i=0; i< N; i++)
- v_[i] = val;
-
-#endif
- }
-
-
-
- void destroy()
- {
- /* do nothing, if no memory has been previously allocated */
- if (v_ == NULL) return ;
-
- /* if we are here, then matrix was previously allocated */
- delete [] (v_);
-
- v_ = NULL;
- vm1_ = NULL;
- }
-
-
- public:
-
- // access
-
- iterator begin() { return v_;}
- iterator end() { return v_ + n_; }
- const iterator begin() const { return v_;}
- const iterator end() const { return v_ + n_; }
-
- // destructor
-
- ~Vector()
- {
- destroy();
- }
-
- // constructors
-
- Vector() : v_(0), vm1_(0), n_(0) {};
-
- Vector(const Vector<T> &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<N; i++)
- ins >> v_[i];
- }
-
-
- // methods
- //
- Vector<T>& newsize(Subscript N)
- {
- if (n_ == N) return *this;
-
- destroy();
- initialize(N);
-
- return *this;
- }
-
-
- // assignments
- //
- Vector<T>& operator=(const Vector<T> &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<T>& 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 <class T>
-std::ostream& operator<<(std::ostream &s, const Vector<T> &A)
-{
- Subscript N=A.dim();
-
- s << N << "\n";
-
- for (Subscript i=0; i<N; i++)
- s << A[i] << " " << "\n";
- s << "\n";
-
- return s;
-}
-
-template <class T>
-std::istream & operator>>(std::istream &s, Vector<T> &A)
-{
-
- Subscript N;
-
- s >> N;
-
- if ( !(N == A.size() ))
- {
- A.newsize(N);
- }
-
-
- for (Subscript i=0; i<N; i++)
- s >> A[i];
-
-
- return s;
-}
-
-// *******************[ basic matrix algorithms ]***************************
-
-
-template <class T>
-Vector<T> operator+(const Vector<T> &A,
- const Vector<T> &B)
-{
- Subscript N = A.dim();
-
- assert(N==B.dim());
-
- Vector<T> tmp(N);
- Subscript i;
-
- for (i=0; i<N; i++)
- tmp[i] = A[i] + B[i];
-
- return tmp;
-}
-
-template <class T>
-Vector<T> operator-(const Vector<T> &A,
- const Vector<T> &B)
-{
- Subscript N = A.dim();
-
- assert(N==B.dim());
-
- Vector<T> tmp(N);
- Subscript i;
-
- for (i=0; i<N; i++)
- tmp[i] = A[i] - B[i];
-
- return tmp;
-}
-
-template <class T>
-Vector<T> operator*(const Vector<T> &A,
- const Vector<T> &B)
-{
- Subscript N = A.dim();
-
- assert(N==B.dim());
-
- Vector<T> tmp(N);
- Subscript i;
-
- for (i=0; i<N; i++)
- tmp[i] = A[i] * B[i];
-
- return tmp;
-}
-
-
-template <class T>
-T dot_prod(const Vector<T> &A, const Vector<T> &B)
-{
- Subscript N = A.dim();
- assert(N == B.dim());
-
- Subscript i;
- T sum = 0;
-
- for (i=0; i<N; i++)
- sum += A[i] * B[i];
-
- return sum;
-}
-
-} /* namespace TNT */
-
-#endif
-// TNT_VEC_H
diff --git a/lib/include/tnt/tnt_version.h b/lib/include/tnt/tnt_version.h
deleted file mode 100644
index 047e7d3..0000000
--- a/lib/include/tnt/tnt_version.h
+++ /dev/null
@@ -1,39 +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_VERSION_H
-#define TNT_VERSION_H
-
-
-//---------------------------------------------------------------------
-// current version
-//---------------------------------------------------------------------
-
-
-#define TNT_MAJOR_VERSION '1'
-#define TNT_MINOR_VERSION '2'
-#define TNT_SUBMINOR_VERSION '6'
-#define TNT_VERSION_STRING "1.2.6"
-
-
-
-
-
-#endif
-// TNT_VERSION_H
--
cgit v1.2.3