diff options
Diffstat (limited to 'lib/include/tnt')
26 files changed, 7184 insertions, 0 deletions
diff --git a/lib/include/tnt/jama_cholesky.h b/lib/include/tnt/jama_cholesky.h new file mode 100644 index 0000000..371d2f7 --- /dev/null +++ b/lib/include/tnt/jama_cholesky.h @@ -0,0 +1,258 @@ +#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 new file mode 100644 index 0000000..8e3572f --- /dev/null +++ b/lib/include/tnt/jama_eig.h @@ -0,0 +1,1034 @@ +#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 new file mode 100644 index 0000000..e95b433 --- /dev/null +++ b/lib/include/tnt/jama_lu.h @@ -0,0 +1,323 @@ +#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 new file mode 100644 index 0000000..98d37f5 --- /dev/null +++ b/lib/include/tnt/jama_qr.h @@ -0,0 +1,326 @@ +#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 new file mode 100644 index 0000000..72ce3a7 --- /dev/null +++ b/lib/include/tnt/jama_svd.h @@ -0,0 +1,543 @@ +#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 new file mode 100644 index 0000000..92463e0 --- /dev/null +++ b/lib/include/tnt/tnt.h @@ -0,0 +1,64 @@ +/* +* +* 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 new file mode 100644 index 0000000..858df57 --- /dev/null +++ b/lib/include/tnt/tnt_array1d.h @@ -0,0 +1,278 @@ +/* +* +* 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 new file mode 100644 index 0000000..683e0e2 --- /dev/null +++ b/lib/include/tnt/tnt_array1d_utils.h @@ -0,0 +1,230 @@ +/* +* +* 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 new file mode 100644 index 0000000..c791575 --- /dev/null +++ b/lib/include/tnt/tnt_array2d.h @@ -0,0 +1,315 @@ +/* +* +* 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 new file mode 100644 index 0000000..7041ed3 --- /dev/null +++ b/lib/include/tnt/tnt_array2d_utils.h @@ -0,0 +1,287 @@ +/* +* +* 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 new file mode 100644 index 0000000..c210d2e --- /dev/null +++ b/lib/include/tnt/tnt_array3d.h @@ -0,0 +1,296 @@ +/* +* +* 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 new file mode 100644 index 0000000..5acdc1d --- /dev/null +++ b/lib/include/tnt/tnt_array3d_utils.h @@ -0,0 +1,236 @@ + + +#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 new file mode 100644 index 0000000..5ff4c48 --- /dev/null +++ b/lib/include/tnt/tnt_cmat.h @@ -0,0 +1,580 @@ +/* +* +* 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 new file mode 100644 index 0000000..ad3bba0 --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array1d.h @@ -0,0 +1,267 @@ +/* +* +* 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 new file mode 100644 index 0000000..b037b17 --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array1d_utils.h @@ -0,0 +1,242 @@ +/* +* +* 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 new file mode 100644 index 0000000..f307536 --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array2d.h @@ -0,0 +1,225 @@ +/* +* +* 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 new file mode 100644 index 0000000..bb68673 --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array2d_utils.h @@ -0,0 +1,236 @@ +/* +* +* 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 new file mode 100644 index 0000000..e51affb --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array3d.h @@ -0,0 +1,223 @@ +/* +* +* 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 new file mode 100644 index 0000000..a13a275 --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array3d_utils.h @@ -0,0 +1,249 @@ +/* +* +* 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 new file mode 100644 index 0000000..5a67eb5 --- /dev/null +++ b/lib/include/tnt/tnt_i_refvec.h @@ -0,0 +1,243 @@ +/* +* +* 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 new file mode 100644 index 0000000..f9c1c91 --- /dev/null +++ b/lib/include/tnt/tnt_math_utils.h @@ -0,0 +1,34 @@ +#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 new file mode 100644 index 0000000..0d4fde1 --- /dev/null +++ b/lib/include/tnt/tnt_sparse_matrix_csr.h @@ -0,0 +1,103 @@ +/* +* +* 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 new file mode 100644 index 0000000..8dc5d23 --- /dev/null +++ b/lib/include/tnt/tnt_stopwatch.h @@ -0,0 +1,95 @@ +/* +* +* 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 new file mode 100644 index 0000000..d8fe120 --- /dev/null +++ b/lib/include/tnt/tnt_subscript.h @@ -0,0 +1,54 @@ +/* +* +* 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 new file mode 100644 index 0000000..3455d79 --- /dev/null +++ b/lib/include/tnt/tnt_vec.h @@ -0,0 +1,404 @@ +/* +* +* 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 new file mode 100644 index 0000000..047e7d3 --- /dev/null +++ b/lib/include/tnt/tnt_version.h @@ -0,0 +1,39 @@ +/* +* +* 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 |