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
#	User	Rev	Content
1	gezelter	1336	/*
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