Index: src/sas/simulation/iqPy/libiqPy/tnt/tnt_array2d_utils.h
===================================================================
--- src/sas/simulation/iqPy/libiqPy/tnt/tnt_array2d_utils.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/simulation/iqPy/libiqPy/tnt/tnt_array2d_utils.h	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -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
