Subversion Repositories gelsvn

Rev

Rev 46 | Rev 382 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 46 Rev 89
1
#ifndef __CGLA_ARITHMATFLOAT_H__
1
#ifndef __CGLA_ARITHMATFLOAT_H__
2
#define __CGLA_ARITHMATFLOAT_H__
2
#define __CGLA_ARITHMATFLOAT_H__
3
 
3
 
4
#include <vector>
4
#include <vector>
5
#include <iostream>
5
#include <iostream>
6
 
6
 
7
#include "CGLA.h"
7
#include "CGLA.h"
8
 
8
 
9
 
9
 
10
namespace CGLA 
10
namespace CGLA 
11
{
11
{
12
 
12
 
13
  /** Basic class template for matrices.
13
  /** \brief Basic class template for matrices.
14
		
14
		
15
  In this template a matrix is defined as an array of vectors. This may
15
  In this template a matrix is defined as an array of vectors. This may
16
  not in all cases be the most efficient but it has the advantage that 
16
  not in all cases be the most efficient but it has the advantage that 
17
  it is possible to use the double subscripting notation:
17
  it is possible to use the double subscripting notation:
18
		
18
		
19
  T x = m[i][j]
19
  T x = m[i][j]
20
 
20
 
21
  This template should be used through inheritance just like the 
21
  This template should be used through inheritance just like the 
22
  vector template */
22
  vector template */
23
  template <class VVT, class HVT, class MT, unsigned int ROWS>
23
  template <class VVT, class HVT, class MT, unsigned int ROWS>
24
    class ArithMatFloat
24
    class ArithMatFloat
25
    { 
25
    { 
26
    public:
26
    public:
27
 
27
 
28
      /// Horizontal vector type
28
      /// Horizontal vector type
29
      typedef HVT HVectorType;
29
      typedef HVT HVectorType;
30
 
30
 
31
      /// Vertical vector type
31
      /// Vertical vector type
32
      typedef VVT VVectorType;
32
      typedef VVT VVectorType;
33
 
33
 
34
      /// The type of a matrix element
34
      /// The type of a matrix element
35
      typedef typename HVT::ScalarType ScalarType;
35
      typedef typename HVT::ScalarType ScalarType;
36
		
36
		
37
    protected:
37
    protected:
38
 
38
 
39
      /// The actual contents of the matrix.
39
      /// The actual contents of the matrix.
40
      HVT data[ROWS];
40
      HVT data[ROWS];
41
 
41
 
42
    protected:
42
    protected:
43
 
43
 
44
      /// Construct 0 matrix
44
      /// Construct 0 matrix
45
      ArithMatFloat() 
45
      ArithMatFloat() 
46
	{
46
	{
47
#ifndef NDEBUG
47
#ifndef NDEBUG
48
	  std::fill_n(data, ROWS, HVT(CGLA_INIT_VALUE));
48
	  std::fill_n(data, ROWS, HVT(CGLA_INIT_VALUE));
49
#endif
49
#endif
50
	}
50
	}
51
 
51
 
52
      /// Construct a matrix where all entries are the same.
52
      /// Construct a matrix where all entries are the same.
53
      explicit ArithMatFloat(ScalarType x)
53
      explicit ArithMatFloat(ScalarType x)
54
	{
54
	{
55
	  std::fill_n(data, ROWS, HVT(x));
55
	  std::fill_n(data, ROWS, HVT(x));
56
	}
56
	}
57
 
57
 
58
      /// Construct a matrix where all rows are the same.
58
      /// Construct a matrix where all rows are the same.
59
      explicit ArithMatFloat(HVT _a)
59
      explicit ArithMatFloat(HVT _a)
60
	{
60
	{
61
	  std::fill_n(data, ROWS, _a);
61
	  std::fill_n(data, ROWS, _a);
62
	}
62
	}
63
 
63
 
64
      /// Construct a matrix with two rows.
64
      /// Construct a matrix with two rows.
65
      ArithMatFloat(HVT _a, HVT _b)
65
      ArithMatFloat(HVT _a, HVT _b)
66
	{
66
	{
67
	  assert(ROWS==2);
67
	  assert(ROWS==2);
68
	  data[0] = _a;
68
	  data[0] = _a;
69
	  data[1] = _b;
69
	  data[1] = _b;
70
	}
70
	}
71
 
71
 
72
      /// Construct a matrix with three rows.
72
      /// Construct a matrix with three rows.
73
      ArithMatFloat(HVT _a, HVT _b, HVT _c)
73
      ArithMatFloat(HVT _a, HVT _b, HVT _c)
74
	{
74
	{
75
	  assert(ROWS==3);
75
	  assert(ROWS==3);
76
	  data[0] = _a;
76
	  data[0] = _a;
77
	  data[1] = _b;
77
	  data[1] = _b;
78
	  data[2] = _c;
78
	  data[2] = _c;
79
	}
79
	}
80
 
80
 
81
      /// Construct a matrix with four rows.
81
      /// Construct a matrix with four rows.
82
      ArithMatFloat(HVT _a, HVT _b, HVT _c, HVT _d)
82
      ArithMatFloat(HVT _a, HVT _b, HVT _c, HVT _d)
83
	{
83
	{
84
	  assert(ROWS==4);
84
	  assert(ROWS==4);
85
	  data[0] = _a;
85
	  data[0] = _a;
86
	  data[1] = _b;
86
	  data[1] = _b;
87
	  data[2] = _c;
87
	  data[2] = _c;
88
	  data[3] = _d;
88
	  data[3] = _d;
89
	}
89
	}
90
		
90
		
91
    public:
91
    public:
92
 
92
 
93
      /// Get vertical dimension of matrix 
93
      /// Get vertical dimension of matrix 
94
      static unsigned int get_v_dim() {return VVT::get_dim();}
94
      static unsigned int get_v_dim() {return VVT::get_dim();}
95
 
95
 
96
      /// Get horizontal dimension of matrix
96
      /// Get horizontal dimension of matrix
97
      static unsigned int get_h_dim() {return HVT::get_dim();}
97
      static unsigned int get_h_dim() {return HVT::get_dim();}
98
 
98
 
99
 
99
 
100
      /** Get const pointer to data array.
100
      /** Get const pointer to data array.
101
	  This function may be useful when interfacing with some other API 
101
	  This function may be useful when interfacing with some other API 
102
	  such as OpenGL (TM). */
102
	  such as OpenGL (TM). */
103
      const ScalarType* get() const 
103
      const ScalarType* get() const 
104
	{
104
	{
105
	  return data[0].get();
105
	  return data[0].get();
106
	}
106
	}
107
 
107
 
108
      /** Get pointer to data array.
108
      /** Get pointer to data array.
109
	  This function may be useful when interfacing with some other API 
109
	  This function may be useful when interfacing with some other API 
110
	  such as OpenGL (TM). */
110
	  such as OpenGL (TM). */
111
      ScalarType* get()
111
      ScalarType* get()
112
	{
112
	{
113
	  return data[0].get();
113
	  return data[0].get();
114
	}
114
	}
115
 
115
 
116
      //----------------------------------------------------------------------
116
      //----------------------------------------------------------------------
117
      // index operators
117
      // index operators
118
      //----------------------------------------------------------------------
118
      //----------------------------------------------------------------------
119
 
119
 
120
      /// Const index operator. Returns i'th row of matrix.
120
      /// Const index operator. Returns i'th row of matrix.
121
      const HVT& operator [] ( unsigned int i ) const
121
      const HVT& operator [] ( unsigned int i ) const
122
	{
122
	{
123
	  assert(i<ROWS);
123
	  assert(i<ROWS);
124
	  return data[i];
124
	  return data[i];
125
	}
125
	}
126
 
126
 
127
      /// Non-const index operator. Returns i'th row of matrix.
127
      /// Non-const index operator. Returns i'th row of matrix.
128
      HVT& operator [] ( unsigned int i ) 
128
      HVT& operator [] ( unsigned int i ) 
129
	{
129
	{
130
	  assert(i<ROWS);
130
	  assert(i<ROWS);
131
	  return data[i];
131
	  return data[i];
132
	}
132
	}
133
 
133
 
134
      //----------------------------------------------------------------------
134
      //----------------------------------------------------------------------
135
 
135
 
136
      /// Equality operator. 
136
      /// Equality operator. 
137
      bool operator==(const MT& v) const 
137
      bool operator==(const MT& v) const 
138
	{
138
	{
139
	  return std::inner_product(data, &data[ROWS], &v[0], true,
139
	  return std::inner_product(data, &data[ROWS], &v[0], true,
140
				    std::logical_and<bool>(), std::equal_to<HVT>());
140
				    std::logical_and<bool>(), std::equal_to<HVT>());
141
	}
141
	}
142
 
142
 
143
      /// Inequality operator.
143
      /// Inequality operator.
144
      bool operator!=(const MT& v) const 
144
      bool operator!=(const MT& v) const 
145
	{
145
	{
146
	  return !(*this==v);
146
	  return !(*this==v);
147
	}
147
	}
148
 
148
 
149
      //----------------------------------------------------------------------
149
      //----------------------------------------------------------------------
150
 
150
 
151
      /// Multiply scalar onto matrix. All entries are multiplied by scalar.
151
      /// Multiply scalar onto matrix. All entries are multiplied by scalar.
152
      const MT operator * (ScalarType k) const
152
      const MT operator * (ScalarType k) const
153
	{
153
	{
154
	  MT v_new;
154
	  MT v_new;
155
	  std::transform(data, &data[ROWS], &v_new[0], std::bind2nd(std::multiplies<HVT>(), k));
155
	  std::transform(data, &data[ROWS], &v_new[0], std::bind2nd(std::multiplies<HVT>(), k));
156
	  return v_new;
156
	  return v_new;
157
	}
157
	}
158
 
158
 
159
      /// Divide all entries in matrix by scalar.
159
      /// Divide all entries in matrix by scalar.
160
      const MT operator / (ScalarType k) const
160
      const MT operator / (ScalarType k) const
161
	{
161
	{
162
	  MT v_new;
162
	  MT v_new;
163
	  std::transform(data, &data[ROWS], &v_new[0], std::bind2nd(std::divides<HVT>(), k));
163
	  std::transform(data, &data[ROWS], &v_new[0], std::bind2nd(std::divides<HVT>(), k));
164
	  return v_new;      
164
	  return v_new;      
165
	}
165
	}
166
 
166
 
167
      /// Assignment multiplication of matrix by scalar.
167
      /// Assignment multiplication of matrix by scalar.
168
      const MT& operator *=(ScalarType k) 
168
      const MT& operator *=(ScalarType k) 
169
	{
169
	{
170
	  std::transform(data, &data[ROWS], data, std::bind2nd(std::multiplies<HVT>(), k));
170
	  std::transform(data, &data[ROWS], data, std::bind2nd(std::multiplies<HVT>(), k));
171
	  return static_cast<const MT&>(*this);
171
	  return static_cast<const MT&>(*this);
172
	}
172
	}
173
 
173
 
174
      /// Assignment division of matrix by scalar.
174
      /// Assignment division of matrix by scalar.
175
      const MT& operator /=(ScalarType k) 
175
      const MT& operator /=(ScalarType k) 
176
	{ 
176
	{ 
177
	  std::transform(data, &data[ROWS], data, std::bind2nd(std::divides<HVT>(), k));
177
	  std::transform(data, &data[ROWS], data, std::bind2nd(std::divides<HVT>(), k));
178
	  return static_cast<const MT&>(*this);
178
	  return static_cast<const MT&>(*this);
179
	}
179
	}
180
 
180
 
181
      //----------------------------------------------------------------------
181
      //----------------------------------------------------------------------
182
 
182
 
183
      /// Add two matrices. 
183
      /// Add two matrices. 
184
      const MT operator + (const MT& m1) const
184
      const MT operator + (const MT& m1) const
185
	{
185
	{
186
	  MT v_new;
186
	  MT v_new;
187
	  std::transform(data, &data[ROWS], &m1[0], &v_new[0], std::plus<HVT>());
187
	  std::transform(data, &data[ROWS], &m1[0], &v_new[0], std::plus<HVT>());
188
	  return v_new;
188
	  return v_new;
189
	}
189
	}
190
 
190
 
191
      /// Subtract two matrices.
191
      /// Subtract two matrices.
192
      const MT operator - (const MT& m1) const
192
      const MT operator - (const MT& m1) const
193
	{
193
	{
194
	  MT v_new;
194
	  MT v_new;
195
	  std::transform(data, &data[ROWS], &m1[0], &v_new[0], std::minus<HVT>());
195
	  std::transform(data, &data[ROWS], &m1[0], &v_new[0], std::minus<HVT>());
196
	  return v_new;
196
	  return v_new;
197
	}
197
	}
198
 
198
 
199
      /// Assigment addition of matrices.
199
      /// Assigment addition of matrices.
200
      const MT& operator +=(const MT& v) 
200
      const MT& operator +=(const MT& v) 
201
	{
201
	{
202
	  std::transform(data, &data[ROWS], &v[0], data, std::plus<HVT>());
202
	  std::transform(data, &data[ROWS], &v[0], data, std::plus<HVT>());
203
	  return static_cast<const MT&>(*this);
203
	  return static_cast<const MT&>(*this);
204
	}
204
	}
205
 
205
 
206
      /// Assigment subtraction of matrices.
206
      /// Assigment subtraction of matrices.
207
      const MT& operator -=(const MT& v) 
207
      const MT& operator -=(const MT& v) 
208
	{
208
	{
209
	  std::transform(data, &data[ROWS], &v[0], data, std::minus<HVT>());
209
	  std::transform(data, &data[ROWS], &v[0], data, std::minus<HVT>());
210
	  return static_cast<const MT&>(*this);
210
	  return static_cast<const MT&>(*this);
211
	}
211
	}
212
 
212
 
213
      //----------------------------------------------------------------------
213
      //----------------------------------------------------------------------
214
 
214
 
215
      /// Negate matrix.
215
      /// Negate matrix.
216
      const MT operator - () const
216
      const MT operator - () const
217
	{
217
	{
218
	  MT v_new;
218
	  MT v_new;
219
	  std::transform(data, &data[ROWS], &v_new[0], std::negate<HVT>());
219
	  std::transform(data, &data[ROWS], &v_new[0], std::negate<HVT>());
220
	  return v_new;
220
	  return v_new;
221
	}
221
	}
222
    };
222
    };
223
 
223
 
224
  /// Multiply scalar onto matrix
224
  /// Multiply scalar onto matrix
225
  template <class VVT, class HVT, class MT, unsigned int ROWS>
225
  template <class VVT, class HVT, class MT, unsigned int ROWS>
226
    inline const MT operator * (double k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
226
    inline const MT operator * (double k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
227
    {
227
    {
228
      return v * k;
228
      return v * k;
229
    }
229
    }
230
 
230
 
231
  /// Multiply scalar onto matrix
231
  /// Multiply scalar onto matrix
232
  template <class VVT, class HVT, class MT, unsigned int ROWS>
232
  template <class VVT, class HVT, class MT, unsigned int ROWS>
233
    inline const MT operator * (float k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
233
    inline const MT operator * (float k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
234
    {
234
    {
235
      return v * k;
235
      return v * k;
236
    }
236
    }
237
 
237
 
238
  /// Multiply scalar onto matrix
238
  /// Multiply scalar onto matrix
239
  template <class VVT, class HVT, class MT, unsigned int ROWS>
239
  template <class VVT, class HVT, class MT, unsigned int ROWS>
240
    inline const MT operator * (int k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
240
    inline const MT operator * (int k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
241
    {
241
    {
242
      return v * k;
242
      return v * k;
243
    }
243
    }
244
 
244
 
245
  /// Multiply vector onto matrix 
245
  /// Multiply vector onto matrix 
246
  template <class VVT, class HVT, class MT, unsigned int ROWS>
246
  template <class VVT, class HVT, class MT, unsigned int ROWS>
247
    inline VVT operator*(const ArithMatFloat<VVT,HVT,MT,ROWS>& m,const HVT& v) 
247
    inline VVT operator*(const ArithMatFloat<VVT,HVT,MT,ROWS>& m,const HVT& v) 
248
    {
248
    {
249
      VVT v2;
249
      VVT v2;
250
      for(unsigned int i=0;i<ROWS;i++) v2[i] = dot(m[i], v);
250
      for(unsigned int i=0;i<ROWS;i++) v2[i] = dot(m[i], v);
251
      return v2;
251
      return v2;
252
    }
252
    }
253
 
253
 
254
 
254
 
255
#ifndef WIN32
255
#ifndef WIN32
256
  /** Multiply two arbitrary matrices. 
256
  /** Multiply two arbitrary matrices. 
257
      In principle, this function could return a matrix, but in general
257
      In principle, this function could return a matrix, but in general
258
      the new matrix will be of a type that is different from either of
258
      the new matrix will be of a type that is different from either of
259
      the two matrices that are multiplied together. We do not want to 
259
      the two matrices that are multiplied together. We do not want to 
260
      return an ArithMatFloat - so it seems best to let the return value be
260
      return an ArithMatFloat - so it seems best to let the return value be
261
      a reference arg.
261
      a reference arg.
262
		
262
		
263
      This template can only be instantiated if the dimensions of the
263
      This template can only be instantiated if the dimensions of the
264
      matrices match -- i.e. if the multiplication can actually be
264
      matrices match -- i.e. if the multiplication can actually be
265
      carried out. This is more type safe than the win32 version below.
265
      carried out. This is more type safe than the win32 version below.
266
  */
266
  */
267
 
267
 
268
  template <class VVT, class HVT, 
268
  template <class VVT, class HVT, 
269
    class HV1T, class VV2T,
269
    class HV1T, class VV2T,
270
    class MT1, class MT2, class MT,
270
    class MT1, class MT2, class MT,
271
    unsigned int ROWS1, unsigned int ROWS2>
271
    unsigned int ROWS1, unsigned int ROWS2>
272
    inline void mul(const ArithMatFloat<VVT,HV1T,MT1,ROWS1>& m1,
272
    inline void mul(const ArithMatFloat<VVT,HV1T,MT1,ROWS1>& m1,
273
		    const ArithMatFloat<VV2T,HVT,MT2,ROWS2>& m2,
273
		    const ArithMatFloat<VV2T,HVT,MT2,ROWS2>& m2,
274
		    ArithMatFloat<VVT,HVT,MT,ROWS1>& m)
274
		    ArithMatFloat<VVT,HVT,MT,ROWS1>& m)
275
    {
275
    {
276
      unsigned int cols = ArithMatFloat<VVT,HVT,MT,ROWS1>::get_h_dim();
276
      unsigned int cols = ArithMatFloat<VVT,HVT,MT,ROWS1>::get_h_dim();
277
      for(unsigned int i=0;i<ROWS1;i++)
277
      for(unsigned int i=0;i<ROWS1;i++)
278
	for(unsigned int j=0;j<cols;j++)
278
	for(unsigned int j=0;j<cols;j++)
279
	  {
279
	  {
280
	    m[i][j] = 0;
280
	    m[i][j] = 0;
281
	    for(unsigned int k=0;k<ROWS2;k++)
281
	    for(unsigned int k=0;k<ROWS2;k++)
282
	      m[i][j] += m1[i][k] * m2[k][j]; 
282
	      m[i][j] += m1[i][k] * m2[k][j]; 
283
	  }
283
	  }
284
    }
284
    }
285
 
285
 
286
 
286
 
287
  /** Transpose. See the discussion on mul if you are curious as to why
287
  /** Transpose. See the discussion on mul if you are curious as to why
288
      I don't simply return the transpose. */
288
      I don't simply return the transpose. */
289
  template <class VVT, class HVT, class M1T, class M2T, unsigned int ROWS, unsigned int COLS>
289
  template <class VVT, class HVT, class M1T, class M2T, unsigned int ROWS, unsigned int COLS>
290
    inline void transpose(const ArithMatFloat<VVT,HVT,M1T,ROWS>& m,
290
    inline void transpose(const ArithMatFloat<VVT,HVT,M1T,ROWS>& m,
291
			  ArithMatFloat<HVT,VVT,M2T,COLS>& m_new)
291
			  ArithMatFloat<HVT,VVT,M2T,COLS>& m_new)
292
    {
292
    {
293
      for(unsigned int i=0;i<M2T::get_v_dim();++i)
293
      for(unsigned int i=0;i<M2T::get_v_dim();++i)
294
	for(unsigned int j=0;j<M2T::get_h_dim();++j)
294
	for(unsigned int j=0;j<M2T::get_h_dim();++j)
295
	  m_new[i][j] = m[j][i];
295
	  m_new[i][j] = m[j][i];
296
    }
296
    }
297
 
297
 
298
#else
298
#else
299
 
299
 
300
  //----------------- win32 -------------------------------
300
  //----------------- win32 -------------------------------
301
  // Visual studio is not good at deducing the args. to these template functions.
301
  // Visual studio is not good at deducing the args. to these template functions.
302
  // This means that you can call the two functions below with 
302
  // This means that you can call the two functions below with 
303
  // matrices of wrong dimension.
303
  // matrices of wrong dimension.
304
 
304
 
305
  template <class M1, class M2, class M>
305
  template <class M1, class M2, class M>
306
    inline void mul(const M1& m1, const M2& m2, M& m)
306
    inline void mul(const M1& m1, const M2& m2, M& m)
307
    {
307
    {
308
      unsigned int cols = M::get_h_dim();
308
      unsigned int cols = M::get_h_dim();
309
      unsigned int rows1 = M1::get_v_dim();
309
      unsigned int rows1 = M1::get_v_dim();
310
      unsigned int rows2 = M2::get_v_dim();
310
      unsigned int rows2 = M2::get_v_dim();
311
 
311
 
312
      for(unsigned int i=0;i<rows1;++i)
312
      for(unsigned int i=0;i<rows1;++i)
313
	for(unsigned int j=0;j<cols;++j)
313
	for(unsigned int j=0;j<cols;++j)
314
	  {
314
	  {
315
	    m[i][j] = 0;
315
	    m[i][j] = 0;
316
	    for(unsigned int k=0;k<rows2;++k)
316
	    for(unsigned int k=0;k<rows2;++k)
317
	      m[i][j] += m1[i][k] * m2[k][j];
317
	      m[i][j] += m1[i][k] * m2[k][j];
318
	  }
318
	  }
319
    }
319
    }
320
 
320
 
321
 
321
 
322
  /** Transpose. See the discussion on mul if you are curious as to why
322
  /** Transpose. See the discussion on mul if you are curious as to why
323
      I don't simply return the transpose. */
323
      I don't simply return the transpose. */
324
  template <class M1, class M2>
324
  template <class M1, class M2>
325
    inline void transpose(const M1& m1, M2& m2)
325
    inline void transpose(const M1& m1, M2& m2)
326
    {
326
    {
327
      for(unsigned int i=0;i<M2::get_v_dim();++i)
327
      for(unsigned int i=0;i<M2::get_v_dim();++i)
328
	for(unsigned int j=0;j<M2::get_h_dim();++j)
328
	for(unsigned int j=0;j<M2::get_h_dim();++j)
329
	  m2[i][j] = m1[j][i];
329
	  m2[i][j] = m1[j][i];
330
    }
330
    }
331
 
331
 
332
#endif
332
#endif
333
 
333
 
334
  /** Compute the outer product of a and b: a * transpose(b). This is 
334
  /** Compute the outer product of a and b: a * transpose(b). This is 
335
      a matrix with a::rows and b::columns. */
335
      a matrix with a::rows and b::columns. */
336
  template <class VVT, class HVT, class MT, unsigned int ROWS>
336
  template <class VVT, class HVT, class MT, unsigned int ROWS>
337
    void outer_product(const VVT& a, const HVT& b, 
337
    void outer_product(const VVT& a, const HVT& b, 
338
		       ArithMatFloat<VVT,HVT,MT,ROWS>& m)
338
		       ArithMatFloat<VVT,HVT,MT,ROWS>& m)
339
    {
339
    {
340
      unsigned int R = VVT::get_dim();
340
      unsigned int R = VVT::get_dim();
341
      unsigned int C = HVT::get_dim();
341
      unsigned int C = HVT::get_dim();
342
      for(unsigned int i=0;i<R;++i)
342
      for(unsigned int i=0;i<R;++i)
343
	for(unsigned int j=0;j<C;++j)
343
	for(unsigned int j=0;j<C;++j)
344
	  {
344
	  {
345
	    m[i][j] = a[i] * b[j];
345
	    m[i][j] = a[i] * b[j];
346
	  }
346
	  }
347
    }
347
    }
348
 
348
 
349
  /** Compute the outer product of a and b using an arbitrary 
349
  /** Compute the outer product of a and b using an arbitrary 
350
      binary operation: op(a, transpose(b)). This is 
350
      binary operation: op(a, transpose(b)). This is 
351
      a matrix with a::rows and b::columns. */
351
      a matrix with a::rows and b::columns. */
352
  template <class VVT, class HVT, class MT, int ROWS, class BinOp>
352
  template <class VVT, class HVT, class MT, int ROWS, class BinOp>
353
    void outer_product(const VVT& a, const HVT& b, 
353
    void outer_product(const VVT& a, const HVT& b, 
354
		       ArithMatFloat<VVT,HVT,MT,ROWS>& m, BinOp op)
354
		       ArithMatFloat<VVT,HVT,MT,ROWS>& m, BinOp op)
355
    {
355
    {
356
      int R = VVT::get_dim();
356
      int R = VVT::get_dim();
357
      int C = HVT::get_dim();
357
      int C = HVT::get_dim();
358
      for(int i=0;i<R;++i)
358
      for(int i=0;i<R;++i)
359
	for(int j=0;j<C;++j)
359
	for(int j=0;j<C;++j)
360
	  {
360
	  {
361
	    m[i][j] = op(a[i], b[j]);
361
	    m[i][j] = op(a[i], b[j]);
362
	  }
362
	  }
363
    }
363
    }
364
 
364
 
365
  /** Copy a matrix to another matrix, cell by cell.
365
  /** Copy a matrix to another matrix, cell by cell.
366
      This conversion that takes a const matrix as first argument
366
      This conversion that takes a const matrix as first argument
367
      (source) and a non-const matrix as second argument
367
      (source) and a non-const matrix as second argument
368
      (destination). The contents of the first matrix is simply copied
368
      (destination). The contents of the first matrix is simply copied
369
      to the second matrix. 
369
      to the second matrix. 
370
			
370
			
371
      However, if the first matrix is	larger than the second,
371
      However, if the first matrix is	larger than the second,
372
      the cells outside the range of the destination are simply not
372
      the cells outside the range of the destination are simply not
373
      copied. If the destination is larger, the cells outside the 
373
      copied. If the destination is larger, the cells outside the 
374
      range of the source matrix are not touched.
374
      range of the source matrix are not touched.
375
 
375
 
376
      An obvious use of this function is to copy a 3x3 rotation matrix
376
      An obvious use of this function is to copy a 3x3 rotation matrix
377
      into a 4x4 transformation matrix.
377
      into a 4x4 transformation matrix.
378
  */
378
  */
379
 
379
 
380
  template <class M1, class M2>
380
  template <class M1, class M2>
381
    void copy_matrix(const M1& inmat, M2& outmat)
381
    void copy_matrix(const M1& inmat, M2& outmat)
382
    {
382
    {
383
      const unsigned int R = s_min(inmat.get_v_dim(), outmat.get_v_dim());
383
      const unsigned int R = s_min(inmat.get_v_dim(), outmat.get_v_dim());
384
      const unsigned int C = s_min(inmat.get_h_dim(), outmat.get_h_dim());
384
      const unsigned int C = s_min(inmat.get_h_dim(), outmat.get_h_dim());
385
      for(unsigned int i=0;i<R;++i)
385
      for(unsigned int i=0;i<R;++i)
386
	for(unsigned int j=0;j<C;++j)
386
	for(unsigned int j=0;j<C;++j)
387
	  outmat[i][j] = inmat[i][j];
387
	  outmat[i][j] = inmat[i][j];
388
    }
388
    }
389
 
389
 
390
  /** Put to operator */
390
  /** Put to operator */
391
  template <class VVT, class HVT, class MT, unsigned int ROWS>
391
  template <class VVT, class HVT, class MT, unsigned int ROWS>
392
    inline std::ostream& 
392
    inline std::ostream& 
393
    operator<<(std::ostream&os, const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
393
    operator<<(std::ostream&os, const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
394
    {
394
    {
395
      os << "[\n";
395
      os << "[\n";
396
      for(unsigned int i=0;i<ROWS;i++) os << "  " << m[i] << "\n";
396
      for(unsigned int i=0;i<ROWS;i++) os << "  " << m[i] << "\n";
397
      os << "]\n";
397
      os << "]\n";
398
      return os;
398
      return os;
399
    }
399
    }
400
 
400
 
401
  /** Get from operator */
401
  /** Get from operator */
402
  template <class VVT, class HVT, class MT, unsigned int ROWS>
402
  template <class VVT, class HVT, class MT, unsigned int ROWS>
403
    inline std::istream& operator>>(std::istream&is, 
403
    inline std::istream& operator>>(std::istream&is, 
404
				    const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
404
				    const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
405
    {
405
    {
406
      for(unsigned int i=0;i<ROWS;i++) is>>m[i];
406
      for(unsigned int i=0;i<ROWS;i++) is>>m[i];
407
      return is;
407
      return is;
408
    }
408
    }
409
}
409
}
410
#endif
410
#endif
411
 
411