Subversion Repositories gelsvn

Rev

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

Rev 43 Rev 45
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
  /** 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, 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 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 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 [] ( 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 [] ( 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, 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, 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, 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, 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(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
    int ROWS1, 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
      int cols = ArithMatFloat<VVT,HVT,MT,ROWS1>::get_h_dim();
276
      unsigned int cols = ArithMatFloat<VVT,HVT,MT,ROWS1>::get_h_dim();
277
      for(int i=0;i<ROWS1;i++)
277
      for(unsigned int i=0;i<ROWS1;i++)
278
	for(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(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, int ROWS, 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(int i=0;i<M2T::get_v_dim();++i)
293
      for(unsigned int i=0;i<M2T::get_v_dim();++i)
294
	for(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
      int cols = M::get_h_dim();
308
      unsigned int cols = M::get_h_dim();
309
      int rows1 = M1::get_v_dim();
309
      unsigned int rows1 = M1::get_v_dim();
310
      int rows2 = M2::get_v_dim();
310
      unsigned int rows2 = M2::get_v_dim();
311
 
311
 
312
      for(int i=0;i<rows1;++i)
312
      for(unsigned int i=0;i<rows1;++i)
313
	for(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(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(int i=0;i<M2::get_v_dim();++i)
327
      for(unsigned int i=0;i<M2::get_v_dim();++i)
328
	for(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, 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
      int R = VVT::get_dim();
340
      unsigned int R = VVT::get_dim();
341
      int C = HVT::get_dim();
341
      unsigned int C = HVT::get_dim();
342
      for(int i=0;i<R;++i)
342
      for(unsigned int i=0;i<R;++i)
343
	for(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
  /** Copy a matrix to another matrix, cell by cell.
349
  /** Copy a matrix to another matrix, cell by cell.
350
      This conversion that takes a const matrix as first argument
350
      This conversion that takes a const matrix as first argument
351
      (source) and a non-const matrix as second argument
351
      (source) and a non-const matrix as second argument
352
      (destination). The contents of the first matrix is simply copied
352
      (destination). The contents of the first matrix is simply copied
353
      to the second matrix. 
353
      to the second matrix. 
354
			
354
			
355
      However, if the first matrix is	larger than the second,
355
      However, if the first matrix is	larger than the second,
356
      the cells outside the range of the destination are simply not
356
      the cells outside the range of the destination are simply not
357
      copied. If the destination is larger, the cells outside the 
357
      copied. If the destination is larger, the cells outside the 
358
      range of the source matrix are not touched.
358
      range of the source matrix are not touched.
359
 
359
 
360
      An obvious use of this function is to copy a 3x3 rotation matrix
360
      An obvious use of this function is to copy a 3x3 rotation matrix
361
      into a 4x4 transformation matrix.
361
      into a 4x4 transformation matrix.
362
  */
362
  */
363
 
363
 
364
  template <class M1, class M2>
364
  template <class M1, class M2>
365
    void copy_matrix(const M1& inmat, M2& outmat)
365
    void copy_matrix(const M1& inmat, M2& outmat)
366
    {
366
    {
367
      const int R = s_min(inmat.get_v_dim(), outmat.get_v_dim());
367
      const unsigned int R = s_min(inmat.get_v_dim(), outmat.get_v_dim());
368
      const int C = s_min(inmat.get_h_dim(), outmat.get_h_dim());
368
      const unsigned int C = s_min(inmat.get_h_dim(), outmat.get_h_dim());
369
      for(int i=0;i<R;++i)
369
      for(unsigned int i=0;i<R;++i)
370
	for(int j=0;j<C;++j)
370
	for(unsigned int j=0;j<C;++j)
371
	  outmat[i][j] = inmat[i][j];
371
	  outmat[i][j] = inmat[i][j];
372
    }
372
    }
373
 
373
 
374
  /** Put to operator */
374
  /** Put to operator */
375
  template <class VVT, class HVT, class MT, int ROWS>
375
  template <class VVT, class HVT, class MT, unsigned int ROWS>
376
    inline std::ostream& 
376
    inline std::ostream& 
377
    operator<<(std::ostream&os, const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
377
    operator<<(std::ostream&os, const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
378
    {
378
    {
379
      os << "[\n";
379
      os << "[\n";
380
      for(int i=0;i<ROWS;i++) os << "  " << m[i] << "\n";
380
      for(unsigned int i=0;i<ROWS;i++) os << "  " << m[i] << "\n";
381
      os << "]\n";
381
      os << "]\n";
382
      return os;
382
      return os;
383
    }
383
    }
384
 
384
 
385
  /** Get from operator */
385
  /** Get from operator */
386
  template <class VVT, class HVT, class MT, int ROWS>
386
  template <class VVT, class HVT, class MT, unsigned int ROWS>
387
    inline std::istream& operator>>(std::istream&is, 
387
    inline std::istream& operator>>(std::istream&is, 
388
				    const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
388
				    const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
389
    {
389
    {
390
      for(int i=0;i<ROWS;i++) is>>m[i];
390
      for(unsigned int i=0;i<ROWS;i++) is>>m[i];
391
      return is;
391
      return is;
392
    }
392
    }
393
}
393
}
394
#endif
394
#endif
395
 
395