GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GMatrixSymmetric.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GMatrixSymmetric.cpp - Symmetric 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 GMatrixSymmetric.cpp
23  * @brief Symmetric 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 "GTools.hpp"
33 #include "GException.hpp"
34 #include "GVector.hpp"
35 #include "GMatrix.hpp"
36 #include "GMatrixSparse.hpp"
37 #include "GMatrixSymmetric.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 #define G_CONSTRUCTOR "GMatrixSymmetric::GMatrixSymmetric(int&, int&)"
41 #define G_MATRIX "GMatrixSymmetric::GMatrixSymmetric(GMatrix&)"
42 #define G_SPARSEMATRIX "GMatrixSymmetric::GMatrixSymmetric(GSparseMatrix&)"
43 #define G_OP_ADD "GMatrixSymmetric::operator+=(GMatrixSymmetric&)"
44 #define G_OP_SUB "GMatrixSymmetric::operator-=(GMatrixSymmetric&)"
45 #define G_OP_MUL_VEC "GMatrixSymmetric::operator*(GVector&)"
46 #define G_OP_MAT_MUL "GMatrixSymmetric::operator*=(GMatrixSymmetric&)"
47 #define G_AT "GMatrixSymmetric::at(int&, int&)"
48 #define G_EXTRACT_ROW "GMatrixSymmetric::row(int&)"
49 #define G_SET_ROW "GMatrixSymmetric::row(int&, GVector&)"
50 #define G_EXTRACT_COLUMN "GMatrixSymmetric::column(int&)"
51 #define G_SET_COLUMN "GMatrixSymmetric::column(int&, GVector&)"
52 #define G_ADD_TO_ROW "GMatrixSymmetric::add_to_row(int&, GVector&)"
53 #define G_ADD_TO_COLUMN "GMatrixSymmetric::add_to_column(int&, GVector&)"
54 #define G_CHOL_DECOMP "GMatrixSymmetric::cholesky_decompose(int&)"
55 #define G_CHOL_SOLVE "GMatrixSymmetric::cholesky_solver(GVector&, int&)"
56 #define G_CHOL_INVERT "GMatrixSymmetric::cholesky_invert(int&)"
57 #define G_COPY_MEMBERS "GMatrixSymmetric::copy_members(GMatrixSymmetric&)"
58 #define G_ALLOC_MEMBERS "GMatrixSymmetric::alloc_members(int&, int&)"
59 
60 
61 /*==========================================================================
62  = =
63  = Constructors/destructors =
64  = =
65  ==========================================================================*/
66 
67 /***********************************************************************//**
68  * @brief Void constructor
69  *
70  * Constructs an empty matrix (i.e. a matrix without any elements; the number
71  * of rows and columns of the matrix will be zero).
72  *
73  * @todo Check the entire matrix class to see whether handling with an
74  * empty matrix may lead to some conflicts.
75  ***************************************************************************/
77 {
78  // Initialise class members for clean destruction
79  init_members();
80 
81  // Return
82  return;
83 }
84 
85 /***********************************************************************//**
86  * @brief Matrix constructor
87  *
88  * @param[in] rows Number of rows [>0].
89  * @param[in] columns Number of columns [>0].
90  *
91  * @exception GException::empty
92  * Specified number of rows or columns is not valid.
93  ***************************************************************************/
94 GMatrixSymmetric::GMatrixSymmetric(const int& rows, const int& columns) :
95  GMatrixBase()
96 {
97  // Continue only if matrix is valid
98  if (rows > 0 && columns > 0) {
99 
100  // Initialise class members for clean destruction
101  init_members();
102 
103  // Allocate matrix memory
104  alloc_members(rows, columns);
105 
106  }
107  else {
109  }
110 
111  // Return
112  return;
113 }
114 
115 
116 /***********************************************************************//**
117  * @brief GMatrix to GMatrixSymmetric storage class convertor
118  *
119  * @param[in] matrix General matrix (GMatrix).
120  *
121  * @exception GException::matrix_not_symmetric
122  * Matrix is not symmetric.
123  *
124  * Converts a general matrix into the symmetric storage class. If the input
125  * matrix is not symmetric, an exception is thrown.
126  ***************************************************************************/
128 {
129  // Initialise class members for clean destruction
130  init_members();
131 
132  // Allocate matrix memory
133  alloc_members(matrix.rows(), matrix.columns());
134 
135  // Fill matrix
136  for (int col = 0; col < m_cols; ++col) {
137  for (int row = col; row < m_rows; ++row) {
138  double value_ll = matrix(row,col);
139  double value_ur = matrix(col,row);
140  if (value_ll != value_ur) {
142  matrix.rows(),
143  matrix.columns());
144  }
145  (*this)(row, col) = matrix(row, col);
146  }
147  }
148 
149  // Return
150  return;
151 }
152 
153 
154 /***********************************************************************//**
155  * @brief GMatrixSparse to GMatrixSymmetric storage class convertor
156  *
157  * @param[in] matrix Sparse matrix (GMatrixSparse).
158  *
159  * @exception GException::matrix_not_symmetric
160  * Sparse matrix is not symmetric.
161  *
162  * Converts a sparse matrix into the symmetric storage class. If the input
163  * matrix is not symmetric, an exception is thrown.
164  ***************************************************************************/
166 {
167  // Initialise class members for clean destruction
168  init_members();
169 
170  // Allocate matrix memory
171  alloc_members(matrix.rows(), matrix.columns());
172 
173  // Fill matrix
174  for (int col = 0; col < m_cols; ++col) {
175  for (int row = col; row < m_rows; ++row) {
176  double value_ll = matrix(row,col);
177  double value_ur = matrix(col,row);
178  if (value_ll != value_ur) {
180  matrix.rows(),
181  matrix.columns());
182  }
183  (*this)(row, col) = matrix(row, col);
184  }
185  }
186 
187  // Return
188  return;
189 }
190 
191 
192 /***********************************************************************//**
193  * @brief Copy constructor
194  *
195  * @param[in] matrix Matrix.
196  ***************************************************************************/
198  GMatrixBase(matrix)
199 {
200  // Initialise private members for clean destruction
201  init_members();
202 
203  // Copy members
204  copy_members(matrix);
205 
206  // Return
207  return;
208 }
209 
210 
211 /***********************************************************************//**
212  * @brief Destructor
213  ***************************************************************************/
215 {
216  // Free members
217  free_members();
218 
219  // Return
220  return;
221 }
222 
223 
224 /*==========================================================================
225  = =
226  = Operators =
227  = =
228  ==========================================================================*/
229 
230 /***********************************************************************//**
231  * @brief Matrix assignment operator
232  *
233  * @param[in] matrix Matrix.
234  * @return Matrix.
235  *
236  * Assigns the content of another matrix to the actual matrix instance.
237  ***************************************************************************/
239 {
240  // Execute only if object is not identical
241  if (this != &matrix) {
242 
243  // Assign base class members. Note that this method will also perform
244  // the allocation of the matrix memory and the copy of the matrix
245  // attributes.
246  this->GMatrixBase::operator=(matrix);
247 
248  // Free members
249  free_members();
250 
251  // Initialise members
252  init_members();
253 
254  // Copy members
255  copy_members(matrix);
256 
257  } // endif: object was not identical
258 
259  // Return this object
260  return *this;
261 }
262 
263 
264 /***********************************************************************//**
265  * @brief Value assignment operator
266  *
267  * @param[in] value Value.
268  * @return Matrix.
269  *
270  * Assigns the specified @p value to all elements of the matrix.
271  ***************************************************************************/
273 {
274  // Assign value
275  double* ptr = m_data;
276  for (int i = 0; i < m_elements; ++i) {
277  *ptr++ = value;
278  }
279 
280  // Return this object
281  return *this;
282 }
283 
284 
285 /***********************************************************************//**
286  * @brief Return reference to matrix element
287  *
288  * @param[in] row Matrix row [0,...,rows()-1].
289  * @param[in] column Matrix column [0,...,columns()-1].
290  * @return Reference to matrix element.
291  ***************************************************************************/
292 double& GMatrixSymmetric::operator()(const int& row, const int& column)
293 {
294  // Get element index
295  int inx = (row >= column) ? m_colstart[column]+(row-column)
296  : m_colstart[row]+(column-row);
297 
298  // Return element
299  return m_data[inx];
300 }
301 
302 
303 /***********************************************************************//**
304  * @brief Return reference to matrix element (const version)
305  *
306  * @param[in] row Matrix row [0,...,rows()-1].
307  * @param[in] column Matrix column [0,...,columns()-1].
308  * @return Reference to matrix element.
309  ***************************************************************************/
310 const double& GMatrixSymmetric::operator()(const int& row,
311  const int& column) const
312 {
313  // Get element index
314  int inx = (row >= column) ? m_colstart[column]+(row-column)
315  : m_colstart[row]+(column-row);
316 
317  // Return element
318  return m_data[inx];
319 }
320 
321 
322 /***********************************************************************//**
323  * @brief Vector multiplication
324  *
325  * @param[in] vector Vector.
326  *
327  * @exception GException::matrix_vector_mismatch
328  * Vector length differs from number of columns in matrix.
329  *
330  * This method performs a vector multiplication of a matrix. The vector
331  * multiplication will produce a vector. The matrix multiplication can only
332  * be performed when the number of matrix columns is equal to the length of
333  * the vector.
334  ***************************************************************************/
336 {
337  // Raise an exception if the matrix and vector dimensions are not compatible
338  if (m_cols != vector.size()) {
340  m_rows, m_cols);
341  }
342 
343  // Perform vector multiplication
344  GVector result(m_rows);
345  for (int row = 0; row < m_rows; ++row) {
346  double sum = 0.0;
347  for (int col = 0; col < m_cols; ++col) {
348  sum += (*this)(row,col) * vector[col];
349  }
350  result[row] = sum;
351  }
352 
353  // Return result
354  return result;
355 }
356 
357 
358 /***********************************************************************//**
359  * @brief Negate matrix elements
360  *
361  * @return Matrix with negated elements.
362  *
363  * Returns a matrix where each element has been replaced by its negative
364  * element.
365  ***************************************************************************/
367 {
368  // Copy matrix
369  GMatrixSymmetric matrix = *this;
370 
371  // Take the absolute value of all matrix elements
372  double* ptr = matrix.m_data;
373  for (int i = 0; i < matrix.m_elements; ++i, ++ptr) {
374  *ptr = -(*ptr);
375  }
376 
377  // Return matrix
378  return matrix;
379 }
380 
381 
382 /***********************************************************************//**
383  * @brief Unary matrix addition operator
384  *
385  * @param[in] matrix Matrix.
386  *
387  * @exception GException::matrix_mismatch
388  * Incompatible matrix size.
389  *
390  * This method performs a matrix addition. The operation can only succeed
391  * when the dimensions of both matrices are identical.
392  ***************************************************************************/
394 {
395  // Raise an exception if the matrix dimensions are not compatible
396  if (m_rows != matrix.m_rows || m_cols != matrix.m_cols) {
398  m_rows, m_cols,
399  matrix.m_rows, matrix.m_cols);
400  }
401 
402  // Add all matrix elements
403  const double* src = matrix.m_data;
404  double* dst = m_data;
405  for (int i = 0; i < m_elements; ++i) {
406  *dst++ += *src++;
407  }
408 
409  // Return result
410  return *this;
411 }
412 
413 
414 /***********************************************************************//**
415  * @brief Unary matrix scalar addition operator
416  *
417  * @param[in] scalar Scalar.
418  *
419  * Adds a @p scalar to each matrix element.
420  ***************************************************************************/
422 {
423  // Add all matrix elements
424  double* dst = m_data;
425  for (int i = 0; i < m_elements; ++i) {
426  *dst++ += scalar;
427  }
428 
429  // Return result
430  return *this;
431 }
432 
433 
434 /***********************************************************************//**
435  * @brief Unary matrix subtraction operator
436  *
437  * @param[in] matrix Matrix.
438  *
439  * @exception GException::matrix_mismatch
440  * Incompatible matrix size.
441  *
442  * This method performs a matrix addition. The operation can only succeed
443  * when the dimensions of both matrices are identical.
444  ***************************************************************************/
446 {
447  // Raise an exception if the matrix dimensions are not compatible
448  if (m_rows != matrix.m_rows || m_cols != matrix.m_cols) {
450  m_rows, m_cols,
451  matrix.m_rows, matrix.m_cols);
452  }
453 
454  // Subtract all matrix elements
455  const double* src = matrix.m_data;
456  double* dst = m_data;
457  for (int i = 0; i < m_elements; ++i) {
458  *dst++ -= *src++;
459  }
460 
461  // Return result
462  return *this;
463 }
464 
465 
466 /***********************************************************************//**
467  * @brief Unary matrix scalar subtraction operator
468  *
469  * @param[in] scalar Scalar.
470  *
471  * Subtracts a @p scalar from each matrix element.
472  ***************************************************************************/
474 {
475  // Add all matrix elements
476  double* dst = m_data;
477  for (int i = 0; i < m_elements; ++i) {
478  *dst++ -= scalar;
479  }
480 
481  // Return result
482  return *this;
483 }
484 
485 
486 /*==========================================================================
487  = =
488  = Public methods =
489  = =
490  ==========================================================================*/
491 
492 /***********************************************************************//**
493  * @brief Clear matrix
494  ***************************************************************************/
496 {
497  // Free members
498  free_members();
499 
500  // Initialise private members
501  init_members();
502 
503  // Return
504  return;
505 }
506 
507 
508 /***********************************************************************//**
509  * @brief Clone matrix
510  *
511  * @return Pointer to deep copy of matrix.
512  ***************************************************************************/
514 {
515  // Clone matrix
516  return new GMatrixSymmetric(*this);
517 }
518 
519 
520 /***********************************************************************//**
521  * @brief Return reference to matrix element
522  *
523  * @param[in] row Matrix row [0,...,rows()-1].
524  * @param[in] column Matrix column [0,...,columns()-1].
525  * @return Reference to matrix element.
526  *
527  * @exception GException::out_of_range
528  * Row or column index out of range.
529  ***************************************************************************/
530 double& GMatrixSymmetric::at(const int& row, const int& column)
531 {
532  // Raise exception if row or column index is out of range
533  if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) {
534  throw GException::out_of_range(G_AT, row, column, m_rows, m_cols);
535  }
536 
537  // Get element index
538  int inx = (row >= column) ? m_colstart[column]+(row-column)
539  : m_colstart[row]+(column-row);
540 
541  // Return element
542  return m_data[inx];
543 }
544 
545 
546 /***********************************************************************//**
547  * @brief Return reference to matrix element (const version)
548  *
549  * @param[in] row Matrix row [0,...,rows()-1].
550  * @param[in] column Matrix column [0,...,columns()-1].
551  * @return Reference to matrix element.
552  *
553  * @exception GException::out_of_range
554  * Row or column index out of range.
555  ***************************************************************************/
556 const double& GMatrixSymmetric::at(const int& row, const int& column) const
557 {
558  // Raise exception if row or column index is out of range
559  if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) {
560  throw GException::out_of_range(G_AT, row, column, m_rows, m_cols);
561  }
562 
563  // Get element index
564  int inx = (row >= column) ? m_colstart[column]+(row-column)
565  : m_colstart[row]+(column-row);
566 
567  // Return element
568  return m_data[inx];
569 }
570 
571 
572 /***********************************************************************//**
573  * @brief Extract row as vector from matrix
574  *
575  * @param[in] row Row to be extracted (starting from 0).
576  *
577  * @exception GException::out_of_range
578  * Invalid row index specified.
579  *
580  * This method extracts a matrix row into a vector.
581  ***************************************************************************/
582 GVector GMatrixSymmetric::row(const int& row) const
583 {
584  // Raise an exception if the row index is invalid
585  #if defined(G_RANGE_CHECK)
586  if (row < 0 || row >= m_rows) {
588  }
589  #endif
590 
591  // Create result vector
592  GVector result(m_cols);
593 
594  // Extract row into vector
595  for (int col = 0; col < m_cols; ++col) {
596  result[col] = (*this)(row,col);
597  }
598 
599  // Return vector
600  return result;
601 }
602 
603 
604 /***********************************************************************//**
605  * @brief Set row in matrix
606  *
607  * @todo To be implemented.
608  ***************************************************************************/
609 void GMatrixSymmetric::row(const int& row, const GVector& vector)
610 {
611  // Raise an exception if the row index is invalid
612  #if defined(G_RANGE_CHECK)
613  if (row < 0 || row >= m_rows) {
614  throw GException::out_of_range(G_SET_ROW, row, 0, m_rows-1);
615  }
616  #endif
617 
618  // Return
619  return;
620 }
621 
622 
623 /***********************************************************************//**
624  * @brief Extract column as vector from matrix
625  *
626  * @param[in] column Column index [0,...,columns()-1].
627  *
628  * @exception GException::out_of_range
629  * Invalid row index specified.
630  *
631  * This method extracts a matrix column into a vector.
632  ***************************************************************************/
633 GVector GMatrixSymmetric::column(const int& column) const
634 {
635  // Raise an exception if the column index is invalid
636  #if defined(G_RANGE_CHECK)
637  if (column < 0 || column >= m_cols) {
638  throw GException::out_of_range(G_EXTRACT_COLUMN, column, 0, m_cols-1);
639  }
640  #endif
641 
642  // Create result vector
643  GVector result(m_rows);
644 
645  // Extract column into vector
646  for (int row = 0; row < m_rows; ++row) {
647  result[row] = (*this)(row, column);
648  }
649 
650  // Return vector
651  return result;
652 }
653 
654 
655 /***********************************************************************//**
656  * @brief Set matrix column from vector
657  *
658  * @param[in] column Column index [0,...,columns()-1].
659  * @param[in] vector Vector.
660  *
661  * @exception GException::out_of_range
662  * Invalid column index specified.
663  * @exception GException::matrix_vector_mismatch
664  * Matrix dimension mismatches the vector size.
665  *
666  * Inserts the content of a vector into a matrix column. Any previous
667  * content in the matrix column will be overwritted.
668  ***************************************************************************/
669 void GMatrixSymmetric::column(const int& column, const GVector& vector)
670 {
671  // Raise an exception if the column index is invalid
672  #if defined(G_RANGE_CHECK)
673  if (column < 0 || column >= m_cols) {
674  throw GException::out_of_range(G_SET_COLUMN, column, 0, m_cols-1);
675  }
676  #endif
677 
678  // Raise an exception if the matrix and vector dimensions are not
679  // compatible
680  if (m_rows != vector.size()) {
682  m_rows, m_cols);
683  }
684 
685  // Insert column into vector
686  for (int row = 0; row < m_rows; ++row) {
687  (*this)(row, column) = vector[row];
688  }
689 
690  // Return
691  return;
692 }
693 
694 
695 /***********************************************************************//**
696  * @brief Add row to matrix elements
697  *
698  * @todo To be implemented.
699  ***************************************************************************/
700 void GMatrixSymmetric::add_to_row(const int& row, const GVector& vector)
701 {
702  // Raise an exception if the row index is invalid
703  #if defined(G_RANGE_CHECK)
704  if (row < 0 || row >= m_rows) {
705  throw GException::out_of_range(G_ADD_TO_ROW, row, 0, m_rows-1);
706  }
707  #endif
708 
709  // Return
710  return;
711 }
712 
713 
714 /***********************************************************************//**
715  * @brief Add vector column into matrix
716  *
717  * @param[in] column Column index [0,...,columns()-1].
718  * @param[in] vector Vector.
719  *
720  * @exception GException::out_of_range
721  * Invalid column index specified.
722  * @exception GException::matrix_vector_mismatch
723  * Matrix dimension mismatches the vector size.
724  *
725  * Adds the content of a vector to a matrix column.
726  ***************************************************************************/
727 void GMatrixSymmetric::add_to_column(const int& column, const GVector& vector)
728 {
729  // Raise an exception if the column index is invalid
730  #if defined(G_RANGE_CHECK)
731  if (column < 0 || column >= m_cols) {
732  throw GException::out_of_range(G_ADD_TO_COLUMN, column, 0, m_cols-1);
733  }
734  #endif
735 
736  // Raise an exception if the matrix and vector dimensions are not
737  // compatible
738  if (m_rows != vector.size()) {
740  m_rows, m_cols);
741  }
742 
743  // Insert column into vector
744  for (int row = 0; row < m_rows; ++row) {
745  (*this)(row, column) += vector[row];
746  }
747 
748  // Return
749  return;
750 }
751 
752 
753 /***********************************************************************//**
754  * @brief Return inverted matrix
755  *
756  * @return Inverted matrix.
757  *
758  * Returns inverse of matrix. Inversion is done for the moment using Cholesky
759  * decomposition. This does not work on any kind of matrix.
760  *
761  * @todo Specify in documentation for which kind of matrix the method works.
762  ***************************************************************************/
764 {
765  // Allocate result matrix
766  GMatrixSymmetric matrix(m_cols, m_rows);
767 
768  // Invert matrix
769  matrix.cholesky_invert(true);
770 
771  // Return matrix
772  return matrix;
773 }
774 
775 
776 /***********************************************************************//**
777  * @brief Solves linear matrix equation
778  *
779  * @param[in] vector Solution vector.
780  *
781  * Solves the linear equation
782  *
783  * \f[M \times {\tt solution} = {\tt vector} \f]
784  *
785  * where \f$M\f$ is the matrix, \f${\tt vector}\f$ is the result, and
786  * \f${\tt solution}\f$ is the solution. Solving is done using Cholesky
787  * decomposition. This does not work on any kind of matrix.
788  *
789  * @todo Specify in documentation for which kind of matrix the method works.
790  ***************************************************************************/
792 {
793  // Get Cholesky decomposition of matrix
794  GMatrixSymmetric decomposition = cholesky_decompose(true);
795 
796  // Solve linear equation
797  GVector result = decomposition.cholesky_solver(vector);
798 
799  // Return result
800  return result;
801 }
802 
803 
804 /***********************************************************************//**
805  * @brief Return absolute of matrix
806  *
807  * @return Absolute of matrix
808  *
809  * Returns matrix where all elements of the matrix have been replaced by
810  * their absolute values.
811  ***************************************************************************/
813 {
814  // Allocate result matrix
815  GMatrixSymmetric matrix(m_rows, m_cols);
816 
817  // Take the absolute value of all matrix elements
818  double* src = m_data;
819  double* dst = matrix.m_data;
820  for (int i = 0; i < m_elements; ++i) {
821  *dst++ = std::abs(*src++);
822  }
823 
824  // Return matrix
825  return matrix;
826 }
827 
828 
829 /***********************************************************************//**
830  * @brief Determine fill of matrix
831  *
832  * @return Matrix fill (between 0 and 1).
833  *
834  * The fill of a matrix is defined as the number of non-zero elements
835  * devided by the total number of matrix elements.
836  ***************************************************************************/
837 double GMatrixSymmetric::fill(void) const
838 {
839  // Determine the number of zero elements
840  int zero = 0;
841  for (int col = 0, i = 0; col < m_cols; ++col) {
842  if (m_data[i++] == 0.0) { // Count diag. once
843  zero++;
844  }
845  for (int row = col+1; row < m_rows; ++row) {
846  if (m_data[i++] == 0.0) { // Count off-diag. twice
847  zero +=2;
848  }
849  }
850  }
851 
852  // Return the fill
853  return (1.0-double(zero)/double(m_elements));
854 }
855 
856 
857 /***********************************************************************//**
858  * @brief Sum matrix elements
859  *
860  * @return Sum of all matrix elements.
861  ***************************************************************************/
862 double GMatrixSymmetric::sum(void) const
863 {
864  // Initialise matrix sums (diagonal and off-diagonal)
865  double diag = 0.0;
866  double off_diag = 0.0;
867 
868  // Calulate sum over diagonal elements
869  for (int row = 0; row < m_rows; ++row) {
870  diag += m_data[m_colstart[row]];
871  }
872 
873  // Calulate sum over off-diagonal elements
874  for (int row = 0; row < m_rows; ++row) {
875  for (int col = row+1; col < m_cols; ++col) {
876  off_diag += m_data[m_colstart[row]+(col-row)];
877  }
878  }
879 
880  // Calculate total
881  double result = diag + 2.0 * off_diag;
882 
883  // Return result
884  return result;
885 }
886 
887 
888 /***********************************************************************//**
889  * @brief Extract lower triangle of matrix
890  *
891  * This method extracts the lower triangle of a matrix into another matrix.
892  * (including the diagonal elements).
893  * All remaining matrix elements will be zero.
894  ***************************************************************************/
896 {
897  // Define result matrix
898  GMatrix result(m_rows, m_cols);
899 
900  // Extract all elements
901  for (int row = 0; row < m_rows; ++row) {
902  for (int col = 0; col <= row; ++col) {
903  result(row,col) = m_data[m_colstart[col]+(row-col)];
904  }
905  }
906 
907  // Return result
908  return result;
909 }
910 
911 
912 /***********************************************************************//**
913  * @brief Extract upper triangle of matrix
914  *
915  * This method extracts the upper triangle of a matrix into another matrix.
916  * (including the diagonal elements).
917  * All remaining matrix elements will be zero.
918  ***************************************************************************/
920 {
921  // Define result matrix
922  GMatrix result(m_rows, m_cols);
923 
924  // Extract all elements
925  for (int row = 0; row < m_rows; ++row) {
926  for (int col = row; col < m_cols; ++col) {
927  result(row,col) = m_data[m_colstart[row]+(col-row)];
928  }
929  }
930 
931  // Return result
932  return result;
933 }
934 
935 
936 /***********************************************************************//**
937  * @brief Return Cholesky decomposition
938  *
939  * @param[in] compress Use zero-row/column compression (defaults to true).
940  * @return Cholesky decomposition of matrix
941  *
942  * @exception GException::matrix_not_pos_definite
943  * Matrix is not positive definite.
944  * @exception GException::matrix_zero
945  * All matrix elements are zero.
946  *
947  * Returns the Cholesky decomposition of a sparse matrix. The decomposition
948  * is stored within a GMatrixSymmetric object.
949  *
950  * The method is inspired by the algorithm found in Numerical Recipes.
951  * The decomposition, which is a matrix occupying only the lower triange,
952  * is stored in the elements of the symmetric matrix. To visualise the
953  * matrix one has to use 'lower_triangle()' to extract the relevant part.
954  * Case A operates on a full matrix, Case B operates on a (logically)
955  * compressed matrix where zero rows/columns have been removed.
956  ***************************************************************************/
958 {
959  // Create copy of matrix
960  GMatrixSymmetric matrix = *this;
961 
962  // Set-up incides of non zero rows if matrix compression is requested
963  if (compress) {
964  matrix.set_inx();
965  }
966 
967  // Check if zero-row/col compression is needed
968  int no_zeros = ((compress && (matrix.m_num_inx == matrix.m_rows)) || !compress);
969 
970  // Case A: no zero-row/col compression needed
971  if (no_zeros) {
972 
973  // Loop over upper triangle (col >= row)
974  double diag = 0.0;
975  for (int row = 0; row < matrix.m_rows; ++row) {
976  double* ptr = matrix.m_data + matrix.m_colstart[row];
977  for (int col = row; col < matrix.m_cols; ++col, ++ptr) {
978  // sum = M(row,col)
979  double sum = *ptr;
980  for (int k = 0; k < row; ++k) {
981  int offset = matrix.m_colstart[k] - k; // is always positive
982  sum -= matrix.m_data[offset+row] * matrix.m_data[offset+col]; // sum -= M(row,k)*M(col,k)
983  }
984  if (row == col) {
985  if (sum <= 0.0) {
987  }
988  *ptr = std::sqrt(sum); // M(row,row) = sqrt(sum)
989  diag = 1.0/(*ptr);
990  }
991  else {
992  *ptr = sum*diag; // M(row,col) = sum/M(row,row)
993  }
994  }
995  }
996  } // endif: there were no zero rows/cols in matrix
997 
998  // Case B: zero-row/col compression needed
999  else if (matrix.m_num_inx > 0) {
1000 
1001  // Allocate loop variables and pointers
1002  int row;
1003  int col;
1004  int k;
1005  int* row_ptr;
1006  int* col_ptr;
1007  int* k_ptr;
1008 
1009  // Loop over upper triangle (col >= row)
1010  double diag = 0.0;
1011  for (row = 0, row_ptr = matrix.m_inx; row < matrix.m_num_inx; ++row, ++row_ptr) {
1012  double* ptr_0 = matrix.m_data + matrix.m_colstart[*row_ptr] - *row_ptr;
1013  for (col = row, col_ptr = matrix.m_inx + row; col < matrix.m_num_inx; ++col, ++col_ptr) {
1014  double* ptr = ptr_0 + *col_ptr;
1015  double sum = *ptr; // sum = M(row,col)
1016  for (k = 0, k_ptr = matrix.m_inx; k < row; ++k, ++k_ptr) {
1017  int offset = matrix.m_colstart[*k_ptr] - *k_ptr; // is always positive
1018  sum -= matrix.m_data[offset+*row_ptr] *
1019  matrix.m_data[offset+*col_ptr];
1020  // sum -= M(row,k)*M(col,k)
1021  }
1022  if (*row_ptr == *col_ptr) {
1023  if (sum <= 0.0) {
1024  throw GException::matrix_not_pos_definite(G_CHOL_DECOMP, *row_ptr, sum);
1025  }
1026  *ptr = std::sqrt(sum); // M(row,row) = sqrt(sum)
1027  diag = 1.0/(*ptr);
1028  }
1029  else {
1030  *ptr = sum*diag; // M(row,col) = sum/M(row,row)
1031  }
1032  }
1033  }
1034  } // endelse: zero-row/col compression needed
1035 
1036  // Case C: all matrix elements are zero
1037  else {
1039  }
1040 
1041  // Return matrix
1042  return matrix;
1043 }
1044 
1045 
1046 /***********************************************************************//**
1047  * @brief Cholesky solver
1048  *
1049  * @param[in] vector Vector for which should be solved.
1050  * @param[in] compress Use zero-row/column compression (default: true).
1051  *
1052  * @exception GException::matrix_vector_mismatch
1053  * Matrix and vector do not match.
1054  * @exception GException::matrix_zero
1055  * All matrix elements are zero.
1056  *
1057  * Solves the linear equation A*x=b using a Cholesky decomposition of A.
1058  * This function is to be applied on a decomposition GMatrixSymmetric matrix
1059  * that is produced by 'cholesky_decompose'. Case A operates on a full
1060  * matrix, Case B on a zero rows/columns (logically) compressed matrix.
1061  ***************************************************************************/
1063  const bool& compress) const
1064 {
1065  // Raise an exception if the matrix and vector dimensions are not compatible
1066  if (m_rows != vector.size()) {
1068  m_rows, m_cols);
1069  }
1070 
1071  // Allocate result vector
1072  GVector x(m_rows);
1073 
1074  // Check if zero-row/col compression is needed
1075  int no_zeros = ((compress && (m_num_inx == m_rows)) || !compress);
1076 
1077  // Case A: no zero-row/col compression needed
1078  if (no_zeros) {
1079 
1080  // Solve L*y=b, storing y in x (row>k)
1081  for (int row = 0; row < m_rows; ++row) {
1082  double sum = vector[row];
1083  for (int k = 0; k < row; ++k) {
1084  sum -= m_data[m_colstart[k]+(row-k)] * x[k]; // sum -= M(row,k) * x(k)
1085  }
1086  x[row] = sum/m_data[m_colstart[row]]; // x(row) = sum/M(row,row)
1087  }
1088 
1089  // Solve trans(L)*x=y (k>row)
1090  for (int row = m_rows-1; row >= 0; --row) {
1091  double sum = x[row];
1092  double* ptr = m_data + m_colstart[row] + 1;
1093  for (int k = row+1; k < m_rows; ++k) {
1094  sum -= *ptr++ * x[k]; // sum -= M(k,row) * x(k)
1095  }
1096  x[row] = sum/m_data[m_colstart[row]]; // x(row) = sum/M(row,row)
1097  }
1098  } // endif: no zero-row/col compression needed
1099 
1100  // Case B: zero-row/col compression needed
1101  else if (m_num_inx > 0) {
1102 
1103  // Allocate loop variables and pointers
1104  int row;
1105  int k;
1106  int* row_ptr;
1107  int* k_ptr;
1108 
1109  // Solve L*y=b, storing y in x (row>k)
1110  for (row = 0, row_ptr = m_inx; row < m_num_inx; ++row, ++row_ptr) {
1111  double sum = vector[*row_ptr];
1112  double* ptr = m_data + *row_ptr;
1113  for (k = 0, k_ptr = m_inx; k < row; ++k, ++k_ptr) {
1114  sum -= *(ptr + m_colstart[*k_ptr] - *k_ptr) * x[*k_ptr]; // sum -= M(row,k) * x(k)
1115  }
1116  x[*row_ptr] = sum/m_data[m_colstart[*row_ptr]]; // x(row) = sum/M(row,row)
1117  }
1118 
1119  // Solve trans(L)*x=y (k>row)
1120  for (row = m_num_inx-1, row_ptr = m_inx+m_num_inx-1; row >= 0; --row, --row_ptr) {
1121  double sum = x[*row_ptr];
1122  double* ptr_diag = m_data + m_colstart[*row_ptr];
1123  double* ptr = ptr_diag - *row_ptr;
1124  for (k = row+1, k_ptr = m_inx+row+1; k < m_num_inx; ++k, ++k_ptr) {
1125  sum -= *(ptr + *k_ptr) * x[*k_ptr]; // sum -= M(k,row) * x(k)
1126  }
1127  x[*row_ptr] = sum/(*ptr_diag); // x(row) = sum/M(row,row)
1128  }
1129  } // endelse: zero-row/col compression needed
1130 
1131  // Case C: all matrix elements are zero
1132  else {
1134  }
1135 
1136  // Return result vector
1137  return x;
1138 }
1139 
1140 
1141 /***********************************************************************//**
1142  * @brief Invert matrix using a Cholesky decomposition
1143  *
1144  * @param[in] compress Use zero-row/column compression (defaults to true).
1145  * @return Inverted matrix.
1146  *
1147  * @exception GException::matrix_zero
1148  * All matrix elements are zero.
1149  *
1150  * Inverts the matrix using a Cholesky decomposition.
1151  *
1152  * The method distinguish two cases. Case A operates on a full matrix while
1153  * Case B operates on a (logically) compressed matrix where all zero
1154  * rows/columns are skipped.
1155  ***************************************************************************/
1157 {
1158  // Generate Cholesky decomposition of matrix
1159  GMatrixSymmetric matrix = cholesky_decompose(compress);
1160 
1161  // Check if zero-row/col compression is needed
1162  int no_zeros = ((compress && (matrix.m_num_inx == matrix.m_rows)) || !compress);
1163 
1164  // Case A: no zero-row/col compression needed
1165  if (no_zeros) {
1166 
1167  // Generate inverse of Cholesky decomposition (col>row)
1168  for (int row = 0; row < matrix.m_rows; ++row) {
1169 
1170  // M(row,row) = 1/M(row,row)
1171  double* ptr = matrix.m_data + matrix.m_colstart[row];
1172  *ptr = 1.0/(*ptr);
1173 
1174  for (int col = row+1; col < matrix.m_cols; ++col) {
1175 
1176  // sum -= M(col,k)*M(k,row)
1177  double sum = 0.0;
1178  double* ptr1 = matrix.m_data + col - row;
1179  double* ptr2 = ptr;
1180  for (int k = row; k < col; ++k) {
1181  sum -= *(ptr1-- + matrix.m_colstart[k]) * *ptr2++;
1182  }
1183 
1184  // M(col,row) = sum/M(col,col)
1185  *(ptr+col-row) = sum/matrix.m_data[matrix.m_colstart[col]];
1186  }
1187  }
1188 
1189  // Matrix multiplication (col>=row)
1190  for (int row = 0; row < matrix.m_rows; ++row) {
1191  double* ptr = matrix.m_data + matrix.m_colstart[row];
1192  for (int col = row; col < matrix.m_cols; ++col) {
1193  // sum += M(row,k)*M(k,col)
1194  double sum = 0.0;
1195  double* ptr1 = ptr + col - row;
1196  double* ptr2 = matrix.m_data + matrix.m_colstart[col];
1197  for (int k = col; k < matrix.m_cols; ++k) {
1198  sum += *ptr1++ * *ptr2++;
1199  }
1200  // M(row,col) = sum
1201  *(ptr+col-row) = sum;
1202  }
1203  }
1204  } // endif: no zero-row/col compression needed
1205 
1206  // Case B: zero-row/col compression needed
1207  else if (matrix.m_num_inx > 0) {
1208 
1209  // Allocate loop variables and pointers
1210  int row;
1211  int col;
1212  int k;
1213  int* row_ptr;
1214  int* col_ptr;
1215  int* k_ptr;
1216 
1217  // Generate inverse of Cholesky decomposition (col>row)
1218  for (row = 0, row_ptr = matrix.m_inx;
1219  row < matrix.m_num_inx; ++row, ++row_ptr) {
1220 
1221  // M(row,row) = 1/M(row,row)
1222  double* ptr_diag = matrix.m_data + matrix.m_colstart[*row_ptr];
1223  double* ptr_2 = ptr_diag - *row_ptr;
1224  *ptr_diag = 1.0/(*ptr_diag);
1225 
1226  for (col = row+1, col_ptr = matrix.m_inx+row+1;
1227  col < matrix.m_num_inx; ++col, ++col_ptr) {
1228 
1229  // sum -= M(col,k)*M(k,row)
1230  double sum = 0.0;
1231  double* ptr_1 = matrix.m_data + *col_ptr;
1232  for (k = row, k_ptr = matrix.m_inx+row;
1233  k < col; ++k, ++k_ptr) {
1234  sum -= *(ptr_1 + matrix.m_colstart[*k_ptr] - *k_ptr) *
1235  *(ptr_2 + *k_ptr);
1236  }
1237 
1238  // M(col,row) = sum/M(col,col)
1239  *(ptr_2 + *col_ptr) =
1240  sum/matrix.m_data[matrix.m_colstart[*col_ptr]];
1241  }
1242  }
1243 
1244  // Matrix multiplication (col>=row)
1245  for (row = 0, row_ptr = matrix.m_inx;
1246  row < matrix.m_num_inx; ++row, ++row_ptr) {
1247  double* ptr_diag = matrix.m_data + matrix.m_colstart[*row_ptr];
1248  double* ptr_1 = ptr_diag - *row_ptr;
1249 
1250  for (col = row, col_ptr = matrix.m_inx+row;
1251  col < matrix.m_num_inx; ++col, ++col_ptr) {
1252 
1253  // sum += M(row,k)*M(k,col)
1254  double sum = 0.0;
1255  double* ptr_2 = matrix.m_data + matrix.m_colstart[*col_ptr] - *col_ptr;
1256  for (k = col, k_ptr = matrix.m_inx+col;
1257  k < matrix.m_num_inx; ++k, ++k_ptr) {
1258  sum += *(ptr_1 + *k_ptr) * *(ptr_2 + *k_ptr);
1259  }
1260 
1261  // M(row,col) = sum
1262  *(ptr_1 + *col_ptr) = sum;
1263  }
1264  }
1265  } // endelse: zero-row/col compression needed
1266 
1267  // Case C: all matrix elements are zero
1268  else {
1270  }
1271 
1272  // Return matrix
1273  return matrix;
1274 }
1275 
1276 
1277 /***********************************************************************//**
1278  * @brief Print matrix
1279  *
1280  * @param[in] chatter Chattiness (defaults to NORMAL).
1281  * @return String containing matrix information
1282  ***************************************************************************/
1283 std::string GMatrixSymmetric::print(const GChatter& chatter) const
1284 {
1285  // Initialise result string
1286  std::string result;
1287 
1288  // Append header
1289  result.append("=== GMatrixSymmetric ===");
1290 
1291  // Continue only if chatter is not silent
1292  if (chatter != SILENT) {
1293 
1294  // Append information
1295  result.append("\n"+gammalib::parformat("Number of rows"));
1296  result.append(gammalib::str(m_rows));
1297  if (m_rowsel != NULL) {
1298  result.append(" (compressed "+gammalib::str(m_num_rowsel)+")");
1299  }
1300  result.append("\n"+gammalib::parformat("Number of columns"));
1301  result.append(gammalib::str(m_cols));
1302  if (m_colsel != NULL) {
1303  result.append(" (compressed "+gammalib::str(m_num_colsel)+")");
1304  }
1305  result.append("\n"+gammalib::parformat("Number of elements"));
1306  result.append(gammalib::str(m_elements));
1307  result.append("\n"+gammalib::parformat("Number of allocated cells"));
1308  result.append(gammalib::str(m_alloc));
1309 
1310  // Append elements and compression schemes
1311  result.append(print_elements(chatter));
1312  result.append(print_row_compression(chatter));
1313  result.append(print_col_compression(chatter));
1314 
1315  } // endif: chatter was not silent
1316 
1317  // Return result
1318  return result;
1319 }
1320 
1321 
1322 /*==========================================================================
1323  = =
1324  = Private methods =
1325  = =
1326  ==========================================================================*/
1327 
1328 /***********************************************************************//**
1329  * @brief Initialise class mambers
1330  ***************************************************************************/
1332 {
1333  // Initialise members
1334  m_num_inx = 0;
1335  m_inx = NULL;
1336 
1337  // Return
1338  return;
1339 }
1340 
1341 
1342 /***********************************************************************//**
1343  * @brief Copy class members
1344  *
1345  * @param[in] matrix Matrix.
1346  ***************************************************************************/
1348 {
1349  // Copy attributes
1350  m_num_inx = matrix.m_num_inx;
1351 
1352  // Copy index selection
1353  if (m_cols > 0) {
1354  m_inx = new int[m_cols];
1355  for (int i = 0; i < m_cols; ++i) {
1356  m_inx[i] = matrix.m_inx[i];
1357  }
1358  }
1359 
1360  // Return
1361  return;
1362 }
1363 
1364 
1365 /***********************************************************************//**
1366  * @brief Delete class members
1367  ***************************************************************************/
1369 {
1370  // De-allocate only if memory has indeed been allocated by derived class
1371  if (m_inx != NULL) delete [] m_inx;
1372 
1373  // Return
1374  return;
1375 }
1376 
1377 
1378 /***********************************************************************//**
1379  * @brief Allocates matrix memory
1380  *
1381  * @param[in] rows Number of rows.
1382  * @param[in] columns Number of columns.
1383  *
1384  * @exception GException::matrix_not_symmetric
1385  * Matrix is not symmetric.
1386  *
1387  * This method is the main constructor code that allocates and initialises
1388  * memory for matrix elements. The method assumes that no memory has been
1389  * allocated for the matrix elements, the column start index array and the
1390  * index array.
1391  * The method allocates the memory for matrix elements, the column start
1392  * indices and the index array, sets all matrix elements to 0.0, and sets
1393  * the column start indices. The content of the index array is undefined.
1394  *
1395  * @todo Verify if the index array m_inx should be initialized.
1396  ***************************************************************************/
1397 void GMatrixSymmetric::alloc_members(const int& rows, const int& columns)
1398 {
1399  // Determine number of physical elements in matrix
1400  int elements = rows*(rows+1)/2;
1401 
1402  // Throw exception if number of rows and columns is not identical
1403  if (rows != columns) {
1404  throw GException::matrix_not_symmetric(G_ALLOC_MEMBERS, rows, columns);
1405  }
1406 
1407  // Continue only if number of elements is positive
1408  if (elements > 0) {
1409 
1410  // Free any existing memory
1411  if (m_data != NULL) delete [] m_data;
1412  if (m_colstart != NULL) delete [] m_colstart;
1413  if (m_inx != NULL) delete [] m_inx;
1414 
1415  // Allocate matrix array and column start index array.
1416  m_data = new double[elements];
1417  m_colstart = new int[columns+1];
1418  m_inx = new int[columns];
1419 
1420  // Store matrix size (logical and physical)
1421  m_rows = rows;
1422  m_cols = columns;
1423  m_elements = elements;
1424  m_alloc = elements;
1425 
1426  // Set-up column start indices
1427  m_colstart[0] = 0;
1428  int offset = rows;
1429  for (int col = 1; col <= m_cols; ++col) {
1430  m_colstart[col] = m_colstart[col-1] + offset--;
1431  }
1432 
1433  // Initialise matrix elements to 0.0
1434  for (int i = 0; i < m_elements; ++i) {
1435  m_data[i] = 0.0;
1436  }
1437 
1438  } // endif: number of elements was positive
1439 
1440  // Return
1441  return;
1442 }
1443 
1444 
1445 /***********************************************************************//**
1446  * @brief Set index selection
1447  *
1448  * Determines the non-zero rows/columns in matrix and set up index array
1449  * that points to these rows/columns. This array is used for compressed
1450  * matrix calculations (Case B in the above routines).
1451  ***************************************************************************/
1453 {
1454  // Allocate loop variables and pointers
1455  int row;
1456  int col;
1457 
1458  // Find all nonzero rows/cols
1459  m_num_inx = 0;
1460  for (row = 0; row < m_rows; ++row) {
1461 
1462  // Check first for zero diagonal. If we have one then check the rest of
1463  // the row
1464  if (m_data[m_colstart[row]] == 0.0) {
1465  for (col = 0; col < row; ++col) {
1466  if (m_data[m_colstart[col]+(row-col)] != 0.0) {
1467  break;
1468  }
1469  }
1470  // Found a non-zero element
1471  if (col < row) {
1472  m_inx[m_num_inx++] = row;
1473  }
1474  else {
1475  for (col = row+1; col < m_cols; ++col) {
1476  if (m_data[m_colstart[row]+(col-row)] != 0.0) {
1477  break;
1478  }
1479  }
1480  // Found a non-zero element
1481  if (col < m_cols) {
1482  m_inx[m_num_inx++] = row;
1483  }
1484  }
1485  }
1486  else {
1487  m_inx[m_num_inx++] = row;
1488  }
1489  }
1490 
1491  // Return
1492  return;
1493 }
1494 
1495 
1496 /*==========================================================================
1497  = =
1498  = Friend functions =
1499  = =
1500  ==========================================================================*/
1501 
1502 /***********************************************************************//**
1503  * @brief Binary matrix multiplication operator
1504  *
1505  * @param[in] matrix Matrix to be multiplied.
1506  *
1507  * @exception GException::matrix_mismatch
1508  * Incompatible matrix size.
1509  *
1510  * This method performs a matrix multiplication. Since the product of two
1511  * symmetric matrices is not necessarily symmetric, this method returns a
1512  * GMatrix object.
1513  *
1514  * The operation can only succeed when the dimensions of both matrices are
1515  * compatible.
1516  ***************************************************************************/
1518 {
1519  // Raise an exception if the matrix dimensions are not compatible
1520  if (m_cols != matrix.m_rows) {
1522  m_rows, m_cols,
1523  matrix.m_rows, matrix.m_cols);
1524  }
1525 
1526  // Allocate result matrix
1527  GMatrix result(m_rows, matrix.m_cols);
1528 
1529  // Loop over all elements of result matrix
1530  for (int row = 0; row < m_rows; ++row) {
1531  for (int col = 0; col < matrix.m_cols; ++col) {
1532  double sum = 0.0;
1533  for (int i = 0; i < m_cols; ++i) {
1534  sum += (*this)(row,i) * matrix(i,col);
1535  }
1536  result(row,col) = sum;
1537  }
1538  }
1539 
1540  // Return result
1541  return result;
1542 }
GMatrix extract_lower_triangle(void) const
Extract lower triangle of matrix.
GMatrixSymmetric cholesky_invert(const bool &compress=true) const
Invert matrix using a Cholesky decomposition.
#define G_ALLOC_MEMBERS
#define G_AT
void free_members(void)
Delete class members.
#define G_SPARSEMATRIX
Abstract matrix base class interface definition.
virtual ~GMatrixSymmetric(void)
Destructor.
std::string print_row_compression(const GChatter &chatter=NORMAL) const
Print row compression scheme if it exists.
GMatrixSymmetric(void)
Void constructor.
#define G_CHOL_INVERT
Sparse matrix class interface definition.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
Generic matrix class definition.
#define G_MATRIX
#define G_SET_COLUMN
virtual void add_to_row(const int &row, const GVector &vector)
Add row to matrix elements.
#define G_OP_SUB
int m_rows
Number of rows.
#define G_ADD_TO_COLUMN
Symmetric matrix class interface definition.
virtual double & operator()(const int &row, const int &column)
Return reference to matrix element.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print matrix.
virtual void clear(void)
Clear matrix.
Gammalib tools definition.
#define G_SET_ROW
#define G_EXTRACT_COLUMN
int m_num_inx
Number of indices in array.
GMatrixSymmetric invert(void) const
Return inverted matrix.
int m_num_colsel
Number of selected columns (for compressed decomposition)
virtual GVector column(const int &column) const
Extract column as vector from matrix.
GVector cholesky_solver(const GVector &vector, const bool &compress=true) const
Cholesky solver.
virtual GMatrixSymmetric & operator-=(const GMatrixSymmetric &matrix)
Unary matrix subtraction operator.
virtual GMatrixSymmetric & operator+=(const GMatrixSymmetric &matrix)
Unary matrix addition operator.
virtual GMatrixSymmetric * clone(void) const
Clone matrix.
#define G_ADD_TO_ROW
#define G_CHOL_SOLVE
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
#define G_OP_MUL_VEC
std::string print_elements(const GChatter &chatter=NORMAL, const int &num=15) const
Print all matrix elements.
double * m_data
Matrix data.
const int & rows(void) const
Return number of matrix rows.
int m_cols
Number of columns.
virtual GMatrixSymmetric operator-(void) const
Negate matrix elements.
virtual GMatrixSymmetric & operator=(const GMatrixSymmetric &matrix)
Matrix assignment operator.
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.
#define G_CHOL_DECOMP
#define G_OP_MAT_MUL
Symmetric matrix class definition.
virtual GVector operator*(const GVector &v) const
Vector multiplication.
GMatrix extract_upper_triangle(void) const
Extract upper triangle of matrix.
void init_members(void)
Initialise class mambers.
int m_elements
Number of elements stored in matrix.
virtual double sum(void) const
Sum matrix elements.
virtual double fill(void) const
Determine fill of matrix.
void copy_members(const GMatrixSymmetric &matrix)
Copy class members.
Sparse matrix class definition.
GMatrixSymmetric abs(void) const
Return absolute of matrix.
GMatrixSymmetric cholesky_decompose(const bool &compress=true) const
Return Cholesky decomposition.
GVector solve(const GVector &vector) const
Solves linear matrix equation.
virtual void add_to_column(const int &column, const GVector &vector)
Add vector column into matrix.
Exception handler interface definition.
virtual GVector row(const int &row) const
Extract row as vector from matrix.
virtual GMatrixBase & operator=(const GMatrixBase &matrix)
Assignment operator.
virtual double & at(const int &row, const int &column)
Return reference to matrix element.
const int & size(void) const
Return size of vector.
Definition: GVector.hpp:180
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_inx
Index array of non-zero rows/columns.
int * m_colstart
Column start indices (m_cols+1)
void alloc_members(const int &rows, const int &columns)
Allocates matrix memory.
int * m_rowsel
Row selection (for compressed decomposition)
#define G_CONSTRUCTOR
const int & columns(void) const
Return number of matrix columns.
#define G_EXTRACT_ROW
#define G_OP_ADD
int * m_colsel
Column selection (for compressed decomposition)
void set_inx(void)
Set index selection.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413