GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GMatrix.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GMatrix.cpp - General matrix class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2006-2017 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GMatrix.cpp
23  * @brief General matrix class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GVector.hpp"
36 #include "GMatrix.hpp"
37 #include "GMatrixSparse.hpp"
38 #include "GMatrixSymmetric.hpp"
39 
40 /* __ Method name definitions ____________________________________________ */
41 #define G_CONSTRUCTOR "GMatrix::GMatrix(int&, int&)"
42 #define G_OP_MUL_VEC "GMatrix::operator*(GVector&)"
43 #define G_OP_ADD "GMatrix::operator+=(GMatrix&)"
44 #define G_OP_SUB "GMatrix::operator-=(GMatrix&)"
45 #define G_OP_MAT_MUL "GMatrix::operator*=(GMatrix&)"
46 #define G_AT "GMatrix::at(int&, int&)"
47 #define G_EXTRACT_ROW "GMatrix::row(int&)"
48 #define G_SET_ROW "GMatrix::row(int&, GVector&)"
49 #define G_EXTRACT_COLUMN "GMatrix::column(int&)"
50 #define G_SET_COLUMN "GMatrix::column(int&, GVector&)"
51 #define G_ADD_TO_ROW "GMatrix::add_to_row(int&, GVector&)"
52 #define G_ADD_TO_COLUMN "GMatrix::add_to_column(int&, GVector&)"
53 #define G_INVERT "GMatrix::invert()"
54 #define G_SOLVE "GMatrix::solve(GVector&)"
55 #define G_EXTRACT_LOWER "GMatrix::extract_lower_triangle()"
56 #define G_EXTRACT_UPPER "GMatrix::extract_upper_triangle()"
57 
58 
59 /*==========================================================================
60  = =
61  = Constructors/destructors =
62  = =
63  ==========================================================================*/
64 
65 /***********************************************************************//**
66  * @brief Void matrix constructor
67  *
68  * Constructs an empty matrix (i.e. a matrix without any elements; the number
69  * of rows and columns of the matrix will be zero).
70  *
71  * @todo Verify that the class is save against empty matrix objects.
72  ***************************************************************************/
74 {
75  // Initialise class members
76  init_members();
77 
78  // Return
79  return;
80 }
81 
82 
83 /***********************************************************************//**
84  * @brief Matrix constructor
85  *
86  * @param[in] rows Number of rows [>0].
87  * @param[in] columns Number of columns [>0].
88  *
89  * @exception GException::empty
90  * Specified number of rows or columns is not valid.
91  *
92  * Constructs a matrix with the specified number of @p rows and @p columns.
93  * The memory for all matrix elements will be allocated and all matrix
94  * elements will be initialized to 0.
95  ***************************************************************************/
96 GMatrix::GMatrix(const int& rows, const int& columns) : GMatrixBase()
97 {
98  // Continue only if matrix is valid
99  if (rows > 0 && columns > 0) {
100 
101  // Initialise class members for clean destruction
102  init_members();
103 
104  // Allocate matrix memory
105  alloc_members(rows, columns);
106 
107  }
108  else {
110  }
111 
112  // Return
113  return;
114 }
115 
116 
117 /***********************************************************************//**
118  * @brief Copy constructor
119  *
120  * @param[in] matrix Matrix.
121  ***************************************************************************/
122 GMatrix::GMatrix(const GMatrix& matrix) : GMatrixBase(matrix)
123 {
124  // Initialise class members for clean destruction
125  init_members();
126 
127  // Copy members
128  copy_members(matrix);
129 
130  // Return
131  return;
132 }
133 
134 
135 /***********************************************************************//**
136  * @brief GMatrixSymmetric to GMatrix storage class convertor
137  *
138  * @param[in] matrix Symmetric matrix (GMatrixSymmetric).
139  *
140  * Constructs a GMatrix object from a symmetric matrix object of type
141  * GMatrixSymmetric. Since the result is a generic matrix, the constructor
142  * will succeed in all cases.
143  ***************************************************************************/
145 {
146  // Initialise class members for clean destruction
147  init_members();
148 
149  // Allocate matrix memory
150  alloc_members(matrix.rows(), matrix.columns());
151 
152  // Fill matrix. We benefit here from the symmetry of the matrix and
153  // loop only over the lower triangle of the symmetric matrix to perform
154  // the fill.
155  for (int col = 0; col < m_cols; ++col) {
156  for (int row = col; row < m_rows; ++row) {
157  double value = matrix(row, col);
158  (*this)(row, col) = value;
159  (*this)(col, row) = value;
160  }
161  }
162 
163  // Return
164  return;
165 }
166 
167 
168 /***********************************************************************//**
169  * @brief GMatrixSparse to GMatrix storage class convertor
170  *
171  * @param[in] matrix Sparse matrix (GMatrixSparse).
172  *
173  * Constructs a GMatrix object from a sparse matrix object of type
174  * GMatrixSparse. Since the result is a generic matrix, the constructor
175  * will succeed in all cases.
176  ***************************************************************************/
178 {
179  // Initialise class members for clean destruction
180  init_members();
181 
182  // Construct matrix
183  alloc_members(matrix.rows(), matrix.columns());
184 
185  // Fill matrix
186  for (int col = 0; col < m_cols; ++col) {
187  for (int row = 0; row < m_rows; ++row) {
188  (*this)(row, col) = matrix(row, col);
189  }
190  }
191 
192  // Return
193  return;
194 }
195 
196 
197 /***********************************************************************//**
198  * @brief Destructor
199  ***************************************************************************/
201 {
202  // Free members
203  free_members();
204 
205  // Return
206  return;
207 }
208 
209 
210 /*==========================================================================
211  = =
212  = Operators =
213  = =
214  ==========================================================================*/
215 
216 /***********************************************************************//**
217  * @brief Matrix assignment operator
218  *
219  * @param[in] matrix Matrix.
220  * @return Matrix.
221  *
222  * Assigns the content of another matrix to the actual matrix instance.
223  ***************************************************************************/
225 {
226  // Execute only if object is not identical
227  if (this != &matrix) {
228 
229  // Assign base class members. Note that this method will also perform
230  // the allocation of the matrix memory and the copy of the matrix
231  // attributes.
232  this->GMatrixBase::operator=(matrix);
233 
234  // Free derived class members
235  free_members();
236 
237  // Initialise derived class members
238  init_members();
239 
240  // Copy derived class members
241  copy_members(matrix);
242 
243  } // endif: object was not identical
244 
245  // Return this object
246  return *this;
247 }
248 
249 
250 /***********************************************************************//**
251  * @brief Value assignment operator
252  *
253  * @param[in] value Value.
254  * @return Matrix.
255  *
256  * Assigns the specified @p value to all elements of the matrix.
257  ***************************************************************************/
258 GMatrix& GMatrix::operator=(const double& value)
259 {
260  // Assign value
261  double* ptr = m_data;
262  for (int i = 0; i < m_elements; ++i) {
263  *ptr++ = value;
264  }
265 
266  // Return this object
267  return *this;
268 }
269 
270 
271 /***********************************************************************//**
272  * @brief Vector multiplication
273  *
274  * @param[in] vector Vector.
275  *
276  * @exception GException::matrix_vector_mismatch
277  * Vector length differs from number of columns in matrix.
278  *
279  * This method performs a vector multiplication of a matrix. The vector
280  * multiplication will produce a vector. The matrix multiplication can only
281  * be performed when the number of matrix columns is identical to the length
282  * of the vector.
283  ***************************************************************************/
284 GVector GMatrix::operator*(const GVector& vector) const
285 {
286  // Raise an exception if the matrix and vector dimensions are not compatible
287  if (m_cols != vector.size()) {
289  m_rows, m_cols);
290  }
291 
292  // Perform vector multiplication
293  GVector result(m_rows);
294  for (int row = 0; row < m_rows; ++row) {
295  double sum = 0.0;
296  const double* src = m_data + row;
297  for (int col = 0; col < m_cols; ++col, src += m_rows) {
298  sum += *src * vector[col];
299  }
300  result[row] = sum;
301  }
302 
303  // Return result
304  return result;
305 }
306 
307 
308 /***********************************************************************//**
309  * @brief Negate matrix elements
310  *
311  * @return Matrix with negated elements.
312  *
313  * Returns a matrix where each element has been replaced by its negative
314  * element.
315  ***************************************************************************/
317 {
318  // Copy matrix
319  GMatrix matrix = *this;
320 
321  // Take the absolute value of all matrix elements
322  double* ptr = matrix.m_data;
323  for (int i = 0; i < matrix.m_elements; ++i, ++ptr) {
324  *ptr = -(*ptr);
325  }
326 
327  // Return matrix
328  return matrix;
329 }
330 
331 
332 /***********************************************************************//**
333  * @brief Unary matrix addition operator
334  *
335  * @param[in] matrix Matrix.
336  *
337  * @exception GException::matrix_mismatch
338  * Incompatible matrix size.
339  *
340  * This method performs a matrix addition. The operation can only succeed
341  * when the dimensions of both matrices are identical.
342  ***************************************************************************/
344 {
345  // Raise an exception if the matrix dimensions are not compatible
346  if (m_rows != matrix.m_rows || m_cols != matrix.m_cols) {
348  m_rows, m_cols,
349  matrix.m_rows, matrix.m_cols);
350  }
351 
352  // Add all matrix elements
353  const double* src = matrix.m_data;
354  double* dst = m_data;
355  for (int i = 0; i < m_elements; ++i) {
356  *dst++ += *src++;
357  }
358 
359  // Return result
360  return *this;
361 }
362 
363 
364 /***********************************************************************//**
365  * @brief Unary matrix scalar addition operator
366  *
367  * @param[in] scalar Scalar.
368  *
369  * Adds a @p scalar to each matrix element.
370  ***************************************************************************/
371 GMatrix& GMatrix::operator+=(const double& scalar)
372 {
373  // Add all matrix elements
374  double* dst = m_data;
375  for (int i = 0; i < m_elements; ++i) {
376  *dst++ += scalar;
377  }
378 
379  // Return result
380  return *this;
381 }
382 
383 
384 /***********************************************************************//**
385  * @brief Unary matrix subtraction operator
386  *
387  * @param[in] matrix Matrix.
388  *
389  * @exception GException::matrix_mismatch
390  * Incompatible matrix size.
391  *
392  * This method performs a matrix addition. The operation can only succeed
393  * when the dimensions of both matrices are identical.
394  ***************************************************************************/
396 {
397  // Raise an exception if the matrix dimensions are not compatible
398  if (m_rows != matrix.m_rows || m_cols != matrix.m_cols) {
400  m_rows, m_cols,
401  matrix.m_rows, matrix.m_cols);
402  }
403 
404  // Subtract all matrix elements
405  const double* src = matrix.m_data;
406  double* dst = m_data;
407  for (int i = 0; i < m_elements; ++i) {
408  *dst++ -= *src++;
409  }
410 
411  // Return result
412  return *this;
413 }
414 
415 
416 /***********************************************************************//**
417  * @brief Unary matrix scalar subtraction operator
418  *
419  * @param[in] scalar Scalar.
420  *
421  * Subtracts a @p scalar from each matrix element.
422  ***************************************************************************/
423 GMatrix& GMatrix::operator-=(const double& scalar)
424 {
425  // Add all matrix elements
426  double* dst = m_data;
427  for (int i = 0; i < m_elements; ++i) {
428  *dst++ -= scalar;
429  }
430 
431  // Return result
432  return *this;
433 }
434 
435 
436 /***********************************************************************//**
437  * @brief Unary matrix multiplication operator
438  *
439  * @param[in] matrix Matrix.
440  *
441  * @exception GException::matrix_mismatch
442  * Incompatible matrix size.
443  *
444  * This method performs a matrix multiplication. The operation can only
445  * succeed when the dimensions of both matrices are compatible.
446  *
447  * In case of rectangular matrices the result matrix does not change and
448  * the operation is performed inplace. For the general case the result
449  * matrix changes in size and for simplicity a new matrix is allocated to
450  * hold the result.
451  ***************************************************************************/
453 {
454  // Raise an exception if the matrix dimensions are not compatible
455  if (m_cols != matrix.m_rows) {
457  m_rows, m_cols,
458  matrix.m_rows, matrix.m_cols);
459  }
460 
461  // Case A: Matrices are rectangular, so perform 'inplace' multiplication
462  if (m_rows == m_cols) {
463  for (int row = 0; row < m_rows; ++row) {
464  GVector v_row = this->row(row);
465  for (int col = 0; col < m_cols; ++col) {
466  double sum = 0.0;
467  for (int i = 0; i < m_cols; ++i) {
468  sum += v_row[i] * matrix(i,col);
469  }
470  (*this)(row,col) = sum;
471  }
472  }
473  }
474 
475  // Case B: Matrices are not rectangular, so we cannot work inplace
476  else {
477 
478  // Allocate result matrix
479  GMatrix result(m_rows, matrix.m_cols);
480 
481  // Loop over all elements of result matrix
482  for (int row = 0; row < m_rows; ++row) {
483  for (int col = 0; col < matrix.m_cols; ++col) {
484  double sum = 0.0;
485  for (int i = 0; i < m_cols; ++i) {
486  sum += (*this)(row,i) * matrix(i,col);
487  }
488  result(row,col) = sum;
489  }
490  }
491 
492  // Assign result
493  *this = result;
494 
495  } // endelse: matrices were not rectangular
496 
497  // Return result
498  return *this;
499 }
500 
501 
502 /*==========================================================================
503  = =
504  = Public methods =
505  = =
506  ==========================================================================*/
507 
508 /***********************************************************************//**
509  * @brief Clear matrix
510  ***************************************************************************/
511 void GMatrix::clear(void)
512 {
513  // Free members
514  free_members();
515 
516  // Initialise private members
517  init_members();
518 
519  // Return
520  return;
521 }
522 
523 
524 /***********************************************************************//**
525  * @brief Clone matrix
526  *
527  * @return Pointer to deep copy of matrix.
528  ***************************************************************************/
530 {
531  // Clone matrix
532  return new GMatrix(*this);
533 }
534 
535 
536 /***********************************************************************//**
537  * @brief Return reference to matrix element
538  *
539  * @param[in] row Matrix row [0,...,rows()-1].
540  * @param[in] column Matrix column [0,...,columns()-1].
541  * @return Reference to matrix element.
542  *
543  * @exception GException::out_of_range
544  * Row or column index out of range.
545  *
546  * Returns a reference to the matrix element at @p row and @p column.
547  * Verifies the validity of the @p row and @p column argument.
548  ***************************************************************************/
549 double& GMatrix::at(const int& row, const int& column)
550 {
551  // Raise exception if row or column index is out of range
552  if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) {
553  throw GException::out_of_range(G_AT, row, column, m_rows, m_cols);
554  }
555 
556  // Return element
557  return (m_data[m_colstart[column]+row]);
558 }
559 
560 
561 /***********************************************************************//**
562  * @brief Return reference to matrix element (const version)
563  *
564  * @param[in] row Matrix row [0,...,rows()-1].
565  * @param[in] column Matrix column [0,...,columns()-1].
566  * @return Const reference to matrix element.
567  *
568  * @exception GException::out_of_range
569  * Row or column index out of range.
570  *
571  * Returns a const reference to the matrix element at @p row and @p column.
572  * Verifies the validity of the @p row and @p column argument.
573  ***************************************************************************/
574 const double& GMatrix::at(const int& row, const int& column) const
575 {
576  // Raise exception if row or column index is out of range
577  if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) {
578  throw GException::out_of_range(G_AT, row, column, m_rows, m_cols);
579  }
580 
581  // Return element
582  return (m_data[m_colstart[column]+row]);
583 }
584 
585 
586 
587 /***********************************************************************//**
588  * @brief Extract row as vector from matrix
589  *
590  * @param[in] row Matrix row [0,...,rows()-1].
591  * @return Vector of matrix row elements (columns() elements).
592  *
593  * @exception GException::out_of_range
594  * Invalid row index specified.
595  *
596  * Extracts one @p row of the matrix into a vector. The vector will contain
597  * columns() elements.
598  ***************************************************************************/
599 GVector GMatrix::row(const int& row) const
600 {
601  // Raise an exception if the row index is invalid
602  #if defined(G_RANGE_CHECK)
603  if (row < 0 || row >= m_rows) {
605  }
606  #endif
607 
608  // Create result vector
609  GVector result(m_cols);
610 
611  // Extract row into vector. Recall that the matrix elements are stored
612  // column wise, hence we have to increment the element counter i by
613  // the number of rows to go into the next matrix column.
614  for (int col = 0, i = row; col < m_cols; ++col, i += m_rows) {
615  result[col] = m_data[i];
616  }
617 
618  // Return vector
619  return result;
620 }
621 
622 
623 /***********************************************************************//**
624  * @brief Set row in matrix
625  *
626  * @param[in] row Matrix row [0,...,rows()-1].
627  * @param[in] vector Vector of matrix row elements (columns() elements).
628  *
629  * @exception GException::out_of_range
630  * Invalid row index specified.
631  * @exception GException::matrix_vector_mismatch
632  * Vector does not match the matrix dimensions.
633  *
634  * Sets the elements from a vector as the elements of a matrix row. The
635  * length of the vector must be identical to the number of columns in the
636  * matrix.
637  ***************************************************************************/
638 void GMatrix::row(const int& row, const GVector& vector)
639 {
640  // Raise an exception if the row index is invalid
641  #if defined(G_RANGE_CHECK)
642  if (row < 0 || row >= m_rows) {
643  throw GException::out_of_range(G_SET_ROW, row, 0, m_rows-1);
644  }
645  #endif
646 
647  // Raise an exception if the matrix and vector dimensions are not
648  // compatible
649  if (m_cols != vector.size()) {
651  m_rows, m_cols);
652  }
653 
654  // Set row from vector. Recall that the matrix elements are stored
655  // column wise, hence we have to increment the element counter i by
656  // the number of rows to go into the next matrix column.
657  for (int col = 0, i = row; col < m_cols; ++col, i += m_rows) {
658  m_data[i] = vector[col];
659  }
660 
661  // Return
662  return;
663 }
664 
665 
666 /***********************************************************************//**
667  * @brief Extract column as vector from matrix
668  *
669  * @param[in] column Column index [0,...,columns()-1].
670  *
671  * @exception GException::out_of_range
672  * Invalid row index specified.
673  *
674  * This method extracts a matrix column into a vector.
675  ***************************************************************************/
676 GVector GMatrix::column(const int& column) const
677 {
678  // Raise an exception if the column index is invalid
679  #if defined(G_RANGE_CHECK)
680  if (column < 0 || column >= m_cols) {
681  throw GException::out_of_range(G_EXTRACT_COLUMN, column, 0, m_cols-1);
682  }
683  #endif
684 
685  // Create result vector
686  GVector result(m_rows);
687 
688  // Continue only if we have elements
689  if (m_elements > 0) {
690 
691  // Get start of data in matrix
692  int i = m_colstart[column];
693 
694  // Extract column into vector
695  for (int row = 0; row < m_rows; ++row) {
696  result[row] = m_data[i++];
697  }
698 
699  } // endif: matrix had elements
700 
701  // Return vector
702  return result;
703 }
704 
705 
706 /***********************************************************************//**
707  * @brief Set matrix column from vector
708  *
709  * @param[in] column Column index [0,...,columns()-1].
710  * @param[in] vector Vector.
711  *
712  * @exception GException::out_of_range
713  * Invalid column index specified.
714  * @exception GException::matrix_vector_mismatch
715  * Matrix dimension mismatches the vector size.
716  *
717  * Inserts the content of a vector into a matrix column. Any previous
718  * content in the matrix column will be overwritted.
719  ***************************************************************************/
720 void GMatrix::column(const int& column, const GVector& vector)
721 {
722  // Raise an exception if the column index is invalid
723  #if defined(G_RANGE_CHECK)
724  if (column < 0 || column >= m_cols) {
725  throw GException::out_of_range(G_SET_COLUMN, column, 0, m_cols-1);
726  }
727  #endif
728 
729  // Raise an exception if the matrix and vector dimensions are not
730  // compatible
731  if (m_rows != vector.size()) {
733  m_rows, m_cols);
734  }
735 
736  // Continue only if we have elements
737  if (m_elements > 0) {
738 
739  // Get start index of data in matrix
740  int i = m_colstart[column];
741 
742  // Insert column into vector
743  for (int row = 0; row < m_rows; ++row) {
744  m_data[i++] = vector[row];
745  }
746 
747  } // endif: matrix had elements
748 
749  // Return
750  return;
751 }
752 
753 
754 /***********************************************************************//**
755  * @brief Add row to matrix elements
756  *
757  * @param[in] row Matrix row [0,...,rows()-1].
758  * @param[in] vector Vector of matrix row elements (columns() elements).
759  *
760  * @exception GException::out_of_range
761  * Invalid row index specified.
762  * @exception GException::matrix_vector_mismatch
763  * Vector does not match the matrix dimensions.
764  ***************************************************************************/
765 void GMatrix::add_to_row(const int& row, const GVector& vector)
766 {
767  // Raise an exception if the row index is invalid
768  #if defined(G_RANGE_CHECK)
769  if (row < 0 || row >= m_rows) {
770  throw GException::out_of_range(G_ADD_TO_ROW, row, 0, m_rows-1);
771  }
772  #endif
773 
774  // Raise an exception if the matrix and vector dimensions are not
775  // compatible
776  if (m_cols != vector.size()) {
778  m_rows, m_cols);
779  }
780 
781  // Add row from vector. Recall that the matrix elements are stored
782  // column wise, hence we have to increment the element counter i by
783  // the number of rows to go into the next matrix column.
784  for (int col = 0, i = row; col < m_cols; ++col, i += m_rows) {
785  m_data[i] += vector[col];
786  }
787 
788  // Return
789  return;
790 }
791 
792 
793 /***********************************************************************//**
794  * @brief Add vector column into matrix
795  *
796  * @param[in] column Column index [0,...,columns()-1].
797  * @param[in] vector Vector.
798  *
799  * @exception GException::out_of_range
800  * Invalid column index specified.
801  * @exception GException::matrix_vector_mismatch
802  * Matrix dimension mismatches the vector size.
803  *
804  * Adds the content of a vector to a matrix column.
805  ***************************************************************************/
806 void GMatrix::add_to_column(const int& column, const GVector& vector)
807 {
808  // Raise an exception if the column index is invalid
809  #if defined(G_RANGE_CHECK)
810  if (column < 0 || column >= m_cols) {
811  throw GException::out_of_range(G_ADD_TO_COLUMN, column, 0, m_cols-1);
812  }
813  #endif
814 
815  // Raise an exception if the matrix and vector dimensions are not
816  // compatible
817  if (m_rows != vector.size()) {
819  m_rows, m_cols);
820  }
821 
822  // Continue only if we have elements
823  if (m_elements > 0) {
824 
825  // Get start index of data in matrix
826  int i = m_colstart[column];
827 
828  // Add column into vector
829  for (int row = 0; row < m_rows; ++row) {
830  m_data[i++] += vector[row];
831  }
832 
833  } // endif: matrix had elements
834 
835  // Return
836  return;
837 }
838 
839 
840 /***********************************************************************//**
841  * @brief Return transposed matrix
842  *
843  * @return Transposed matrix.
844  *
845  * Returns transposed matrix of the matrix.
846  ***************************************************************************/
848 {
849  // Allocate result matrix
850  GMatrix matrix(m_cols, m_rows);
851 
852  // Transpose matrix
853  for (int row = 0; row < m_rows; ++row) {
854  for (int col = 0; col < m_cols; ++col) {
855  matrix(col, row) = (*this)(row, col);
856  }
857  }
858 
859  // Return matrix
860  return matrix;
861 }
862 
863 
864 /***********************************************************************//**
865  * @brief Return inverted matrix
866  *
867  * @return Inverted matrix.
868  *
869  * @exception GException::feature_not_implemented
870  * Feature not yet implemented.
871  *
872  * Returns inverse of matrix.
873  *
874  * @todo Needs to be implemented.
875  ***************************************************************************/
877 {
878  // Allocate result matrix
879  GMatrix matrix(m_cols, m_rows);
880 
881  // Throw exception
883 
884  // Return matrix
885  return matrix;
886 }
887 
888 
889 /***********************************************************************//**
890  * @brief Solves linear matrix equation
891  *
892  * @param[in] vector Solution vector.
893  *
894  * @exception GException::feature_not_implemented
895  * Feature not yet implemented.
896  *
897  * Solves the linear equation
898  *
899  * \f[M \times {\tt solution} = {\tt vector} \f]
900  *
901  * where \f$M\f$ is the matrix, \f${\tt vector}\f$ is the result, and
902  * \f${\tt solution}\f$ is the solution.
903  *
904  * @todo Needs to be implemented.
905  ***************************************************************************/
906 GVector GMatrix::solve(const GVector& vector) const
907 {
908  // Allocate result vector
909  GVector result;
910 
911  // Throw exception
913 
914  // Return vector
915  return result;
916 }
917 
918 
919 /***********************************************************************//**
920  * @brief Return absolute of matrix
921  *
922  * @return Absolute of matrix
923  *
924  * Returns matrix where all elements of the matrix have been replaced by
925  * their absolute values.
926  ***************************************************************************/
928 {
929  // Allocate result matrix
930  GMatrix matrix(m_rows, m_cols);
931 
932  // Take the absolute value of all matrix elements
933  double* src = m_data;
934  double* dst = matrix.m_data;
935  for (int i = 0; i < m_elements; ++i) {
936  *dst++ = std::abs(*src++);
937  }
938 
939  // Return matrix
940  return matrix;
941 }
942 
943 
944 /***********************************************************************//**
945  * @brief Return fill of matrix
946  *
947  * @return Matrix fill (between 0 and 1).
948  *
949  * Returns the fill of the matrix. The fill of a matrix is defined as the
950  * number non-zero elements devided by the number of total elements. By
951  * definition, the fill is comprised in the interval [0,..,1].
952  *
953  * The fill of a matrix with zero elements will be set to 0.
954  ***************************************************************************/
955 double GMatrix::fill(void) const
956 {
957  // Initialise result
958  double result = 0.0;
959 
960  // Continue only if matrix has elements
961  if (m_elements > 0) {
962 
963  // Determine the number of zero elements
964  int zero = 0;
965  for (int i = 0; i < m_elements; ++i) {
966  if (m_data[i] == 0.0) {
967  zero++;
968  }
969  }
970 
971  // Compute fill
972  result = 1.0-double(zero)/double(m_elements);
973 
974  } // endif: there were elements in the matrix
975 
976  // Return result
977  return result;
978 }
979 
980 
981 /***********************************************************************//**
982  * @brief Extract lower triangle of matrix
983  *
984  * @exception GException::matrix_not_square
985  * Matrix is not a square matrix.
986  *
987  * This method extracts the lower triangle of a matrix into another matrix.
988  * (including the diagonal elements).
989  * All remaining matrix elements will be zero.
990  *
991  * Triangular extraction only works for square matrixes.
992  ***************************************************************************/
994 {
995  // Raise an exception if matrix is not square
996  if (m_rows != m_cols) {
998  }
999 
1000  // Define result matrix
1001  GMatrix result(m_rows, m_cols);
1002 
1003  // Extract all elements
1004  for (int col = 0; col < m_cols; ++col) {
1005  int i = m_colstart[col] + col;
1006  for (int row = col; row < m_rows; ++row, ++i) {
1007  result.m_data[i] = m_data[i];
1008  }
1009  }
1010 
1011  // Return result
1012  return result;
1013 }
1014 
1015 
1016 /***********************************************************************//**
1017  * @brief Extract upper triangle of matrix
1018  *
1019  * @exception GException::matrix_not_square
1020  * Matrix is not a square matrix.
1021  *
1022  * This method extracts the upper triangle of a matrix into another matrix.
1023  * (including the diagonal elements).
1024  * All remaining matrix elements will be zero.
1025  *
1026  * Triangular extraction only works for square matrixes.
1027  ***************************************************************************/
1029 {
1030  // Raise an exception if matrix is not squared
1031  if (m_rows != m_cols) {
1033  }
1034 
1035  // Define result matrix
1036  GMatrix result(m_rows, m_cols);
1037 
1038  // Extract all elements
1039  for (int col = 0; col < m_cols; ++col) {
1040  int i = m_colstart[col];
1041  for (int row = 0; row <= col; ++row, ++i) {
1042  result.m_data[i] = m_data[i];
1043  }
1044  }
1045 
1046  // Return result
1047  return result;
1048 }
1049 
1050 
1051 /***********************************************************************//**
1052  * @brief Set Euler rotation matrix around x axis
1053  *
1054  * @param[in] angle Rotation angle (degrees)
1055  ***************************************************************************/
1056 void GMatrix::eulerx(const double& angle)
1057 {
1058  // Free members
1060 
1061  // Initialise members
1063 
1064  // Construct 3*3 matrix
1065  alloc_members(3,3);
1066 
1067  // Compute angles
1068  double arg = angle * gammalib::deg2rad;
1069  double cosangle = std::cos(arg);
1070  double sinangle = std::sin(arg);
1071 
1072  // Set matrix elements
1073  (*this)(0,0) = 1.0;
1074  (*this)(0,1) = 0.0;
1075  (*this)(0,2) = 0.0;
1076  (*this)(1,0) = 0.0;
1077  (*this)(1,1) = cosangle;
1078  (*this)(1,2) = -sinangle;
1079  (*this)(2,0) = 0.0;
1080  (*this)(2,1) = sinangle;
1081  (*this)(2,2) = cosangle;
1082 
1083  // Return
1084  return;
1085 }
1086 
1087 
1088 /***********************************************************************//**
1089  * @brief Set Euler rotation matrix around y axis
1090  *
1091  * @param[in] angle Rotation angle (degrees)
1092  ***************************************************************************/
1093 void GMatrix::eulery(const double& angle)
1094 {
1095  // Free members
1097 
1098  // Initialise members
1100 
1101  // Construct 3*3 matrix
1102  alloc_members(3,3);
1103 
1104  // Compute angles
1105  double arg = angle * gammalib::deg2rad;
1106  double cosangle = std::cos(arg);
1107  double sinangle = std::sin(arg);
1108 
1109  // Set matrix elements
1110  (*this)(0,0) = cosangle;
1111  (*this)(0,1) = 0.0;
1112  (*this)(0,2) = sinangle;
1113  (*this)(1,0) = 0.0;
1114  (*this)(1,1) = 1.0;
1115  (*this)(1,2) = 0.0;
1116  (*this)(2,0) = -sinangle;
1117  (*this)(2,1) = 0.0;
1118  (*this)(2,2) = cosangle;
1119 
1120  // Return
1121  return;
1122 }
1123 
1124 
1125 /***********************************************************************//**
1126  * @brief Set Euler rotation matrix around z axis
1127  *
1128  * @param[in] angle Rotation angle (degrees)
1129  ***************************************************************************/
1130 void GMatrix::eulerz(const double& angle)
1131 {
1132  // Free members
1134 
1135  // Initialise members
1137 
1138  // Construct 3*3 matrix
1139  alloc_members(3,3);
1140 
1141  // Compute angles
1142  double arg = angle * gammalib::deg2rad;
1143  double cosangle = std::cos(arg);
1144  double sinangle = std::sin(arg);
1145 
1146  // Set matrix elements
1147  (*this)(0,0) = cosangle;
1148  (*this)(0,1) = -sinangle;
1149  (*this)(0,2) = 0.0;
1150  (*this)(1,0) = sinangle;
1151  (*this)(1,1) = cosangle;
1152  (*this)(1,2) = 0.0;
1153  (*this)(2,0) = 0.0;
1154  (*this)(2,1) = 0.0;
1155  (*this)(2,2) = 1.0;
1156 
1157  // Return
1158  return;
1159 }
1160 
1161 
1162 /***********************************************************************//**
1163  * @brief Print matrix
1164  *
1165  * @return String containing matrix information.
1166  ***************************************************************************/
1167 std::string GMatrix::print(const GChatter& chatter) const
1168 {
1169  // Initialise result string
1170  std::string result;
1171 
1172  // Continue only if chatter is not silent
1173  if (chatter != SILENT) {
1174 
1175  // Append header
1176  result.append("=== GMatrix ===");
1177 
1178  // Append information
1179  result.append("\n"+gammalib::parformat("Number of rows"));
1180  result.append(gammalib::str(m_rows));
1181  if (m_rowsel != NULL) {
1182  result.append(" (compressed "+gammalib::str(m_num_rowsel)+")");
1183  }
1184  result.append("\n"+gammalib::parformat("Number of columns"));
1185  result.append(gammalib::str(m_cols));
1186  if (m_colsel != NULL) {
1187  result.append(" (compressed "+gammalib::str(m_num_colsel)+")");
1188  }
1189  result.append("\n"+gammalib::parformat("Number of elements"));
1190  result.append(gammalib::str(m_elements));
1191  result.append("\n"+gammalib::parformat("Number of allocated cells"));
1192  result.append(gammalib::str(m_alloc));
1193 
1194  // Append elements and compression schemes
1195  result.append(print_elements(chatter));
1196  result.append(print_row_compression(chatter));
1197  result.append(print_col_compression(chatter));
1198 
1199  } // endif: chatter was not silent
1200 
1201  // Return result
1202  return result;
1203 }
1204 
1205 
1206 /*==========================================================================
1207  = =
1208  = Private methods =
1209  = =
1210  ==========================================================================*/
1211 
1212 /***********************************************************************//**
1213  * @brief Initialise class members
1214  ***************************************************************************/
1216 {
1217  // Return
1218  return;
1219 }
1220 
1221 
1222 /***********************************************************************//**
1223  * @brief Copy class members
1224  *
1225  * @param[in] matrix Matrix.
1226  ***************************************************************************/
1227 void GMatrix::copy_members(const GMatrix& matrix)
1228 {
1229  // Return
1230  return;
1231 }
1232 
1233 
1234 /***********************************************************************//**
1235  * @brief Delete class members
1236  ***************************************************************************/
1238 {
1239  // Return
1240  return;
1241 }
1242 
1243 
1244 /***********************************************************************//**
1245  * @brief Allocate matrix memory
1246  *
1247  * @param[in] rows Number of matrix rows.
1248  * @param[in] columns Number of matrix columns.
1249  *
1250  * Allocates memory for the matrix elements. The method assumes that no
1251  * memory has been allocated for the matrix elements and the column start
1252  * index array. The method allocates the memory for matrix elements and the
1253  * column start indices, sets all matrix elements to 0, and sets the column
1254  * start indices.
1255  *
1256  * If either the number of rows or the number of columns is non-positive, the
1257  * method does nothing.
1258  ***************************************************************************/
1259 void GMatrix::alloc_members(const int& rows, const int& columns)
1260 {
1261  // Determine number of elements to store in matrix
1262  int elements = rows * columns;
1263 
1264  // Continue only if number of elements is positive
1265  if (elements > 0) {
1266 
1267  // Free any existing memory
1268  if (m_data != NULL) delete [] m_data;
1269  if (m_colstart != NULL) delete [] m_colstart;
1270 
1271  // Allocate matrix array and column start index array.
1272  m_data = new double[elements];
1273  m_colstart = new int[columns+1];
1274 
1275  // Store matrix size (logical, storage, allocated)
1276  m_rows = rows;
1277  m_cols = columns;
1278  m_elements = elements;
1279  m_alloc = elements;
1280 
1281  // Set-up column start indices
1282  m_colstart[0] = 0;
1283  for (int col = 1; col <= m_cols; ++col) {
1284  m_colstart[col] = m_colstart[col-1] + m_rows;
1285  }
1286 
1287  // Initialise matrix elements to 0
1288  for (int i = 0; i < m_elements; ++i) {
1289  m_data[i] = 0.0;
1290  }
1291 
1292  } // endif: number of elements was positive
1293 
1294  // Return
1295  return;
1296 }
#define G_AT
Definition: GMatrix.cpp:46
Abstract matrix base class interface definition.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print matrix.
Definition: GMatrix.cpp:1167
std::string print_row_compression(const GChatter &chatter=NORMAL) const
Print row compression scheme if it exists.
void init_members(void)
Initialise class members.
Sparse matrix class interface definition.
GMatrix invert(void) const
Return inverted matrix.
Definition: GMatrix.cpp:876
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
#define G_SET_COLUMN
Definition: GMatrix.cpp:50
virtual double & at(const int &row, const int &column)
Return reference to matrix element.
Definition: GMatrix.cpp:549
#define G_OP_ADD
Definition: GMatrix.cpp:43
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
Generic matrix class definition.
#define G_OP_SUB
Definition: GMatrix.cpp:44
virtual void add_to_row(const int &row, const GVector &vector)
Add row to matrix elements.
Definition: GMatrix.cpp:765
int m_rows
Number of rows.
virtual GMatrix operator-(void) const
Negate matrix elements.
Definition: GMatrix.cpp:316
#define G_ADD_TO_COLUMN
Definition: GMatrix.cpp:52
virtual GVector column(const int &column) const
Extract column as vector from matrix.
Definition: GMatrix.cpp:676
Symmetric matrix class interface definition.
Gammalib tools definition.
#define G_EXTRACT_COLUMN
Definition: GMatrix.cpp:49
#define G_SET_ROW
Definition: GMatrix.cpp:48
int m_num_colsel
Number of selected columns (for compressed decomposition)
void copy_members(const GMatrix &matrix)
Copy class members.
Definition: GMatrix.cpp:1227
virtual void add_to_column(const int &column, const GVector &vector)
Add vector column into matrix.
Definition: GMatrix.cpp:806
#define G_INVERT
Definition: GMatrix.cpp:53
virtual GVector row(const int &row) const
Extract row as vector from matrix.
Definition: GMatrix.cpp:599
void eulery(const double &angle)
Set Euler rotation matrix around y axis.
Definition: GMatrix.cpp:1093
virtual void clear(void)
Clear matrix.
Definition: GMatrix.cpp:511
void init_members(void)
Initialise class members.
Definition: GMatrix.cpp:1215
virtual GMatrix & operator=(const GMatrix &matrix)
Matrix assignment operator.
Definition: GMatrix.cpp:224
#define G_ADD_TO_ROW
Definition: GMatrix.cpp:51
const double deg2rad
Definition: GMath.hpp:43
#define G_EXTRACT_LOWER
Definition: GMatrix.cpp:55
#define G_OP_MUL_VEC
Definition: GMatrix.cpp:42
std::string print_elements(const GChatter &chatter=NORMAL, const int &num=15) const
Print all matrix elements.
#define G_SOLVE
Definition: GMatrix.cpp:54
void eulerx(const double &angle)
Set Euler rotation matrix around x axis.
Definition: GMatrix.cpp:1056
double * m_data
Matrix data.
const int & rows(void) const
Return number of matrix rows.
GMatrix transpose(void) const
Return transposed matrix.
Definition: GMatrix.cpp:847
GVector solve(const GVector &vector) const
Solves linear matrix equation.
Definition: GMatrix.cpp:906
void eulerz(const double &angle)
Set Euler rotation matrix around z axis.
Definition: GMatrix.cpp:1130
#define G_EXTRACT_UPPER
Definition: GMatrix.cpp:56
virtual GVector operator*(const GVector &vector) const
Vector multiplication.
Definition: GMatrix.cpp:284
int m_cols
Number of columns.
virtual GMatrix & operator-=(const GMatrix &matrix)
Unary matrix subtraction operator.
Definition: GMatrix.cpp:395
GChatter
Definition: GTypemaps.hpp:33
int m_alloc
Size of allocated matrix memory.
std::string print_col_compression(const GChatter &chatter=NORMAL) const
Print column compression scheme if it exists.
Vector class interface definition.
virtual double fill(void) const
Return fill of matrix.
Definition: GMatrix.cpp:955
virtual GMatrix & operator*=(const GMatrix &matrix)
Unary matrix multiplication operator.
Definition: GMatrix.cpp:452
#define G_OP_MAT_MUL
Definition: GMatrix.cpp:45
virtual GMatrix * clone(void) const
Clone matrix.
Definition: GMatrix.cpp:529
Symmetric matrix class definition.
void free_members(void)
Delete class members.
Definition: GMatrix.cpp:1237
int m_elements
Number of elements stored in matrix.
GMatrix(void)
Void matrix constructor.
Definition: GMatrix.cpp:73
virtual ~GMatrix(void)
Destructor.
Definition: GMatrix.cpp:200
virtual double sum(void) const
Return matrix element sum.
Definition: GMatrix.hpp:354
#define G_EXTRACT_ROW
Definition: GMatrix.cpp:47
GMatrix extract_lower_triangle(void) const
Extract lower triangle of matrix.
Definition: GMatrix.cpp:993
Sparse matrix class definition.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
GMatrix abs(void) const
Return absolute of matrix.
Definition: GMatrix.cpp:927
Exception handler interface definition.
virtual GMatrixBase & operator=(const GMatrixBase &matrix)
Assignment operator.
const int & size(void) const
Return size of vector.
Definition: GVector.hpp:180
void free_members(void)
Delete class members.
Generic matrix class definition.
Definition: GMatrix.hpp:79
Vector class.
Definition: GVector.hpp:46
int m_num_rowsel
Number of selected rows (for compressed decomposition)
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
int * m_colstart
Column start indices (m_cols+1)
int * m_rowsel
Row selection (for compressed decomposition)
const int & columns(void) const
Return number of matrix columns.
virtual GMatrix & operator+=(const GMatrix &matrix)
Unary matrix addition operator.
Definition: GMatrix.cpp:343
#define G_CONSTRUCTOR
Definition: GMatrix.cpp:41
int * m_colsel
Column selection (for compressed decomposition)
GMatrix extract_upper_triangle(void) const
Extract upper triangle of matrix.
Definition: GMatrix.cpp:1028
void alloc_members(const int &rows, const int &columns)
Allocate matrix memory.
Definition: GMatrix.cpp:1259
Mathematical function definitions.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413