src/math/tnt_array2d_utils.hpp

/*
*
* 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
Revision:	1336
Committed:	Tue Apr 7 20:16:26 2009 UTC (16 years ago) by gezelter
File size:	4613 byte(s)
Log Message:	adding jama and tnt libraries for various linear algebra routines
#	Content
1	/*
2	*
3	* Template Numerical Toolkit (TNT)
4	*
5	* Mathematical and Computational Sciences Division
6	* National Institute of Technology,
7	* Gaithersburg, MD USA
8	*
9	*
10	* This software was developed at the National Institute of Standards and
11	* Technology (NIST) by employees of the Federal Government in the course
12	* of their official duties. Pursuant to title 17 Section 105 of the
13	* United States Code, this software is not subject to copyright protection
14	* and is in the public domain. NIST assumes no responsibility whatsoever for
15	* its use by other parties, and makes no guarantees, expressed or implied,
16	* about its quality, reliability, or any other characteristic.
17	*
18	*/
19
20
21	#ifndef TNT_ARRAY2D_UTILS_H
22	#define TNT_ARRAY2D_UTILS_H
23
24	#include <cstdlib>
25	#include <cassert>
26
27	namespace TNT
28	{
29
30
31	template <class T>
32	std::ostream& operator<<(std::ostream &s, const Array2D<T> &A)
33	{
34	int M=A.dim1();
35	int N=A.dim2();
36
37	s << M << " " << N << "\n";
38
39	for (int i=0; i<M; i++)
40	{
41	for (int j=0; j<N; j++)
42	{
43	s << A[i][j] << " ";
44	}
45	s << "\n";
46	}
47
48
49	return s;
50	}
51
52	template <class T>
53	std::istream& operator>>(std::istream &s, Array2D<T> &A)
54	{
55
56	int M, N;
57
58	s >> M >> N;
59
60	Array2D<T> B(M,N);
61
62	for (int i=0; i<M; i++)
63	for (int j=0; j<N; j++)
64	{
65	s >> B[i][j];
66	}
67
68	A = B;
69	return s;
70	}
71
72
73	template <class T>
74	Array2D<T> operator+(const Array2D<T> &A, const Array2D<T> &B)
75	{
76	int m = A.dim1();
77	int n = A.dim2();
78
79	if (B.dim1() != m \|\| B.dim2() != n )
80	return Array2D<T>();
81
82	else
83	{
84	Array2D<T> C(m,n);
85
86	for (int i=0; i<m; i++)
87	{
88	for (int j=0; j<n; j++)
89	C[i][j] = A[i][j] + B[i][j];
90	}
91	return C;
92	}
93	}
94
95	template <class T>
96	Array2D<T> operator-(const Array2D<T> &A, const Array2D<T> &B)
97	{
98	int m = A.dim1();
99	int n = A.dim2();
100
101	if (B.dim1() != m \|\| B.dim2() != n )
102	return Array2D<T>();
103
104	else
105	{
106	Array2D<T> C(m,n);
107
108	for (int i=0; i<m; i++)
109	{
110	for (int j=0; j<n; j++)
111	C[i][j] = A[i][j] - B[i][j];
112	}
113	return C;
114	}
115	}
116
117
118	template <class T>
119	Array2D<T> operator*(const Array2D<T> &A, const Array2D<T> &B)
120	{
121	int m = A.dim1();
122	int n = A.dim2();
123
124	if (B.dim1() != m \|\| B.dim2() != n )
125	return Array2D<T>();
126
127	else
128	{
129	Array2D<T> C(m,n);
130
131	for (int i=0; i<m; i++)
132	{
133	for (int j=0; j<n; j++)
134	C[i][j] = A[i][j] * B[i][j];
135	}
136	return C;
137	}
138	}
139
140
141
142
143	template <class T>
144	Array2D<T> operator/(const Array2D<T> &A, const Array2D<T> &B)
145	{
146	int m = A.dim1();
147	int n = A.dim2();
148
149	if (B.dim1() != m \|\| B.dim2() != n )
150	return Array2D<T>();
151
152	else
153	{
154	Array2D<T> C(m,n);
155
156	for (int i=0; i<m; i++)
157	{
158	for (int j=0; j<n; j++)
159	C[i][j] = A[i][j] / B[i][j];
160	}
161	return C;
162	}
163	}
164
165
166
167
168
169	template <class T>
170	Array2D<T>& operator+=(Array2D<T> &A, const Array2D<T> &B)
171	{
172	int m = A.dim1();
173	int n = A.dim2();
174
175	if (B.dim1() == m \|\| B.dim2() == n )
176	{
177	for (int i=0; i<m; i++)
178	{
179	for (int j=0; j<n; j++)
180	A[i][j] += B[i][j];
181	}
182	}
183	return A;
184	}
185
186
187
188	template <class T>
189	Array2D<T>& operator-=(Array2D<T> &A, const Array2D<T> &B)
190	{
191	int m = A.dim1();
192	int n = A.dim2();
193
194	if (B.dim1() == m \|\| B.dim2() == n )
195	{
196	for (int i=0; i<m; i++)
197	{
198	for (int j=0; j<n; j++)
199	A[i][j] -= B[i][j];
200	}
201	}
202	return A;
203	}
204
205
206
207	template <class T>
208	Array2D<T>& operator*=(Array2D<T> &A, const Array2D<T> &B)
209	{
210	int m = A.dim1();
211	int n = A.dim2();
212
213	if (B.dim1() == m \|\| B.dim2() == n )
214	{
215	for (int i=0; i<m; i++)
216	{
217	for (int j=0; j<n; j++)
218	A[i][j] *= B[i][j];
219	}
220	}
221	return A;
222	}
223
224
225
226
227
228	template <class T>
229	Array2D<T>& operator/=(Array2D<T> &A, const Array2D<T> &B)
230	{
231	int m = A.dim1();
232	int n = A.dim2();
233
234	if (B.dim1() == m \|\| B.dim2() == n )
235	{
236	for (int i=0; i<m; i++)
237	{
238	for (int j=0; j<n; j++)
239	A[i][j] /= B[i][j];
240	}
241	}
242	return A;
243	}
244
245	/**
246	Matrix Multiply: compute C = A*B, where C[i][j]
247	is the dot-product of row i of A and column j of B.
248
249
250	@param A an (m x n) array
251	@param B an (n x k) array
252	@return the (m x k) array A*B, or a null array (0x0)
253	if the matrices are non-conformant (i.e. the number
254	of columns of A are different than the number of rows of B.)
255
256
257	*/
258	template <class T>
259	Array2D<T> matmult(const Array2D<T> &A, const Array2D<T> &B)
260	{
261	if (A.dim2() != B.dim1())
262	return Array2D<T>();
263
264	int M = A.dim1();
265	int N = A.dim2();
266	int K = B.dim2();
267
268	Array2D<T> C(M,K);
269
270	for (int i=0; i<M; i++)
271	for (int j=0; j<K; j++)
272	{
273	T sum = 0;
274
275	for (int k=0; k<N; k++)
276	sum += A[i][k] * B [k][j];
277
278	C[i][j] = sum;
279	}
280
281	return C;
282
283	}
284
285	} // namespace TNT
286
287	#endif