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