Subversion Repositories gelsvn

Rev

Rev 630 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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