Fawkes API  Fawkes Development Version
matrix.cpp
1 
2 /***************************************************************************
3  * matrix.cpp - A matrix class
4  *
5  * Created: Wed Sep 26 13:54:12 2007
6  * Copyright 2007-2009 Daniel Beck <beck@kbsg.rwth-aachen.de>
7  * 2009 Masrur Doostdar <doostdar@kbsg.rwth-aachen.de>
8  * 2009 Christof Rath <c.rath@student.tugraz.at>
9  *
10  ****************************************************************************/
11 
12 /* This program is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version. A runtime exception applies to
16  * this software (see LICENSE.GPL_WRE file mentioned below for details).
17  *
18  * This program is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU Library General Public License for more details.
22  *
23  * Read the full text in the LICENSE.GPL_WRE file in the doc directory.
24  */
25 
26 #include "matrix.h"
27 #include "vector.h"
28 
29 #include <core/exceptions/software.h>
30 
31 #include <cstdlib>
32 #include <cstdio>
33 #include <algorithm>
34 
35 #ifdef HAVE_OPENCV
36 # include <cstddef>
37 # include <opencv/cv.h>
38 #endif
39 
40 namespace fawkes
41 {
42 
43 /** @class Matrix <geometry/matrix.h>
44  * A general matrix class.
45  * It provides all the operations that are commonly used with a matrix, but has
46  * been optimized with typical robotic applications in mind. That meas especially
47  * that the chose data type is single-precision float and the class has been
48  * optimized for small matrices (up to about 10x10).
49  * @author Daniel Beck
50  * @author Masrur Doostdar
51  * @author Christof Rath
52  */
53 
54 /** @fn inline unsigned int Matrix::num_rows() const
55  * Return the number of rows in the Matrix
56  * @return the number of rows
57  */
58 
59 /** @fn inline unsigned int Matrix::num_cols() const
60  * Return the number of columns in the Matrix
61  * @return the number of columns
62  */
63 
64 /** @fn inline float* Matrix::get_data()
65  * Returns the data pointer
66  * @return the data pointer
67  */
68 
69 /** @fn inline const float* Matrix::get_data() const
70  * Returns the const data pointer
71  * @return the data pointer
72  */
73 
74 /** @fn float Matrix::data(unsigned int row, unsigned int col) const
75  * (Read-only) Access to matrix data without index check.
76  * With this operator it is possible to access a specific
77  * element of the matrix.
78  * Make sure the indices are correct, there is no sanity
79  * check!
80  * @param row the row of the element
81  * @param col the column of the element
82  * @return the value of the specified element
83  */
84 
85  /** @fn float& Matrix::data(unsigned int row, unsigned int col)
86  * (RW) Access to matrix data without index check.
87  * see the read-only access operator for operational details
88  * Make sure the indizes are correct, there is no sanity
89  * check!
90  * @param row the row of the element
91  * @param col the column of the element
92  * @return a reference to the specified element
93  */
94 
95 /** Constructor.
96  * @param num_rows number of rows
97  * @param num_cols number of columns
98  * @param data array containing elements of the matrix in row-by-row-order
99  * @param manage_own_memory if true, the Matrix will manage its memory on its own, else it
100  * will not allocate new memory but works with the provided array
101  */
102 Matrix::Matrix(unsigned int num_rows, unsigned int num_cols,
103  float *data, bool manage_own_memory)
104 {
105  m_num_rows = num_rows;
106  m_num_cols = num_cols;
107  m_num_elements = m_num_rows * m_num_cols;
108 
109  if (!m_num_elements) printf("WTF?\n");
110 
111  if (data == NULL || manage_own_memory)
112  {
113  m_data = (float*) malloc(sizeof(float) * m_num_elements);
114  m_own_memory = true;
115 
116  /* It showed that for arrays up to approx. 1000 elements an optimized for-loop
117  * is faster than a memcpy() call. It is assumed that the same is true for
118  * memset(). */
119  if (data != NULL)
120  {
121  for (unsigned int i = 0; i < m_num_elements; ++i) m_data[i] = data[i];
122  }
123  else
124  {
125  for (unsigned int i = 0; i < m_num_elements; ++i) m_data[i] = 0.f;
126  }
127  }
128  else
129  {
130  m_data = data;
131  m_own_memory = false;
132  }
133 }
134 
135 /** Copy-constructor.
136  * @param tbc matrix to be copied
137  */
139 {
140  m_num_rows = tbc.m_num_rows;
141  m_num_cols = tbc.m_num_cols;
142  m_num_elements = tbc.m_num_elements;
143 
144  m_own_memory = true;
145 
146  m_data = (float*) malloc(sizeof(float) * m_num_elements);
147  for (unsigned int i = 0; i < m_num_elements; ++i) m_data[i] = tbc.m_data[i];
148 }
149 
150 /** Destructor. */
152 {
153  if (m_own_memory) free(m_data);
154 }
155 
156 /** Determines the dimensions of the matrix.
157  * @param num_cols pointer to an unsigned int to where the number of columns is copied to
158  * @param num_rows pointer to an unsigned int to where the number of rows is copied to
159  */
160 void
161 Matrix::size(unsigned int &num_rows, unsigned int &num_cols) const
162 {
163  num_rows = this->num_rows();
164  num_cols = this->num_cols();
165 }
166 
167 
168 /** Sets the diagonal elements to 1.0 and all other to 0.0.
169  * @return a reference to the matrix object
170  */
171 Matrix &
173 {
174  for (unsigned int row = 0; row < num_rows(); ++row)
175  {
176  for (unsigned int col = 0; col < num_cols(); ++col)
177  {
178  data(row, col) = (row == col) ? 1.0 : 0.0;
179  }
180  }
181 
182  return *this;
183 }
184 
185 /** Creates a quadratic matrix with dimension size and sets the diagonal elements to 1.0.
186  * All other elements are set to 0.0.
187  * @param size dimension of the matrix
188  * @param data_buffer if != NULL the given float array will be used as data internal data store
189  * (the object will not perform any memory management in this case)
190  * @return the id matrix object
191  */
192 Matrix
193 Matrix::get_id(unsigned int size, float *data_buffer)
194 {
195  return get_diag(size, 1.f, data_buffer);
196 }
197 
198 /** Creates a quadratic matrix with dimension size and sets the diagonal elements to value.
199  * All other elements are set to 0.0.
200  * @param size dimension of the matrix
201  * @param value of the elements of the main diagonal
202  * @param data_buffer if != NULL the given float array will be used as data internal data store
203  * (the object will not perform any memory management in this case)
204  * @return the diag matrix object
205  */
206 Matrix
207 Matrix::get_diag(unsigned int size, float value, float *data_buffer)
208 {
209  Matrix res(size, size, data_buffer, data_buffer == NULL);
210 
211  if (data_buffer != NULL)
212  {
213  unsigned int diag_elem = 0;
214  for (unsigned int i = 0; i < size * size; ++i)
215  {
216  if (i == diag_elem)
217  {
218  diag_elem += size + 1;
219  data_buffer[i] = value;
220  }
221  else data_buffer[i] = 0.f;
222  }
223  }
224  else for (unsigned int i = 0; i < size; ++i) res.data(i, i) = value;
225 
226  return res;
227 }
228 
229 /** Transposes the matrix.
230  * @return a reference to the matrix object now containing the transposed matrix
231  */
232 Matrix &
234 {
235 #ifdef HAVE_OPENCV
236  if (m_num_cols == m_num_rows)
237  {
238  CvMat cvmat = cvMat(m_num_rows, m_num_cols, CV_32FC1, m_data);
239  cvTranspose(&cvmat, &cvmat);
240 
241  return *this;
242  }
243 #endif
244  if (m_num_cols == m_num_rows) // Perform a in-place transpose
245  {
246  for (unsigned int row = 0; row < m_num_rows - 1; ++row)
247  {
248  for (unsigned int col = row + 1; col < m_num_cols; ++col)
249  {
250  float &a = data(row, col);
251  float &b = data(col, row);
252  float t = a;
253  a = b;
254  b = t;
255  }
256  }
257  }
258  else // Could not find a in-place transpose, so we use a temporary data array
259  {
260  float *new_data = (float*) malloc(sizeof(float) * m_num_elements);
261  float *cur = new_data;
262 
263  for (unsigned int col = 0; col < m_num_cols; ++col)
264  {
265  for (unsigned int row = 0; row < m_num_rows; ++row)
266  {
267  *cur++ = data(row, col);
268  }
269  }
270 
271  unsigned int cols = m_num_cols;
272  m_num_cols = m_num_rows;
273  m_num_rows = cols;
274 
275  if (m_own_memory)
276  {
277  free(m_data);
278  m_data = new_data;
279  }
280  else
281  {
282  for (unsigned int i = 0; i < m_num_elements; ++i) m_data[i] = new_data[i];
283  free(new_data);
284  }
285  }
286 
287  return *this;
288 }
289 
290 /** Computes a matrix that is the transposed of this matrix.
291  * @return a matrix that is the transposed of this matrix
292  */
293 Matrix
295 {
296  Matrix res(m_num_cols, m_num_rows);
297  float *cur = res.get_data();
298 
299  for (unsigned int col = 0; col < m_num_cols; ++col)
300  {
301  for (unsigned int row = 0; row < m_num_rows; ++row)
302  {
303  *cur++ = data(row, col);
304  }
305  }
306  return res;
307 }
308 
309 /** Inverts the matrix.
310  * The algorithm that is implemented for computing the inverse
311  * of the matrix is the Gauss-Jordan-Algorithm. Hereby, the block-
312  * matrix (A|I) consisting of the matrix to be inverted (A) and the
313  * identity matrix (I) is transformed into (I|A^(-1)).
314  * @return a reference to the matrix object which contains now the inverted matrix
315  */
316 Matrix &
318 {
319  if (m_num_rows != m_num_cols)
320  {
321  throw fawkes::Exception("Matrix::invert(): Trying to compute inverse of non-quadratic matrix!");
322  }
323 
324 #ifdef HAVE_OPENCV
325  CvMat cvmat = cvMat(m_num_rows, m_num_cols, CV_32FC1, m_data);
326  cvInv(&cvmat, &cvmat, CV_LU);
327 #else
328  Matrix i = Matrix::get_id(m_num_rows);
329 
330  // for each column...
331  for (unsigned int col = 0; col < m_num_cols; ++col)
332  {
333  // ...multiply the row by the inverse of the element
334  // on the diagonal...
335  float factor = 1.f / data(col, col);
336  i.mult_row(col, factor);
337  this->mult_row(col, factor);
338 
339  // ...and subtract that row multiplied by the elements
340  // in the current column from all other rows.
341  for (unsigned int row = 0; row < m_num_rows; ++row)
342  {
343  if (row != col)
344  {
345  float factor2 = data(row, col);
346  i.sub_row(row, col, factor2);
347  this->sub_row(row, col, factor2);
348  }
349  }
350  }
351 
352  overlay(0, 0, i);
353 #endif
354 
355  return *this;
356 }
357 
358 /** Computes a matrix that is the inverse of this matrix.
359  * @return a matrix that is the inverse of this matrix
360  */
361 Matrix
363 {
364  Matrix res(*this);
365  res.invert();
366 
367  return res;
368 }
369 
370 /** Computes the determinant of the matrix.
371  * @return the determinant
372  */
373 float
374 Matrix::det() const
375 {
376  if (m_num_rows != m_num_cols)
377  {
378  throw fawkes::Exception("Matrix::det(): The determinant can only be calculated for quadratic matrices.");
379  }
380 
381 #ifdef HAVE_OPENCV
382  CvMat cvmat = cvMat(m_num_rows, m_num_cols, CV_32FC1, m_data);
383 
384  return (float)cvDet(&cvmat);
385 #else
386  Matrix tmp_matrix(*this);
387  float result = 1.f;
388 
389  // compute the upper triangular matrix
390  for (unsigned int col = 0; col < m_num_cols; ++col)
391  {
392  float diag_elem = tmp_matrix.data(col, col);
393  result *= diag_elem;
394 
395  // multiply n-th row by m(n,n)^{-1}
396  tmp_matrix.mult_row(col, (1.f / diag_elem));
397  for (unsigned int row = col + 1; row < m_num_rows; ++row)
398  {
399  tmp_matrix.sub_row(row, col, tmp_matrix.data(row, col));
400  }
401  }
402 
403  return result;
404 #endif
405 }
406 
407 /** Returns a submatrix of the matrix.
408  * @param row the row in the original matrix of the top-left element in the submatrix
409  * @param col the column in the original matrix of the top-left element in the submatrix
410  * @param num_rows the number of rows of the submatrix
411  * @param num_cols the number of columns of the submatrix
412  * @return the submatrix
413  */
414 Matrix
415 Matrix::get_submatrix(unsigned int row, unsigned int col,
416  unsigned int num_rows, unsigned int num_cols) const
417 {
418  if ((m_num_rows < row + num_rows) || (m_num_cols < col + num_cols))
419  {
420  throw fawkes::OutOfBoundsException("Matrix::get_submatrix(): The current matrix doesn't contain a submatrix of the requested dimension at the requested position.");
421  }
422 
423  Matrix res(num_rows, num_cols);
424  float *res_data = res.get_data();
425 
426  for (unsigned int r = 0; r < num_rows; ++r)
427  {
428  for (unsigned int c = 0; c < num_cols; ++c)
429  {
430  *res_data++ = data(row + r, col + c);
431  }
432  }
433 
434  return res;
435 }
436 
437 /** Overlays another matrix over this matrix.
438  * @param row the top-most row from which onwards the the elements are
439  * exchanged for corresponding elements in the given matrix
440  * @param col the left-most column from which onwards the the elements
441  * are exchanged for corresponding elements in the given matrix
442  * @param over the matrix to be overlaid
443  */
444 void
445 Matrix::overlay(unsigned int row, unsigned int col, const Matrix &over)
446 {
447  unsigned int max_row = std::min(m_num_rows, over.m_num_rows + row);
448  unsigned int max_col = std::min(m_num_cols, over.m_num_cols + col);
449 
450  for (unsigned int r = row; r < max_row; ++r)
451  {
452  for (unsigned int c = col; c < max_col; ++c)
453  {
454  data(r, c) = over.data(r - row, c - col);
455  }
456  }
457 }
458 
459 /** (Read-only) Access-operator.
460  * With this operator it is possible to access a specific
461  * element of the matrix. (First element is at (0, 0)
462  * @param row the row of the element
463  * @param col the column of the element
464  * @return the value of the specified element
465  */
466 /* Not True: To conform with the mathematical
467  * fashion of specifying the elements of a matrix the top
468  * left element of the matrix is accessed with (1, 1)
469  * (i.e., numeration starts with 1 and not with 0).
470  */
471 float
472 Matrix::operator()(unsigned int row, unsigned int col) const
473 {
474  if (row >= m_num_rows || col >= m_num_cols)
475  {
476  throw fawkes::OutOfBoundsException("Matrix::operator() The requested element is not within the dimension of the matrix.");
477  }
478 
479  return data(row, col);
480 }
481 
482 /** (RW) Access operator.
483  * see the read-only access operator for operational details
484  * @param row the row of the element
485  * @param col the column of the element
486  * @return a reference to the specified element
487  */
488 float &
489 Matrix::operator()(unsigned int row,
490  unsigned int col)
491 {
492  if (row >= m_num_rows || col >= m_num_cols)
493  {
494  throw fawkes::OutOfBoundsException("Matrix::operator() The requested element (%d, %d) is not within the dimension of the %dx%d matrix.");
495  }
496 
497  return data(row, col);
498 }
499 
500 /** Assignment operator.
501  * Copies the data form the rhs Matrix to the lhs Matrix.
502  * @param m the rhs Matrix
503  * @return a reference to this Matrix
504  */
505 Matrix &
507 {
508  if (m_num_elements != m.m_num_elements)
509  {
510  if (!m_own_memory)
511  {
512  throw fawkes::OutOfBoundsException("Matrix::operator=(): The rhs matrix has not the same number of elements. This isn't possible if not self managing memory.");
513  }
514 
515  m_num_elements = m.m_num_elements;
516  free(m_data);
517  m_data = (float*) malloc(sizeof(float) * m_num_elements);
518  }
519 
520  m_num_rows = m.m_num_rows;
521  m_num_cols = m.m_num_cols;
522 
523  for (unsigned int i = 0; i < m_num_elements; ++i) m_data[i] = m.m_data[i];
524 
525  return *this;
526 }
527 
528 /** Matrix multiplication operator.
529  * (Matrix)a.operator*((Matrix)b) computes a * b;
530  * i.e., the 2nd matrix is right-multiplied to the 1st matrix
531  * @param rhs the right-hand-side matrix
532  * @return the product of the two matrices (a * b)
533  */
534 Matrix
535 Matrix::operator*(const Matrix &rhs) const
536 {
537  if (m_num_cols != rhs.m_num_rows)
538  {
539  throw fawkes::Exception("Matrix::operator*(...): Dimension mismatch: a %d x %d matrix can't be multiplied "
540  "with a %d x %d matrix.\n",
541  m_num_rows, m_num_cols, rhs.num_rows(), rhs.num_cols());
542  }
543 
544  unsigned int res_rows = m_num_rows;
545  unsigned int res_cols = rhs.m_num_cols;
546 
547  Matrix res(res_rows, res_cols);
548 
549  for (unsigned int r = 0; r < res_rows; ++r)
550  {
551  for (unsigned int c = 0; c < res_cols; ++c)
552  {
553  float t = 0.0f;
554 
555  for (unsigned int i = 0; i < m_num_cols; ++i)
556  {
557  t += data(r, i) * rhs.data(i, c);
558  }
559 
560  res.data(r, c) = t;
561  }
562  }
563 
564  return res;
565 }
566 
567 /** Combined matrix-multipliation and assignement operator.
568  * @param rhs the right-hand-side Matrix
569  * @return a reference to the Matrix that contains the result of the multiplication
570  */
571 Matrix &
573 {
574  if (m_num_cols != rhs.m_num_rows)
575  {
576  throw fawkes::Exception("Matrix::operator*(...): Dimension mismatch: a %d x %d matrix can't be multiplied "
577  "with a %d x %d matrix.\n",
578  m_num_rows, m_num_cols, rhs.num_rows(), rhs.num_cols());
579  }
580 
581  unsigned int res_rows = m_num_rows;
582  unsigned int res_cols = rhs.m_num_cols;
583  unsigned int res_num_elem = res_rows * res_cols;
584 
585  if (!m_own_memory && (m_num_elements != res_num_elem))
586  {
587  throw fawkes::Exception("Matrix::operator*=(): The resulting matrix has not the same number of elements. This doesn't work if not self managing memory.");
588  }
589 
590  float *new_data = (float*) malloc(sizeof(float) * res_num_elem);
591  float *cur = new_data;
592 
593  for (unsigned int r = 0; r < res_rows; ++r)
594  {
595  for (unsigned int c = 0; c < res_cols; ++c)
596  {
597  float t = 0.0f;
598 
599  for (unsigned int i = 0; i < m_num_cols; ++i)
600  {
601  t += data(r, i) * rhs.data(i, c);
602  }
603 
604  *cur++ = t;
605  }
606  }
607 
608  if (m_own_memory)
609  {
610  free(m_data);
611  m_data = new_data;
612  }
613  else
614  {
615  for (unsigned int i = 0; i < m_num_elements; ++i) m_data[i] = new_data[i];
616  free(new_data);
617  }
618 
619  return *this;
620 }
621 
622 /** Multiply the matrix with given vector.
623  * @param v a vector
624  * @return the result of the matrix-vector multiplication
625  */
626 Vector
627 Matrix::operator*(const Vector &v) const
628 {
629  unsigned int cols = v.size();
630 
631  if (m_num_cols != cols)
632  {
633  throw fawkes::Exception("Matrix::operator*(...): Dimension mismatch: a %d x %d matrix can't be multiplied "
634  "with a vector of length %d.\n", num_rows(), num_cols(), cols);
635  }
636 
637  Vector res(m_num_rows);
638  const float *vector_data = v.data_ptr();
639 
640  for (unsigned int r = 0; r < num_rows(); ++r)
641  {
642  float row_result = 0.f;
643 
644  for (unsigned int c = 0; c < cols; ++c)
645  {
646  row_result += data(r, c) * vector_data[c];
647  }
648  res[r] = row_result;
649  }
650 
651  return res;
652 }
653 
654 /** Mulitply every element of the matrix with the given scalar.
655  * @param f a scalar
656  * @return the result of the multiplication
657  */
658 Matrix
659 Matrix::operator*(const float &f) const
660 {
661  Matrix res(*this);
662  float *data = res.get_data();
663 
664  for (unsigned int i = 0; i < res.m_num_elements; ++i)
665  {
666  data[i] *= f;
667  }
668 
669  return res;
670 }
671 
672 /** Combined scalar multiplication and assignment operator.
673  * @param f a scalar
674  * @return reference to the result
675  */
676 Matrix &
677 Matrix::operator*=(const float &f)
678 {
679  for (unsigned int i = 0; i < m_num_elements; ++i)
680  {
681  m_data[i] *= f;
682  }
683 
684  return *this;
685 }
686 
687 /** Divide every element of the matrix with the given scalar.
688  * @param f a scalar
689  * @return the result of the divison
690  */
691 Matrix
692 Matrix::operator/(const float &f) const
693 {
694  Matrix res(*this);
695  float *data = res.get_data();
696 
697  for (unsigned int i = 0; i < res.m_num_elements; ++i)
698  {
699  data[i] /= f;
700  }
701 
702  return res;
703 }
704 
705 /** Combined scalar division and assignment operator.
706  * @param f a scalar
707  * @return reference to the result
708  */
709 Matrix &
710 Matrix::operator/=(const float &f)
711 {
712  for (unsigned int i = 0; i < m_num_elements; ++i)
713  {
714  m_data[i] /= f;
715  }
716 
717  return *this;
718 }
719 
720 /** Addition operator.
721  * Adds the corresponding elements of the two matrices.
722  * @param rhs the right-hand-side matrix
723  * @return the resulting matrix
724  */
725 Matrix
726 Matrix::operator+(const Matrix &rhs) const
727 {
728  if ((m_num_rows != rhs.m_num_rows) || (m_num_cols != rhs.m_num_cols))
729  {
730  throw fawkes::Exception("Matrix::operator+(...): Dimension mismatch: a %d x %d matrix can't be added to a %d x %d matrix\n",
731  num_rows(), num_cols(), rhs.num_rows(), rhs.num_cols());
732  }
733 
734  Matrix res(*this);
735  const float *rhs_d = rhs.get_data();
736  float *res_d = res.get_data();
737 
738  for (unsigned int i = 0; i < m_num_elements; ++i)
739  {
740  res_d[i] += rhs_d[i];
741  }
742 
743  return res;
744 }
745 
746 /**Add-assign operator.
747  * @param rhs the right-hand-side matrix
748  * @return a reference to the resulting matrix (this)
749  */
750 Matrix &
752 {
753  if ((m_num_rows != rhs.m_num_rows) || (m_num_cols != rhs.m_num_cols))
754  {
755  throw fawkes::Exception("Matrix::operator+(...): Dimension mismatch: a %d x %d matrix can't be added to a %d x %d matrix\n",
756  num_rows(), num_cols(), rhs.num_rows(), rhs.num_cols());
757  }
758 
759  const float *rhs_d = rhs.get_data();
760 
761  for (unsigned int i = 0; i < m_num_elements; ++i)
762  {
763  m_data[i] += rhs_d[i];
764  }
765 
766  return *this;
767 }
768 
769 /** Subtraction operator.
770  * Subtracts the corresponding elements of the two matrices.
771  * @param rhs the right-hand-side matrix
772  * @return the resulting matrix
773  */
774 Matrix
775 Matrix::operator-(const Matrix &rhs) const
776 {
777  if ((num_rows() != rhs.num_rows()) || (num_cols() != rhs.num_cols()))
778  {
779  throw fawkes::Exception("Matrix::operator-(...): Dimension mismatch: a %d x %d matrix can't be subtracted from a %d x %d matrix\n",
780  num_rows(), num_cols(), rhs.num_rows(), rhs.num_cols());
781  }
782 
783  Matrix res(*this);
784 
785  const float *rhs_d = rhs.get_data();
786  float *res_d = res.get_data();
787 
788  for (unsigned int i = 0; i < m_num_elements; ++i)
789  {
790  res_d[i] -= rhs_d[i];
791  }
792 
793  return res;
794 }
795 
796 /**Subtract-assign operator.
797  * @param rhs the right-hand-side matrix
798  * @return a reference to the resulting matrix (this)
799  */
800 Matrix &
802 {
803  if ((num_rows() != rhs.num_rows()) || (num_cols() != rhs.num_cols()))
804  {
805  throw fawkes::Exception("Matrix::operator-=(...): Dimension mismatch: a %d x %d matrix can't be subtracted from a %d x %d matrix\n",
806  num_rows(), num_cols(), rhs.num_rows(), rhs.num_cols());
807  }
808 
809  const float *rhs_d = rhs.get_data();
810 
811  for (unsigned int i = 0; i < m_num_elements; ++i)
812  {
813  m_data[i] -= rhs_d[i];
814  }
815 
816  return *this;
817 }
818 
819 /** Comparison operator.
820  * @param rhs the right-hand-side Matrix
821  * @return true if every element of this matrix is equal to the
822  * corresponding element of the other matrix
823  */
824 bool
825 Matrix::operator==(const Matrix &rhs) const
826 {
827  if ((num_rows() != rhs.num_rows()) || (num_cols() != rhs.num_cols()))
828  {
829  return false;
830  }
831 
832  const float *rhs_d = rhs.get_data();
833 
834  for (unsigned int i = 0; i < m_num_elements; ++i)
835  {
836  if (m_data[i] != rhs_d[i]) return false;
837  }
838 
839  return true;
840 }
841 
842 /** Changes the matrix by multiplying a row with a factor.
843  * @param row the row
844  * @param factor the factor
845  */
846 void
847 Matrix::mult_row(unsigned int row, float factor)
848 {
849  if (row >= m_num_rows)
850  {
851  throw fawkes::OutOfBoundsException("Matrix::mult_row(...)", row, 0, m_num_rows);
852  }
853 
854  for (unsigned int col = 0; col < m_num_cols; ++col)
855  {
856  data(row, col) *= factor;
857  }
858 }
859 
860 /** For two rows A and B and a factor f, A is changed to A - f*B.
861  * @param row_a the row that is changed
862  * @param row_b the row that is substracted from row_a
863  * @param factor the factor by which every element of row_b is multiplied before it is
864  * substracted from row_a
865  */
866 void
867 Matrix::sub_row(unsigned int row_a, unsigned int row_b, float factor)
868 {
869  if (row_a >= m_num_rows)
870  {
871  throw fawkes::OutOfBoundsException("Matrix::sub_row(...) row_a", row_a, 0, m_num_rows);
872  }
873  if (row_b >= m_num_rows)
874  {
875  throw fawkes::OutOfBoundsException("Matrix::sub_row(...) row_b", row_b, 0, m_num_rows);
876  }
877 
878  for (unsigned int col = 0; col < m_num_cols; ++col)
879  {
880  data(row_a, col) -= factor * data(row_b, col);
881  }
882 }
883 
884 /** Print matrix to standard out.
885  * @param name a name that is printed before the content of the matrix (not required)
886  * @param col_sep a string used to separate columns (defaults to '\\t')
887  * @param row_sep a string used to separate rows (defaults to '\\n')
888  */
889 void
890 Matrix::print_info(const char *name, const char *col_sep, const char *row_sep) const
891 {
892  if (name)
893  {
894  printf("%s:\n", name);
895  }
896 
897  for (unsigned int r = 0; r < num_rows(); ++r)
898  {
899  printf((r == 0 ? "[" : " "));
900  for (unsigned int c = 0; c < num_cols(); ++c)
901  {
902  printf("%f", (*this)(r, c));
903  if (c+1 < num_cols())
904  {
905  if (col_sep) printf("%s", col_sep);
906  else printf("\t");
907  }
908  }
909  if (r+1 < num_rows())
910  {
911  if (row_sep) printf("%s", row_sep);
912  else printf("\n");
913  }
914  else printf("]\n\n");
915  }
916 }
917 
918 } // end namespace fawkes
float * data_ptr()
Get pointer to the internal data container.
Definition: vector.h:41
Matrix & id()
Sets the diagonal elements to 1.0 and all other to 0.0.
Definition: matrix.cpp:172
unsigned int num_cols() const
Return the number of columns in the Matrix.
Definition: matrix.h:43
Matrix operator*(const Matrix &rhs) const
Matrix multiplication operator.
Definition: matrix.cpp:535
Matrix & operator-=(const Matrix &rhs)
Subtract-assign operator.
Definition: matrix.cpp:801
Matrix & operator*=(const Matrix &rhs)
Combined matrix-multipliation and assignement operator.
Definition: matrix.cpp:572
Fawkes library namespace.
A simple column vector.
Definition: vector.h:31
A general matrix class.
Definition: matrix.h:33
float det() const
Computes the determinant of the matrix.
Definition: matrix.cpp:374
Matrix operator/(const float &f) const
Divide every element of the matrix with the given scalar.
Definition: matrix.cpp:692
Matrix & operator/=(const float &f)
Combined scalar division and assignment operator.
Definition: matrix.cpp:710
unsigned int num_rows() const
Return the number of rows in the Matrix.
Definition: matrix.h:42
Matrix operator+(const Matrix &rhs) const
Addition operator.
Definition: matrix.cpp:726
void print_info(const char *name=0, const char *col_sep=0, const char *row_sep=0) const
Print matrix to standard out.
Definition: matrix.cpp:890
Matrix & operator+=(const Matrix &rhs)
Add-assign operator.
Definition: matrix.cpp:751
unsigned int size() const
Get the number of elements.
Definition: vector.cpp:111
void overlay(unsigned int row, unsigned int col, const Matrix &m)
Overlays another matrix over this matrix.
Definition: matrix.cpp:445
const float * get_data() const
Returns the const data pointer.
Definition: matrix.h:57
void size(unsigned int &num_rows, unsigned int &num_cols) const
Determines the dimensions of the matrix.
Definition: matrix.cpp:161
static Matrix get_id(unsigned int size, float *data_buffer=0)
Creates a quadratic matrix with dimension size and sets the diagonal elements to 1.0.
Definition: matrix.cpp:193
bool operator==(const Matrix &rhs) const
Comparison operator.
Definition: matrix.cpp:825
Matrix & transpose()
Transposes the matrix.
Definition: matrix.cpp:233
Matrix get_inverse() const
Computes a matrix that is the inverse of this matrix.
Definition: matrix.cpp:362
Base class for exceptions in Fawkes.
Definition: exception.h:36
Matrix & invert()
Inverts the matrix.
Definition: matrix.cpp:317
Matrix get_submatrix(unsigned int row, unsigned int col, unsigned int num_rows, unsigned int num_cols) const
Returns a submatrix of the matrix.
Definition: matrix.cpp:415
float operator()(unsigned int row, unsigned int col) const
(Read-only) Access-operator.
Definition: matrix.cpp:472
Matrix operator-(const Matrix &rhs) const
Subtraction operator.
Definition: matrix.cpp:775
Matrix get_transpose() const
Computes a matrix that is the transposed of this matrix.
Definition: matrix.cpp:294
Matrix & operator=(const Matrix &rhs)
Assignment operator.
Definition: matrix.cpp:506
Index out of bounds.
Definition: software.h:88
~Matrix()
Destructor.
Definition: matrix.cpp:151
Matrix(unsigned int num_rows=0, unsigned int num_cols=0, float *data=0, bool manage_own_memory=true)
Constructor.
Definition: matrix.cpp:102
static Matrix get_diag(unsigned int size, float value, float *data_buffer=0)
Creates a quadratic matrix with dimension size and sets the diagonal elements to value.
Definition: matrix.cpp:207