GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GMatrixSparse.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GMatrixSparse.cpp - Sparse 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 GMatrixSparse.cpp
23  * @brief Sparse matrix class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GVector.hpp"
35 #include "GMatrix.hpp"
36 #include "GMatrixSparse.hpp"
37 #include "GMatrixSymmetric.hpp"
38 #include "GSparseSymbolic.hpp"
39 #include "GSparseNumeric.hpp"
40 
41 /* __ Method name definitions ____________________________________________ */
42 #define G_CONSTRUCTOR "GMatrixSparse::GMatrixSparse(int&, int&, int&)"
43 #define G_OP_MUL_VEC "GMatrixSparse::operator*(GVector&)"
44 #define G_OP_ADD "GMatrixSparse::operator+=(GMatrixSparse&)"
45 #define G_OP_SUB "GMatrixSparse::operator-=(GMatrixSparse&)"
46 #define G_OP_MAT_MUL "GMatrixSparse::operator*=(GMatrixSparse&)"
47 #define G_AT "GMatrixSparse::at(int&, int&)"
48 #define G_EXTRACT_ROW "GMatrixSparse::row(int&)"
49 #define G_SET_ROW "GMatrixSparse::row(int&, GVector&)"
50 #define G_EXTRACT_COLUMN "GMatrixSparse::column(int&)"
51 #define G_SET_COLUMN "GMatrixSparse::column(int&, GVector&)"
52 #define G_SET_COLUMN2 "GMatrixSparse::column(int&, double*, int*, int)"
53 #define G_ADD_TO_COLUMN "GMatrixSparse::add_to_column(int&, GVector&)"
54 #define G_ADD_TO_COLUMN2 "GMatrixSparse::add_to_column(int&, double*,"\
55  " int*, int)"
56 #define G_CHOL_DECOMP "GMatrixSparse::cholesky_decompose(bool)"
57 #define G_CHOL_SOLVE "GMatrixSparse::cholesky_solver(GVector&, bool)"
58 #define G_STACK_INIT "GMatrixSparse::stack_init(int&, int&)"
59 #define G_STACK_PUSH1 "GMatrixSparse::stack_push_column(GVector&, int&)"
60 #define G_STACK_PUSH2 "GMatrixSparse::stack_push_column(double*, int*, int&,"\
61  " int&)"
62 #define G_STACK_FLUSH "GMatrixSparse::stack_flush(void)"
63 #define G_COPY_MEMBERS "GMatrixSparse::copy_members(GMatrixSparse&)"
64 #define G_ALLOC_MEMBERS "GMatrixSparse::alloc_members(int&, int&, int&)"
65 #define G_GET_INDEX "GMatrixSparse::get_index(int&, int&)"
66 #define G_ALLOC "GMatrixSparse::alloc_elements(int&, int&)"
67 #define G_FREE "GMatrixSparse::free_elements(int&, int&)"
68 #define G_REMOVE_ZERO "GMatrixSparse::remove_zero_row_col(void)"
69 #define G_INSERT_ROW "GMatrixSparse::insert_row(int&, GVector&, bool&)"
70 #define G_SYMPERM "cs_symperm(GMatrixSparse*, int*, int&)"
71 #define G_TRANSPOSE "cs_transpose(GMatrixSparse*, int)"
72 
73 /* __ Macros _____________________________________________________________ */
74 #define G_MIN(a,b) (((a) < (b)) ? (a) : (b))
75 #define G_MAX(a,b) (((a) > (b)) ? (a) : (b))
76 
77 /* __ Coding definitions _________________________________________________ */
78 
79 /* __ Debug definitions __________________________________________________ */
80 //#define G_DEBUG_SPARSE_PENDING // Analyse pending values
81 //#define G_DEBUG_SPARSE_INSERTION // Analyse value insertion
82 //#define G_DEBUG_SPARSE_ADDITION // Analyse value addition
83 //#define G_DEBUG_SPARSE_COMPRESSION // Analyse zero row/column compression
84 //#define G_DEBUG_SPARSE_MALLOC // Analyse memory management
85 //#define G_DEBUG_SPARSE_STACK_PUSH // Analyse stack column pushing
86 //#define G_DEBUG_SPARSE_STACK_FLUSH // Analyse stack column flushing
87 
88 
89 /*==========================================================================
90  = =
91  = Constructors/destructors =
92  = =
93  ==========================================================================*/
94 
95 /***********************************************************************//**
96  * @brief Void matrix constructor
97  *
98  * Constructs empty sparse matrix. The number of rows and columns of the
99  * matrix will be zero.
100  ***************************************************************************/
102 {
103  // Initialise class members
104  init_members();
105 
106  // Return
107  return;
108 }
109 
110 
111 
112 /***********************************************************************//**
113  * @brief Matrix constructor
114  *
115  * @param[in] rows Number of rows [>=0].
116  * @param[in] columns Number of columns [>=0].
117  * @param[in] elements Number of allocated elements.
118  *
119  * @exception GException::invalid_argument
120  * Number of rows or columns is negative.
121  *
122  * Constructs sparse matrix of dimension @p rows times @p columns. The
123  * optional @p elements argument specifies the size of the physical
124  * memory that should be allocated. By default, no memory will be allocated
125  * and memory allocation will be performed on-the-fly. If the amount of
126  * required memory is larger than the size specified by @p elements,
127  * additional momeory will be allocated automatically on-the-fly.
128  ***************************************************************************/
130  const int& columns,
131  const int& elements) :
132  GMatrixBase()
133 {
134  // Compile option: raise an exception if number of rows or columns is
135  // negative
136  #if defined(G_RANGE_CHECK)
137  if (rows < 0) {
138  std::string msg = "Number of rows "+gammalib::str(rows)+" is negative. "
139  "Please specify a non-negative number of rows.";
141  }
142  if (columns < 0) {
143  std::string msg = "Number of columns "+gammalib::str(columns)+" is "
144  "negative. Please specify a non-negative number of "
145  "columns.";
147  }
148  #endif
149 
150  // Initialise private members for clean destruction
151  init_members();
152 
153  // Allocate matrix memory
154  alloc_members(rows, columns, elements);
155 
156  // Return
157  return;
158 }
159 
160 
161 /***********************************************************************//**
162  * @brief GMatrix to GMatrixSparse storage class convertor
163  *
164  * @param[in] matrix Generic matrix (GMatrix).
165  *
166  * Constructs a sparse matrix by using the number of rows and columns of
167  * a generic matrix and by assigning the elements of the generic matrix
168  * to the sparse matrix.
169  ***************************************************************************/
171 {
172  // Initialise class members for clean destruction
173  init_members();
174 
175  // Construct matrix
176  alloc_members(matrix.rows(), matrix.columns());
177 
178  // Fill matrix column by column
179  for (int col = 0; col < m_cols; ++col) {
180  GVector vector = matrix.column(col);
181  this->column(col, vector);
182  }
183 
184  // Return
185  return;
186 }
187 
188 
189 /***********************************************************************//**
190  * @brief GMatrixSymmetric to GMatrixSparse storage class convertor
191  *
192  * @param[in] matrix Symmetric matrix (GMatrixSymmetric).
193  *
194  * Constructs a sparse matrix by using the number of rows and columns of
195  * a symmetric matrix and by assigning the elements of the symmetric matrix
196  * to the sparse matrix.
197  ***************************************************************************/
199 {
200  // Initialise class members for clean destruction
201  init_members();
202 
203  // Allocate matrix memory
204  alloc_members(matrix.rows(), matrix.columns());
205 
206  // Fill matrix column by column
207  for (int col = 0; col < m_cols; ++col) {
208  GVector vector = matrix.column(col);
209  this->column(col, vector);
210  }
211 
212  // Return
213  return;
214 }
215 
216 
217 /***********************************************************************//**
218  * @brief Copy constructor
219  *
220  * @param[in] matrix Matrix.
221  *
222  * Constructs matrix by copying an existing matrix.
223  ***************************************************************************/
225 {
226  // Initialise private members for clean destruction
227  init_members();
228 
229  // Copy members. Note that the copy operation does not fill the
230  // pending element.
231  copy_members(matrix);
232 
233  // Return
234  return;
235 }
236 
237 
238 /***********************************************************************//**
239  * @brief Destructor
240  ***************************************************************************/
242 {
243  // Free members
244  free_members();
245 
246  // Return
247  return;
248 }
249 
250 
251 /*==========================================================================
252  = =
253  = Operators =
254  = =
255  ==========================================================================*/
256 
257 /***********************************************************************//**
258  * @brief Matrix assignment operator
259  *
260  * @param[in] matrix Matrix.
261  * @return Matrix.
262  *
263  * Assigns the content of another matrix to the actual matrix instance.
264  ***************************************************************************/
266 {
267  // Execute only if object is not identical
268  if (this != &matrix) {
269 
270  // Assign base class members. Note that this method will also perform
271  // the allocation of the matrix memory and the copy of the matrix
272  // attributes.
273  this->GMatrixBase::operator=(matrix);
274 
275  // Free derived class members
276  free_members();
277 
278  // Initialise derived class members
279  init_members();
280 
281  // Copy derived class members
282  copy_members(matrix);
283 
284  } // endif: object was not identical
285 
286  // Return
287  return *this;
288 }
289 
290 
291 /***********************************************************************//**
292  * @brief Value assignment operator
293  *
294  * @param[in] value Value.
295  * @return Matrix.
296  *
297  * Assigns the specified @p value to all elements of the matrix.
298  ***************************************************************************/
300 {
301  // Continue only if matrix is not empty
302  if (!is_empty()) {
303 
304  // Fill any pending element to have a non-pending state
305  fill_pending();
306 
307  // If value is 0 then simply reinitialize column start indices
308  if (value == 0) {
309 
310  // Initialise column start indices to 0
311  for (int col = 0; col <= m_cols; ++col) {
312  m_colstart[col] = 0;
313  }
314 
315  }
316 
317  // ... otherwise fill column-wise
318  else {
319 
320  // Set column vector
322  column = value;
323 
324  // Column-wise setting
325  for (int col = 0; col < m_cols; ++col) {
326  this->column(col, column);
327  }
328 
329  }
330 
331  } // endif: matrix was not empty
332 
333  // Return this object
334  return *this;
335 }
336 
337 
338 /***********************************************************************//**
339  * @brief Return reference to matrix element
340  *
341  * @param[in] row Matrix row [0,...,rows()-1].
342  * @param[in] column Matrix column [0,...,columns()-1].
343  * @return Reference to matrix element.
344  *
345  * Returns the reference to the matrix element at @p row and @p column. If
346  * the matrix element does not yet exist, a reference to the pending element
347  * with a value of 0.0 is returned.
348  *
349  * If the matrix row or column are invalid an exception is thrown by the
350  * get_index() method.
351  ***************************************************************************/
352 double& GMatrixSparse::operator()(const int& row, const int& column)
353 {
354  // Fill pending element. This will set the value of the pending element
355  // to 0.0.
356  fill_pending();
357 
358  // Get element
359  int inx = get_index(row, column);
360  double* value;
361  if (inx < 0) {
362  value = &m_fill_val;
363  m_fill_row = row;
364  m_fill_col = column;
365  }
366  else {
367  value = &(m_data[inx]);
368  }
369 
370  // Return element
371  return *value;
372 }
373 
374 
375 /***********************************************************************//**
376  * @brief Return reference to matrix element (const version)
377  *
378  * @param[in] row Matrix row [0,...,rows()-1].
379  * @param[in] column Matrix column [0,...,columns()-1].
380  * @return Const reference to matrix element.
381  *
382  * Returns a const reference to the matrix element at @p row and @p column.
383  * If the matrix element does not yet exist, a reference to the zero element
384  * is returned. If the matrix element corresponds to the pending element,
385  * a reference to the pending element is returned. Otherwise, a reference
386  * to the matrix elements is returned.
387  *
388  * If the matrix row or column are invalid an exception is thrown by the
389  * get_index() method.
390  ***************************************************************************/
391 const double& GMatrixSparse::operator()(const int& row,
392  const int& column) const
393 {
394  // Get element. We need here the zero element to return also a pointer
395  // for 0.0 entry that is not stored. Since we have the const version we
396  // don't have to care about modification of this zero value.
397  int inx = get_index(row, column);
398  double* value;
399  if (inx < 0) {
400  value = (double*)&m_zero;
401  }
402  else if (inx == m_elements) {
403  value = (double*)&m_fill_val;
404  }
405  else {
406  value = (double*)&(m_data[inx]);
407  }
408 
409  // Return element
410  return *value;
411 }
412 
413 
414 /***********************************************************************//**
415  * @brief Vector multiplication
416  *
417  * @param[in] vector Vector.
418  *
419  * @exception GException::invalid_argument
420  * Vector size differs from number of columns in matrix.
421  *
422  * This method performs a vector multiplication of a matrix. The vector
423  * multiplication will produce a vector. The matrix multiplication can only
424  * be performed when the number of matrix columns is equal to the length of
425  * the vector.
426  *
427  * The method includes any pending fill.
428  ***************************************************************************/
430 {
431  // Throw exception if the matrix and vector dimensions are not compatible
432  if (m_cols != vector.size()) {
433  std::string msg = "Vector size "+gammalib::str(vector.size())+" does "
434  "not match "+gammalib::str(m_cols)+" matrix columns. "
435  "Please specify a vector of size "+
436  gammalib::str(m_cols)+".";
438  }
439 
440  // Initialise result vector
441  GVector result(m_rows);
442 
443  // Multiply only if there are elements in matrix
444  if (m_elements > 0) {
445 
446  // Perform vector multiplication
447  for (int col = 0; col < m_cols; ++col) {
448  int i_start = m_colstart[col];
449  int i_stop = m_colstart[col+1];
450  double* ptr_data = m_data + i_start;
451  int* ptr_rowinx = m_rowinx + i_start;
452  for (int i = i_start; i < i_stop; ++i) {
453  result[*ptr_rowinx++] += *ptr_data++ * vector[col];
454  }
455  }
456 
457  // If fill is pending then add-in also this element into the product
458  if (m_fill_val != 0.0) {
459  result[m_fill_row] += m_fill_val * vector[m_fill_col];
460  #if defined(G_DEBUG_SPARSE_PENDING)
461  std::cout << G_OP_MUL_VEC << ": pending value " << m_fill_val
462  << " for location (" << m_fill_row << "," << m_fill_col
463  << ") has been used." << std::endl;
464  #endif
465  }
466 
467  } // endif: there were elements in matrix
468 
469  // Return result
470  return result;
471 }
472 
473 
474 /***********************************************************************//**
475  * @brief Negate matrix elements
476  *
477  * @return Matrix with negated elements.
478  ***************************************************************************/
480 {
481  // Copy matrix
482  GMatrixSparse matrix = *this;
483 
484  // Fill pending element
485  matrix.fill_pending();
486 
487  // Negate all matrix elements
488  double* ptr = matrix.m_data;
489  for (int i = 0; i < matrix.m_elements; ++i, ++ptr) {
490  *ptr = -(*ptr);
491  }
492 
493  // Return matrix
494  return matrix;
495 }
496 
497 
498 /***********************************************************************//**
499  * @brief Equalty operator
500  *
501  * @param[in] matrix Matrix.
502  *
503  * This operator checks if two matrices are identical. Two matrices are
504  * considered identical if they have the same dimensions and identicial
505  * elements.
506  *
507  * @todo Implement native sparse code
508  ***************************************************************************/
509 bool GMatrixSparse::operator==(const GMatrixSparse &matrix) const
510 {
511  // Initalise the result to 'equal matrices'
512  bool result = true;
513 
514  // Perform comparison (only if matrix dimensions are identical)
515  if (m_rows == matrix.m_rows && m_cols == matrix.m_cols) {
516 
517  // Loop over all matrix elements
518  for (int row = 0; row < m_rows; ++row) {
519  for (int col = 0; col < m_cols; ++col) {
520  if ((*this)(row,col) != matrix(row,col)) {
521  result = false;
522  break;
523  }
524  }
525  if (!result) {
526  break;
527  }
528  }
529  }
530  else {
531  result = false;
532  }
533 
534  // Return result
535  return result;
536 }
537 
538 
539 /***********************************************************************//**
540  * @brief Non-equality operator
541  *
542  * @param[in] matrix Matrix.
543  *
544  * This operator checks if two matrices are not identical. Two matrices are
545  * considered not identical if they differ in their dimensions or if at
546  * least one element differs.
547  ***************************************************************************/
548 bool GMatrixSparse::operator!=(const GMatrixSparse &matrix) const
549 {
550  // Get negated result of equality operation
551  bool result = !(this->operator==(matrix));
552 
553  // Return result
554  return result;
555 }
556 
557 
558 /***********************************************************************//**
559  * @brief Unary matrix addition operator
560  *
561  * @param[in] matrix Matrix.
562  *
563  * @exception GException::invalid_argument
564  * Incompatible number of matrix rows or columns.
565  *
566  * This method performs a matrix addition. The operation can only succeed
567  * when the dimensions of both matrices are identical.
568  *
569  * @todo Implement native sparse code
570  ***************************************************************************/
572 {
573  // Throw an exception if the matrix dimensions are not compatible
574  if (m_rows != matrix.m_rows) {
575  std::string msg = "Number of matrix rows "+gammalib::str(matrix.m_rows)+
576  " differs from the expected number of "+
577  gammalib::str(m_rows)+" rows. Please specify a "
578  "matrix with "+gammalib::str(m_rows)+" rows.";
580  }
581  if (m_cols != matrix.m_cols) {
582  std::string msg = "Number of matrix columns "+gammalib::str(matrix.m_cols)+
583  " differs from the expected number of "+
584  gammalib::str(m_cols)+" columns. Please specify a "
585  "matrix with "+gammalib::str(m_cols)+" columns.";
587  }
588 
589  // Perform inplace matrix addition using vectors
590  for (int col = 0; col < m_cols; ++col) {
591  GVector v_result = column(col);
592  GVector v_operand = matrix.column(col);
593  v_result += v_operand;
594  column(col, v_result);
595  }
596 
597  // Return result
598  return *this;
599 }
600 
601 
602 /***********************************************************************//**
603  * @brief Unary matrix scalar addition operator
604  *
605  * @param[in] scalar Scalar.
606  *
607  * Adds a @p scalar to each matrix element.
608  ***************************************************************************/
610 {
611  // Perform inplace matrix addition using vectors
612  for (int col = 0; col < m_cols; ++col) {
613  GVector vector = column(col);
614  vector += scalar;
615  column(col, vector);
616  }
617 
618  // Return result
619  return *this;
620 }
621 
622 
623 /***********************************************************************//**
624  * @brief Unary matrix subtraction operator
625  *
626  * @param[in] matrix Matrix.
627  *
628  * @exception GException::invalid_argument
629  * Incompatible matrix size.
630  *
631  * This method performs a matrix addition. The operation can only succeed
632  * when the dimensions of both matrices are identical.
633  *
634  * @todo Implement native sparse code
635  ***************************************************************************/
637 {
638  // Throw an exception if the matrix dimensions are not compatible
639  if (m_rows != matrix.m_rows) {
640  std::string msg = "Number of matrix rows "+gammalib::str(matrix.m_rows)+
641  " differs from the expected number of "+
642  gammalib::str(m_rows)+" rows. Please specify a "
643  "matrix with "+gammalib::str(m_rows)+" rows.";
645  }
646  if (m_cols != matrix.m_cols) {
647  std::string msg = "Number of matrix columns "+gammalib::str(matrix.m_cols)+
648  " differs from the expected number of "+
649  gammalib::str(m_cols)+" columns. Please specify a "
650  "matrix with "+gammalib::str(m_cols)+" columns.";
652  }
653 
654  // Perform inplace matrix subtraction
655  for (int col = 0; col < m_cols; ++col) {
656  GVector v_result = column(col);
657  GVector v_operand = matrix.column(col);
658  v_result -= v_operand;
659  column(col, v_result);
660  }
661 
662  // Return result
663  return *this;
664 }
665 
666 
667 /***********************************************************************//**
668  * @brief Unary matrix scalar subtraction operator
669  *
670  * @param[in] scalar Scalar.
671  *
672  * Subtracts a @p scalar from each matrix element.
673  ***************************************************************************/
675 {
676  // Perform inplace matrix subtraction
677  for (int col = 0; col < m_cols; ++col) {
678  GVector vector = column(col);
679  vector -= scalar;
680  column(col, vector);
681  }
682 
683  // Return result
684  return *this;
685 }
686 
687 
688 /***********************************************************************//**
689  * @brief Unary matrix multiplication operator
690  *
691  * @param[in] matrix Matrix.
692  *
693  * @exception GException::invalid_argument
694  * Number of rows in @p matrix is different from number of
695  * matrix columns.
696  *
697  * This method performs a matrix multiplication. The operation can only
698  * succeed when the dimensions of both matrices are compatible.
699  *
700  * @todo Implement native sparse code
701  ***************************************************************************/
703 {
704  // Throw exception if the matrix dimensions are not compatible
705  if (m_cols != matrix.m_rows) {
706  std::string msg = "Number of "+gammalib::str(m_cols)+" columns in "
707  "the first matrix differs from number of "+
708  gammalib::str(matrix.m_rows)+" rows in the second "
709  "matrix. Please specify a second matrix with "+
710  gammalib::str(m_cols)+" rows.";
712  }
713 
714  // Allocate result matrix
715  GMatrixSparse result(m_rows, matrix.m_cols);
716 
717  // Multiply only if there are elements in both matrices
718  if (m_elements > 0 && matrix.m_elements > 0) {
719 
720  // Loop over all elements of result matrix
721  for (int row = 0; row < m_rows; ++row) {
722  for (int col = 0; col < matrix.m_cols; ++col) {
723  double sum = 0.0;
724  for (int i = 0; i < m_cols; ++i) {
725  sum += (*this)(row,i) * matrix(i,col);
726  }
727  result(row,col) = sum;
728  }
729  }
730 
731  }
732 
733  // Assign result
734  *this = result;
735 
736  // Return result
737  return *this;
738 }
739 
740 
741 /*==========================================================================
742  = =
743  = Public methods =
744  = =
745  ==========================================================================*/
746 
747 /***********************************************************************//**
748  * @brief Clear matrix
749  ***************************************************************************/
751 {
752  // Free members
753  free_members();
754 
755  // Initialise private members
756  init_members();
757 
758  // Return
759  return;
760 }
761 
762 
763 /***********************************************************************//**
764  * @brief Clone matrix
765  *
766  * @return Pointer to deep copy of matrix.
767  ***************************************************************************/
769 {
770  // Clone matrix
771  return new GMatrixSparse(*this);
772 }
773 
774 
775 /***********************************************************************//**
776  * @brief Return reference to matrix element
777  *
778  * @param[in] row Matrix row [0,...,rows()[.
779  * @param[in] column Matrix column [0,...,columns()[.
780  * @return Reference to matrix element.
781  *
782  * @exception GException::out_of_range
783  * Matrix row or column out of range.
784  ***************************************************************************/
785 double& GMatrixSparse::at(const int& row, const int& column)
786 {
787  // Throw exception if row or column index is out of range
788  if (row < 0 || row >= m_rows) {
789  throw GException::out_of_range(G_AT, "Matrix row", row, m_rows);
790  }
791  if (column < 0 || column >= m_cols) {
792  throw GException::out_of_range(G_AT, "Matrix column", column, m_cols);
793  }
794 
795  // Return element
796  return ((*this)(row, column));
797 }
798 
799 
800 /***********************************************************************//**
801  * @brief Return reference to matrix element (const version)
802  *
803  * @param[in] row Matrix row [0,...,rows()[.
804  * @param[in] column Matrix column [0,...,columns()[.
805  * @return Reference to matrix element.
806  *
807  * @exception GException::out_of_range
808  * Matrix row or column out of range.
809  ***************************************************************************/
810 const double& GMatrixSparse::at(const int& row, const int& column) const
811 {
812  // Throw exception if row or column index is out of range
813  if (row < 0 || row >= m_rows) {
814  throw GException::out_of_range(G_AT, "Matrix row", row, m_rows);
815  }
816  if (column < 0 || column >= m_cols) {
817  throw GException::out_of_range(G_AT, "Matrix column", column, m_cols);
818  }
819 
820  // Return element
821  return ((*this)(row, column));
822 }
823 
824 
825 /***********************************************************************//**
826  * @brief Extract row as vector from matrix
827  *
828  * @param[in] row Matrix row [0,...,rows()[.
829  * @return Vector of matrix row elements (columns() elements).
830  *
831  * @exception GException::out_of_range
832  * Invalid matrix row specified.
833  *
834  * This method extracts a matrix row into a vector.
835  ***************************************************************************/
836 GVector GMatrixSparse::row(const int& row) const
837 {
838  // Throw exception if the row index is invalid
839  #if defined(G_RANGE_CHECK)
840  if (row < 0 || row >= m_rows) {
841  throw GException::out_of_range(G_EXTRACT_ROW, "Matrix row", row, m_rows);
842  }
843  #endif
844 
845  // Create result vector
846  GVector result(m_cols);
847 
848  // Loop over all columns to extract data
849  for (int col = 0; col < m_cols; ++col) {
850 
851  // Get the start and stop of the elements
852  int i_start = m_colstart[col];
853  int i_stop = m_colstart[col+1];
854 
855  // If column is not empty then search requested row by bisection
856  if (i_stop > i_start) {
857 
858  // Find row by bisection
859  while ((i_stop - i_start) > 1) {
860  int mid = (i_start+i_stop) / 2;
861  if (m_rowinx[mid] > row) {
862  i_stop = mid;
863  }
864  else {
865  i_start = mid;
866  }
867  }
868 
869  // Copy element if we found the row
870  if (m_rowinx[i_start] == row) {
871  result[col] = m_data[i_start];
872  }
873 
874  } // endif: column was not empty
875 
876  } // endfor: looped over all columns
877 
878  // If there is a pending element then put it in the vector
879  if (m_fill_val != 0.0 && m_fill_row == row) {
880  result[m_fill_col] = m_fill_val;
881  }
882 
883  // Return vector
884  return result;
885 }
886 
887 
888 /***********************************************************************//**
889  * @brief Set row in matrix
890  *
891  * @param[in] row Matrix row [0,...,rows()[.
892  * @param[in] vector Vector of matrix row elements (columns() elements).
893  ***************************************************************************/
894 void GMatrixSparse::row(const int& row, const GVector& vector)
895 {
896  // Insert row
897  insert_row(row, vector, false);
898 
899  // Return
900  return;
901 }
902 
903 
904 /***********************************************************************//**
905  * @brief Extract column as vector from matrix
906  *
907  * @param[in] column Matrix column [0,...,columns()[.
908  * @return Vector of matrix column elements (rows() elements).
909  *
910  * @exception GException::out_of_range
911  * Invalid matrix column specified.
912  *
913  * This method extracts a matrix column into a vector.
914  ***************************************************************************/
915 GVector GMatrixSparse::column(const int& column) const
916 {
917  // Throw exception if the column index is invalid
918  #if defined(G_RANGE_CHECK)
919  if (column < 0 || column >= m_cols) {
920  throw GException::out_of_range(G_EXTRACT_COLUMN, "Matrix column",
921  column, m_cols);
922  }
923  #endif
924 
925  // Create result vector
926  GVector result(m_rows);
927 
928  // Get the start and stop of the elements
929  int i_start = m_colstart[column];
930  int i_stop = m_colstart[column+1];
931 
932  // Extract elements into vector
933  for (int i = i_start; i < i_stop; ++i) {
934  result[m_rowinx[i]] = m_data[i];
935  }
936 
937  // If there is a pending element then put it in the vector
938  if (m_fill_val != 0.0 && m_fill_col == column) {
939  result[m_fill_row] = m_fill_val;
940  }
941 
942  // Return vector
943  return result;
944 }
945 
946 
947 /***********************************************************************//**
948  * @brief Insert vector column into matrix
949  *
950  * @param[in] column Matrix column [0,...,columns()[.
951  * @param[in] vector Vector.
952  *
953  * @exception GException::out_of_range
954  * Invalid matrix column specified.
955  * @exception GException::invalid_argument
956  * Vector size does not match number of matrix rows.
957  *
958  * Inserts the content of a vector into a matrix column. Any previous
959  * content in the matrix column will be overwritted.
960  *
961  * This is the main driver routine to insert data into a matrix. Note that
962  * there is another instance of this function that takes a compressed array.
963  ***************************************************************************/
964 void GMatrixSparse::column(const int& column, const GVector& vector)
965 {
966  // Debug header
967  #if defined(G_DEBUG_SPARSE_INSERTION)
968  std::cout << "GMatrixSparse::column(";
969  std::cout << column << ", [" << vector << "]):" << std::endl;
970  std::cout << " In Data : ";
971  for (int i = 0; i < m_elements; ++i) {
972  std::cout << m_data[i] << " ";
973  }
974  std::cout << std::endl << " In Row .: ";
975  for (int i = 0; i < m_elements; ++i) {
976  std::cout << m_rowinx[i] << " ";
977  }
978  std::cout << std::endl << " In Col .: ";
979  for (int i = 0; i < m_cols+1; ++i) {
980  std::cout << m_colstart[i] << " ";
981  }
982  std::cout << std::endl;
983  #endif
984 
985  // Throw exception if the column is invalid
986  #if defined(G_RANGE_CHECK)
987  if (column < 0 || column >= m_cols) {
988  throw GException::out_of_range(G_SET_COLUMN, "Matrix column",
989  column, m_cols);
990  }
991  #endif
992 
993  // Throw exception if the vector size does not match the number of rows
994  // in the matrix
995  if (m_rows != vector.size()) {
996  std::string msg = "Vector size "+gammalib::str(vector.size())+
997  " does not match "+gammalib::str(m_rows)+" matrix "
998  "rows. Please specify a vector of size "+
999  gammalib::str(m_rows)+".";
1001  }
1002 
1003  // If there is a pending element for this column then delete it since
1004  // the vector overwrites this element
1005  if (m_fill_val != 0.0 && m_fill_col == column) {
1006  #if defined(G_DEBUG_SPARSE_PENDING)
1007  std::cout << G_SET_COLUMN << ": pending value " << m_fill_val <<
1008  " for location (" << m_fill_row << "," << m_fill_col <<
1009  ") became obsolete" << std::endl;
1010  #endif
1011  m_fill_val = 0.0;
1012  m_fill_row = 0;
1013  m_fill_col = 0;
1014  }
1015 
1016  // Determine the number of non-zero elements in the vector
1017  int n_vector = 0;
1018  for (int row = 0; row < m_rows; ++row) {
1019  if (vector[row] != 0.0) {
1020  n_vector++;
1021  }
1022  }
1023 
1024  // Get the start and stop indices of the actual column and compute
1025  // the number of exisiting elements in the column
1026  int i_start = m_colstart[column];
1027  int i_stop = m_colstart[column+1];
1028  int n_exist = i_stop - i_start;
1029 
1030  // Compute the size difference for the new matrix. It is positive if
1031  // the number of non-zero entries in the vector is larger than the
1032  // number of non-zero entries in the matrix (in this case we have to
1033  // increase the matrix size).
1034  int n_diff = n_vector - n_exist;
1035 
1036  // If we need space then allocate it, if we have to much space then free it
1037  if (n_diff > 0) {
1038  alloc_elements(i_start, n_diff);
1039  #if defined(G_DEBUG_SPARSE_INSERTION)
1040  std::cout << " Insert .: " << n_diff << " elements at index " << i_start << std::endl;
1041  #endif
1042  }
1043  else if (n_diff < 0) {
1044  free_elements(i_start, -n_diff);
1045  #if defined(G_DEBUG_SPARSE_INSERTION)
1046  std::cout << " Remove .: " << -n_diff << " elements at index " << i_start << std::endl;
1047  #endif
1048  }
1049 
1050  // Insert the vector elements in the matrix
1051  if (n_vector > 0) {
1052  for (int row = 0, i = i_start; row < m_rows; ++row) {
1053  if (vector[row] != 0.0) {
1054  m_data[i] = vector[row];
1055  m_rowinx[i] = row;
1056  i++;
1057  }
1058  }
1059  }
1060 
1061  // Update column start indices
1062  for (int i = column+1; i <= m_cols; ++i) {
1063  m_colstart[i] += n_diff;
1064  }
1065 
1066  // Debugging: show sparse matrix after insertion
1067  #if defined(G_DEBUG_SPARSE_INSERTION)
1068  std::cout << " Out Data: ";
1069  for (int i = 0; i < m_elements; ++i) {
1070  std::cout << m_data[i] << " ";
1071  }
1072  std::cout << std::endl << " Out Row : ";
1073  for (int i = 0; i < m_elements; ++i) {
1074  std::cout << m_rowinx[i] << " ";
1075  }
1076  std::cout << std::endl << " Out Col : ";
1077  for (int i = 0; i < m_cols+1; ++i) {
1078  std::cout << m_colstart[i] << " ";
1079  }
1080  std::cout << std::endl;
1081  #endif
1082 
1083  // Return
1084  return;
1085 }
1086 
1087 
1088 /***********************************************************************//**
1089  * @brief Insert compressed array into matrix column
1090  *
1091  * @param[in] column Matrix column [0,...,columns()[.
1092  * @param[in] values Compressed array.
1093  * @param[in] rows Row numbers of array.
1094  * @param[in] number Number of elements in array.
1095  *
1096  * @exception GException::out_of_range
1097  * Invalid matrix column specified.
1098  * @exception GException::invalid_argument
1099  * Row numbers are not smaller than number of matrix rows
1100  *
1101  * Inserts the content of a copressed array into a matrix column. Any
1102  * previous content in the matrix column will be overwritted.
1103  *
1104  * This is the main driver routine to insert data into a matrix. Note that
1105  * there is another instance of this function that takes a vector.
1106  ***************************************************************************/
1107 void GMatrixSparse::column(const int& column, const double* values,
1108  const int* rows, int number)
1109 {
1110  // Debug header
1111  #if defined(G_DEBUG_SPARSE_INSERTION)
1112  std::cout << "GMatrixSparse::column(";
1113  std::cout << column << ", values, rows, " << number << "):" << std::endl;
1114  std::cout << " Matrix Data : ";
1115  for (int i = 0; i < m_elements; ++i) {
1116  std::cout << m_data[i] << " ";
1117  }
1118  std::cout << std::endl << " Matrix Row .: ";
1119  for (int i = 0; i < m_elements; ++i) {
1120  std::cout << m_rowinx[i] << " ";
1121  }
1122  std::cout << std::endl << " Matrix Col .: ";
1123  for (int i = 0; i < m_cols+1; ++i) {
1124  std::cout << m_colstart[i] << " ";
1125  }
1126  std::cout << std::endl;
1127  std::cout << " Array Data .: ";
1128  for (int i = 0; i < number; ++i) {
1129  std::cout << values[i] << " ";
1130  }
1131  std::cout << std::endl << " Array Row ..: ";
1132  for (int i = 0; i < number; ++i) {
1133  std::cout << rows[i] << " ";
1134  }
1135  std::cout << std::endl;
1136  #endif
1137 
1138  // Throw exception if the column is invalid
1139  #if defined(G_RANGE_CHECK)
1140  if (column < 0 || column >= m_cols) {
1141  throw GException::out_of_range(G_SET_COLUMN2, "Matrix column",
1142  column, m_cols);
1143  }
1144  #endif
1145 
1146  // Throw exception if the index array seems incompatible with matrix
1147  // dimensions
1148  if (rows[number-1] >= m_rows) {
1149  std::string msg = "Row number "+gammalib::str(rows[number-1])+" for "
1150  "element "+gammalib::str(number)+" is not smaller "
1151  "than number "+gammalib::str(m_rows)+" of matrix "
1152  "rows. Please specify only row numbers smaller than "+
1153  gammalib::str(m_rows)+".";
1155  }
1156 
1157  // If there is a pending element for this column then delete it since
1158  // the vector overwrites this element
1159  if (m_fill_val != 0.0 && m_fill_col == column) {
1160  #if defined(G_DEBUG_SPARSE_PENDING)
1161  std::cout << G_INSERT_COL2 << ": pending value " << m_fill_val <<
1162  " for location (" << m_fill_row << "," << m_fill_col <<
1163  ") became obsolete" << std::endl;
1164  #endif
1165  m_fill_val = 0.0;
1166  m_fill_row = 0;
1167  m_fill_col = 0;
1168  }
1169 
1170  // Get the start and stop indices of the actual column and compute
1171  // the number of exisiting elements in the column
1172  int i_start = m_colstart[column];
1173  int i_stop = m_colstart[column+1];
1174  int n_exist = i_stop - i_start;
1175 
1176  // If the array is empty then make sure that the number of elements is 0
1177  // (we then just delete the existing column)
1178  if (!values || !rows) {
1179  number = 0;
1180  }
1181 
1182  // Compute the size difference for the new matrix. It is positive if
1183  // the number of non-zero entries in the array is larger than the
1184  // number of non-zero entries in the matrix (in this case we have to
1185  // increase the matrix size).
1186  int n_diff = number - n_exist;
1187 
1188  // If we need space then allocate it, if we have to much space then free it
1189  if (n_diff > 0) {
1190  alloc_elements(i_start, n_diff);
1191  #if defined(G_DEBUG_SPARSE_INSERTION)
1192  std::cout << " Insert .: " << n_diff << " elements at index " << i_start << std::endl;
1193  #endif
1194  }
1195  else if (n_diff < 0) {
1196  free_elements(i_start, -n_diff);
1197  #if defined(G_DEBUG_SPARSE_INSERTION)
1198  std::cout << " Remove .: " << -n_diff << " elements at index " << i_start << std::endl;
1199  #endif
1200  }
1201 
1202  // Insert the array elements into the matrix
1203  if (number > 0) {
1204  for (int row = 0, i = i_start; row < number; ++row, ++i) {
1205  m_data[i] = values[row];
1206  m_rowinx[i] = rows[row];
1207  }
1208  }
1209 
1210  // Update column start indices
1211  for (int i = column+1; i <= m_cols; ++i) {
1212  m_colstart[i] += n_diff;
1213  }
1214 
1215  // Debugging: show sparse matrix after insertion
1216  #if defined(G_DEBUG_SPARSE_INSERTION)
1217  std::cout << " Out Data: ";
1218  for (int i = 0; i < m_elements; ++i) {
1219  std::cout << m_data[i] << " ";
1220  }
1221  std::cout << std::endl << " Out Row : ";
1222  for (int i = 0; i < m_elements; ++i) {
1223  std::cout << m_rowinx[i] << " ";
1224  }
1225  std::cout << std::endl << " Out Col : ";
1226  for (int i = 0; i < m_cols+1; ++i) {
1227  std::cout << m_colstart[i] << " ";
1228  }
1229  std::cout << std::endl;
1230  #endif
1231 
1232  // Return
1233  return;
1234 }
1235 
1236 
1237 /***********************************************************************//**
1238  * @brief Add row to matrix elements
1239  *
1240  * @param[in] row Matrix row [0,...,row()[.
1241  * @param[in] vector Vector.
1242  ***************************************************************************/
1243 void GMatrixSparse::add_to_row(const int& row, const GVector& vector)
1244 {
1245  // Add row
1246  insert_row(row, vector, true);
1247 
1248  // Return
1249  return;
1250 }
1251 
1252 
1253 /***********************************************************************//**
1254  * @brief Add vector column into matrix
1255  *
1256  * @param[in] column Matrix column [0,...,columns()[.
1257  * @param[in] vector Vector.
1258  *
1259  * @exception GException::out_of_range
1260  * Invalid matrix column specified.
1261  * @exception GException::invalid_argument
1262  * Matrix dimension mismatches the vector size.
1263  *
1264  * Adds the contect of a vector to a matrix column.
1265  *
1266  * This is the main driver routine to add data to a matrix. It handles both
1267  * normal and stack-based filled. Note that there is another instance of this
1268  * method that takes a compressed array.
1269  ***************************************************************************/
1270 void GMatrixSparse::add_to_column(const int& column, const GVector& vector)
1271 {
1272  // Debug header
1273  #if defined(G_DEBUG_SPARSE_ADDITION)
1274  std::cout << "GMatrixSparse::add_col(";
1275  std::cout << column << ", [" << v << "]):" << std::endl;
1276  std::cout << " In Data : ";
1277  for (int i = 0; i < m_elements; ++i) {
1278  std::cout << m_data[i] << " ";
1279  }
1280  std::cout << std::endl << " In Row .: ";
1281  for (int i = 0; i < m_elements; ++i) {
1282  std::cout << m_rowinx[i] << " ";
1283  }
1284  std::cout << std::endl << " In Col .: ";
1285  for (int i = 0; i < m_cols+1; ++i) {
1286  std::cout << m_colstart[i] << " ";
1287  }
1288  std::cout << std::endl;
1289  #endif
1290 
1291  // Initialise number of non-zero elements to 0
1292  int non_zero = 0;
1293 
1294  // If we have a stack then try to push vector on stack first. Note that
1295  // stack_push_column does its own argument verifications, so to avoid
1296  // double checking we don't do anything before this call ...
1297  if (m_stack_data != NULL) {
1298  non_zero = stack_push_column(vector, column);
1299  if (non_zero == 0) {
1300  return;
1301  }
1302  }
1303 
1304  // ... otherwise check first the arguments and determine the number of
1305  // non-zero elements in vector
1306  else {
1307 
1308  // Throw exception if the column is invalid
1309  #if defined(G_RANGE_CHECK)
1310  if (column < 0 || column >= m_cols) {
1311  throw GException::out_of_range(G_ADD_TO_COLUMN, "Matrix column",
1312  column, m_cols);
1313  }
1314  #endif
1315 
1316  // Throw exception if the vector size does not match the number of rows
1317  // in the matrix
1318  if (m_rows != vector.size()) {
1319  std::string msg = "Vector size "+gammalib::str(vector.size())+
1320  " does not match "+gammalib::str(m_rows)+" matrix "
1321  "rows. Please specify a vector of size "+
1322  gammalib::str(m_rows)+".";
1324  }
1325 
1326  // Determine number of non-zero elements in vector
1327  non_zero = vector.non_zeros();
1328  }
1329 
1330  // Extract vector for column, add elements, and re-insert vector (only if
1331  // vector to insert has non-zeros)
1332  if (non_zero > 0) {
1333 
1334  // Copy input vector
1335  GVector v_column = vector;
1336 
1337  // Add elements to vector
1338  for (int i = m_colstart[column]; i < m_colstart[column+1]; ++i) {
1339  v_column[m_rowinx[i]] += m_data[i];
1340  }
1341 
1342  // If there is a pending element then put it in the vector
1343  if (m_fill_val != 0.0 && m_fill_col == column) {
1344  v_column[m_fill_row] += m_fill_val;
1345  }
1346 
1347  // Set vector into matrix
1348  this->column(column, v_column);
1349 
1350  }
1351 
1352  // Debugging: show sparse matrix after addition
1353  #if defined(G_DEBUG_SPARSE_ADDITION)
1354  std::cout << " Out Data: ";
1355  for (int i = 0; i < m_elements; ++i) {
1356  std::cout << m_data[i] << " ";
1357  }
1358  std::cout << std::endl << " Out Row : ";
1359  for (int i = 0; i < m_elements; ++i) {
1360  std::cout << m_rowinx[i] << " ";
1361  }
1362  std::cout << std::endl << " Out Col : ";
1363  for (int i = 0; i < m_cols+1; ++i) {
1364  std::cout << m_colstart[i] << " ";
1365  }
1366  std::cout << std::endl;
1367  #endif
1368 
1369  // Return
1370  return;
1371 }
1372 
1373 
1374 /***********************************************************************//**
1375  * @brief Add compressed array into matrix column
1376  *
1377  * @param[in] column Matrix column [0,...,columns()[.
1378  * @param[in] values Compressed array.
1379  * @param[in] rows Row indices of array.
1380  * @param[in] number Number of elements in array.
1381  *
1382  * @exception GException::out_of_range
1383  * Invalid matrix column specified.
1384  * @exception GException::invalid_argument
1385  * Matrix dimension mismatches the vector size.
1386  *
1387  * Adds the content of a compressed array into a matrix column.
1388  *
1389  * This is the main driver routine to add data to a matrix. It handles both
1390  * normal and stack-based filled. Note that there is another instance of this
1391  * method that takes a vector.
1392  ***************************************************************************/
1393 void GMatrixSparse::add_to_column(const int& column, const double* values,
1394  const int* rows, int number)
1395 {
1396  // Debug header
1397  #if defined(G_DEBUG_SPARSE_ADDITION)
1398  std::cout << "GMatrixSparse::add_col(";
1399  std::cout << column << ", values, rows, " << number << "):" << std::endl;
1400  std::cout << " Matrix Data : ";
1401  for (int i = 0; i < m_elements; ++i) {
1402  std::cout << m_data[i] << " ";
1403  }
1404  std::cout << std::endl << " Matrix Row .: ";
1405  for (int i = 0; i < m_elements; ++i) {
1406  std::cout << m_rowinx[i] << " ";
1407  }
1408  std::cout << std::endl << " Matrix Col .: ";
1409  for (int i = 0; i < m_cols+1; ++i) {
1410  std::cout << m_colstart[i] << " ";
1411  }
1412  std::cout << std::endl;
1413  std::cout << " Array Data .: ";
1414  for (int i = 0; i < number; ++i) {
1415  std::cout << values[i] << " ";
1416  }
1417  std::cout << std::endl << " Array Row ..: ";
1418  for (int i = 0; i < number; ++i) {
1419  std::cout << rows[i] << " ";
1420  }
1421  std::cout << std::endl;
1422  #endif
1423 
1424  // If we have a stack then try to push elements on stack first. Note that
1425  // stack_push_column does its own argument verifications, so to avoid
1426  // double checking we don't do anything before this call ...
1427  if (m_stack_data != NULL) {
1428  number = stack_push_column(values, rows, number, column);
1429  if (number == 0) {
1430  return;
1431  }
1432  }
1433 
1434  // ... otherwise check the arguments
1435  else {
1436 
1437  // If the array is empty there is nothing to do
1438  if (!values || !rows || (number < 1)) {
1439  return;
1440  }
1441 
1442  // Throw exception if the column is invalid
1443  #if defined(G_RANGE_CHECK)
1444  if (column < 0 || column >= m_cols) {
1445  throw GException::out_of_range(G_ADD_TO_COLUMN2, "Matrix column",
1446  column, m_cols);
1447  }
1448  #endif
1449 
1450  // Throw exception if the index array seems incompatible with matrix
1451  // dimensions
1452  if (rows[number-1] >= m_rows) {
1453  std::string msg = "Row number "+gammalib::str(rows[number-1])+" for "
1454  "element "+gammalib::str(number)+" is not smaller "
1455  "than number "+gammalib::str(m_rows)+" of matrix "
1456  "rows. Please specify only row numbers smaller than "+
1457  gammalib::str(m_rows)+".";
1459  }
1460 
1461  } // endelse: there was no stack
1462 
1463  // Get indices of column in matrix
1464  int i_start = m_colstart[column];
1465  int i_stop = m_colstart[column+1];
1466 
1467  // Case A: the column exists in the matrix, so mix new elements with existing
1468  // data
1469  if (i_start < i_stop) {
1470 
1471  // Fill pending element before the merge (it will be lost otherwise)
1472  fill_pending();
1473 
1474  // Allocate workspace to hold combined column
1475  int wrk_size = number + i_stop - i_start;
1476  double* wrk_double = new double[wrk_size];
1477  int* wrk_int = new int[wrk_size];
1478 
1479  // Mix matrix column with specified data
1480  int num_mix;
1481  mix_column(&(m_data[i_start]), &(m_rowinx[i_start]), i_stop-i_start,
1482  values, rows, number,
1483  wrk_double, wrk_int, &num_mix);
1484 
1485  // Insert mixed column
1486  this->column(column, wrk_double, wrk_int, num_mix);
1487 
1488  // Free workspace
1489  delete [] wrk_int;
1490  delete [] wrk_double;
1491 
1492  } // endif: Case A
1493 
1494  // Case B: the column does not yet exist in the matrix, so just insert it
1495  else {
1496  this->column(column, values, rows, number);
1497  }
1498 
1499  // Debugging: show sparse matrix after insertion
1500  #if defined(G_DEBUG_SPARSE_ADDITION)
1501  std::cout << " Out Data: ";
1502  for (int i = 0; i < m_elements; ++i) {
1503  std::cout << m_data[i] << " ";
1504  }
1505  std::cout << std::endl << " Out Row : ";
1506  for (int i = 0; i < m_elements; ++i) {
1507  std::cout << m_rowinx[i] << " ";
1508  }
1509  std::cout << std::endl << " Out Col : ";
1510  for (int i = 0; i < m_cols+1; ++i) {
1511  std::cout << m_colstart[i] << " ";
1512  }
1513  std::cout << std::endl;
1514  #endif
1515 
1516  // Return
1517  return;
1518 }
1519 
1520 
1521 /***********************************************************************//**
1522  * @brief Return transposed matrix
1523  *
1524  * @return Transposed matrix.
1525  *
1526  * Returns transposed matrix of the matrix.
1527  ***************************************************************************/
1529 {
1530  // Copy matrix
1531  GMatrixSparse matrix = *this;
1532 
1533  // Fill pending element
1534  matrix.fill_pending();
1535 
1536  // Compute the transpose
1537  matrix = cs_transpose(matrix, 1);
1538 
1539  // Return matrix
1540  return matrix;
1541 }
1542 
1543 
1544 /***********************************************************************//**
1545  * @brief Return inverted matrix
1546  *
1547  * @return Inverted matrix.
1548  *
1549  * Returns inverse of matrix. Inversion is done for the moment using Cholesky
1550  * decomposition. This does not work on any kind of matrix.
1551  *
1552  * @todo Specify in documentation for which kind of matrix the method works.
1553  ***************************************************************************/
1555 {
1556  // Copy matrix
1557  GMatrixSparse matrix = *this;
1558 
1559  // Fill pending element
1560  matrix.fill_pending();
1561 
1562  // Invert matrix
1563  matrix = matrix.cholesky_invert(true);
1564 
1565  // Return matrix
1566  return matrix;
1567 }
1568 
1569 
1570 /***********************************************************************//**
1571  * @brief Solves linear matrix equation
1572  *
1573  * @param[in] vector Solution vector.
1574  * @return Solutions of matrix equation.
1575  *
1576  * Solves the linear equation
1577  *
1578  * \f[M \times {\tt solution} = {\tt vector} \f]
1579  *
1580  * where \f$M\f$ is the matrix, \f${\tt vector}\f$ is the result, and
1581  * \f${\tt solution}\f$ is the solution. Solving is done using Cholesky
1582  * decomposition. This does not work on any kind of matrix.
1583  *
1584  * If the matrix is empty and the vector has a zero length the method
1585  * returns an empty vector.
1586  *
1587  * @todo Specify in documentation for which kind of matrix the method works.
1588  ***************************************************************************/
1590 {
1591  // Initialise result with an empty vector
1592  GVector result;
1593 
1594  // Continue only if matrix is not empty or vector size is not zero
1595  if (!is_empty() || vector.size() != 0) {
1596 
1597  // Get Cholesky decomposition of matrix
1598  GMatrixSparse decomposition = cholesky_decompose(true);
1599 
1600  // Solve linear equation
1601  result = decomposition.cholesky_solver(vector);
1602 
1603  }
1604 
1605  // Return result
1606  return result;
1607 }
1608 
1609 
1610 /***********************************************************************//**
1611  * @brief Return absolute of matrix
1612  *
1613  * @return Absolute of matrix
1614  *
1615  * Returns matrix where all elements of the matrix have been replaced by
1616  * their absolute values.
1617  ***************************************************************************/
1619 {
1620  // Copy matrix
1621  GMatrixSparse matrix = *this;
1622 
1623  // Fill pending element
1624  matrix.fill_pending();
1625 
1626  // Take the absolute value of all matrix elements
1627  double* ptr = matrix.m_data;
1628  for (int i = 0; i < matrix.m_elements; ++i, ++ptr) {
1629  *ptr = std::abs(*ptr);
1630  }
1631 
1632  // Return matrix
1633  return matrix;
1634 }
1635 
1636 
1637 /***********************************************************************//**
1638  * @brief Returns fill of matrix
1639  *
1640  * The fill of a matrix is defined as the number non-zero elements devided
1641  * by the number of total elements. By definition, the fill is comprised
1642  * in the interval [0,..,1]. The fill of an undefined matrix is defined to
1643  * be 0.
1644  ***************************************************************************/
1645 double GMatrixSparse::fill(void) const
1646 {
1647  // Initialise result
1648  double result = 0.0;
1649 
1650  // Compute matrix size
1651  int size = m_rows*m_cols;
1652 
1653  // Continue only if matrix has elements
1654  if (size > 0) {
1655 
1656  // Determine number of elements in matrix
1657  int num = (m_fill_val == 0.0) ? m_elements : m_elements + 1;
1658 
1659  // Compute fill
1660  result = double(num) / double(size);
1661 
1662  } // endif: there were elements in matrix
1663 
1664  // Return fill
1665  return result;
1666 }
1667 
1668 
1669 /***********************************************************************//**
1670  * @brief Return minimum matrix element
1671  ***************************************************************************/
1672 double GMatrixSparse::min(void) const
1673 {
1674  // Initialise minimum with fill value
1675  double result = m_fill_val;
1676 
1677  // Search all elements for the smallest one
1678  for (int i = 0; i < m_elements; ++i) {
1679  if (m_data[i] < result) {
1680  result = m_data[i];
1681  }
1682  }
1683 
1684  // If minimum is > 0.0 and there are zero elements then set the minimum
1685  // to 0.0
1686  if ((result > 0.0) && (m_elements < (m_rows*m_cols))) {
1687  result = 0.0;
1688  }
1689 
1690  // Return result
1691  return result;
1692 }
1693 
1694 
1695 /***********************************************************************//**
1696  * @brief Return maximum matrix element
1697  ***************************************************************************/
1698 double GMatrixSparse::max(void) const
1699 {
1700  // Initialise maximum with fill value
1701  double result = m_fill_val;
1702 
1703  // Search all elements for the largest one
1704  for (int i = 0; i < m_elements; ++i) {
1705  if (m_data[i] > result) {
1706  result = m_data[i];
1707  }
1708  }
1709 
1710  // If maximum is < 0.0 and there are zero elements then set the maximum
1711  // to 0.0
1712  if ((result < 0.0) && (m_elements < (m_rows*m_cols))) {
1713  result = 0.0;
1714  }
1715 
1716  // Return result
1717  return result;
1718 }
1719 
1720 
1721 /***********************************************************************//**
1722  * @brief Sum matrix elements
1723  ***************************************************************************/
1724 double GMatrixSparse::sum(void) const
1725 {
1726  // Initialise matrix sum with fill value
1727  double result = m_fill_val;
1728 
1729  // Add all elements
1730  for (int i = 0; i < m_elements; ++i) {
1731  result += m_data[i];
1732  }
1733 
1734  // Return result
1735  return result;
1736 }
1737 
1738 
1739 /***********************************************************************//**
1740  * @brief Return Cholesky decomposition
1741  *
1742  * @param[in] compress Use zero-row/column compression (defaults to true).
1743  * @return Cholesky decomposition of matrix
1744  *
1745  * Returns the Cholesky decomposition of a sparse matrix. The decomposition
1746  * is stored within a GMatrixSparse object.
1747  ***************************************************************************/
1749 {
1750  // Create copy of matrix
1751  GMatrixSparse matrix = *this;
1752 
1753  // Continue only if matrix is not empty
1754  if (matrix.m_rows > 0 && matrix.m_cols > 0) {
1755 
1756  // Save original matrix size
1757  int matrix_rows = matrix.m_rows;
1758  int matrix_cols = matrix.m_cols;
1759 
1760  // Delete any existing symbolic and numeric analysis object and reset
1761  // pointers
1762  if (matrix.m_symbolic != NULL) delete matrix.m_symbolic;
1763  if (matrix.m_numeric != NULL) delete matrix.m_numeric;
1764  matrix.m_symbolic = NULL;
1765  matrix.m_numeric = NULL;
1766 
1767  // Allocate symbolic analysis object
1768  GSparseSymbolic* symbolic = new GSparseSymbolic();
1769 
1770  // Declare numeric analysis object. We don't allocate one since we'll
1771  // throw it away at the end of the function (the L matrix will be copied
1772  // in this object)
1773  GSparseNumeric numeric;
1774 
1775  // Fill pending element into matrix
1776  matrix.fill_pending();
1777 
1778  // Remove rows and columns containing only zeros if matrix compression
1779  // has been selected
1780  if (compress) {
1781  matrix.remove_zero_row_col();
1782  }
1783 
1784  // Ordering an symbolic analysis of matrix. This sets up an array 'pinv'
1785  // which contains the fill-in reducing permutations
1786  symbolic->cholesky_symbolic_analysis(1, matrix);
1787 
1788  // Store symbolic pointer in sparse matrix object
1789  matrix.m_symbolic = symbolic;
1790 
1791  // Perform numeric Cholesky decomposition
1792  numeric.cholesky_numeric_analysis(matrix, *symbolic);
1793 
1794  // Copy L matrix into this object
1795  matrix.free_elements(0, matrix.m_elements);
1796  matrix.alloc_elements(0, numeric.m_L->m_elements);
1797  for (int i = 0; i < matrix.m_elements; ++i) {
1798  matrix.m_data[i] = numeric.m_L->m_data[i];
1799  matrix.m_rowinx[i] = numeric.m_L->m_rowinx[i];
1800  }
1801  for (int col = 0; col <= matrix.m_cols; ++col) {
1802  matrix.m_colstart[col] = numeric.m_L->m_colstart[col];
1803  }
1804 
1805  // Insert zero rows and columns if they have been removed previously.
1806  if (compress) {
1807  matrix.insert_zero_row_col(matrix_rows, matrix_cols);
1808  }
1809 
1810  } // endif: matrix was not empty
1811 
1812  // Return matrix
1813  return matrix;
1814 }
1815 
1816 
1817 /***********************************************************************//**
1818  * @brief Cholesky solver
1819  *
1820  * @param[in] vector Solution vector.
1821  * @param[in] compress Request matrix compression.
1822  *
1823  * @exception GException::invalid_argument
1824  * Matrix and vector do not match.
1825  * @exception GException::invalid_value
1826  * Matrix has not been factorised.
1827  *
1828  * Solves the linear equation A*x=b using a Cholesky decomposition of A.
1829  * This function is to be applied on a GMatrixSparse matrix for which a
1830  * Choleksy factorization has been produced using 'cholesky_decompose'.
1831  ***************************************************************************/
1833  const bool& compress) const
1834 {
1835  // Dump header
1836  #if defined(G_DEBUG_SPARSE_COMPRESSION)
1837  std::cout << "GMatrixSparse::cholesky_solver" << std::endl;
1838  std::cout << " Input vector .....: " << vector << std::endl;
1839  #endif
1840 
1841  // Throw exception if the matrix and vector dimensions are incompatible
1842  if (m_rows != vector.size()) {
1843  std::string msg = "Vector size "+gammalib::str(vector.size())+
1844  " does not match "+gammalib::str(m_rows)+" matrix "
1845  "rows. Please specify a vector of size "+
1846  gammalib::str(m_rows)+".";
1848  }
1849 
1850  // Throw exception if there is no symbolic pointer
1851  if (!m_symbolic) {
1852  std::string msg = "Matrix not factorised. Please call method "
1853  "GMatrixSparse::cholesky_decompose() before calling "
1854  "this method.";
1856  }
1857 
1858  // Throw exception if there is no permutation
1859  if (!m_symbolic->m_pinv) {
1860  std::string msg = "Matrix not factorised. Please call method "
1861  "GMatrixSparse::cholesky_decompose() before calling "
1862  "this method.";
1864  }
1865 
1866  // Flag row and column compression
1867  int row_compressed = (m_rowsel != NULL && m_num_rowsel < m_rows);
1868  int col_compressed = (m_colsel != NULL && m_num_colsel < m_cols);
1869 
1870  // Decide if we need a compression algorithm or not
1871  int no_zero = !(compress && (row_compressed || col_compressed));
1872 
1873  // Allocate vector for permutation and result vector
1874  GVector result(m_cols);
1875 
1876  // Continue only if matrix is not empty
1877  if (m_rows > 0 && m_cols > 0) {
1878 
1879  // Setup pointers to L matrix and x vector
1880  int* Lp = m_colstart;
1881  int* Li = m_rowinx;
1882  double* Lx = m_data;
1883 
1884  // Case A: no zero-row/col compression needed
1885  if (no_zero) {
1886 
1887  // Declare working vector
1888  GVector x(m_rows);
1889 
1890  // Perform inverse vector permutation
1891  for (int i = 0; i < vector.size(); ++i) {
1892  x[m_symbolic->m_pinv[i]] = vector[i];
1893  }
1894 
1895  // Inplace solve L\x=x
1896  for (int col = 0; col < m_cols; ++col) { // loop over columns
1897  x[col] /= Lx[Lp[col]]; // divide by diag.
1898  for (int p = Lp[col]+1; p < Lp[col+1]; p++) { // loop over elements
1899  x[Li[p]] -= Lx[p] * x[col];
1900  }
1901  }
1902 
1903  // Inplace solve L'\x=x
1904  for (int col = m_cols-1; col >= 0; --col) { // loop over columns
1905  for (int p = Lp[col]+1; p < Lp[col+1]; p++) // loop over elements
1906  x[col] -= Lx[p] * x[Li[p]];
1907  x[col] /= Lx[Lp[col]];
1908  }
1909 
1910  // Perform vector permutation
1911  for (int i = 0; i < m_cols; ++i) {
1912  result[i] = x[m_symbolic->m_pinv[i]];
1913  }
1914 
1915  } // endif: Case A
1916 
1917  // Case B: zero-row/column compression requested
1918  else {
1919 
1920  // Allocate row and column mapping arrays
1921  int* row_map = new int[m_rows];
1922  int* col_map = new int[m_cols];
1923 
1924  // Setup row mapping array that maps original matrix rows into compressed
1925  // matrix rows. An entry of -1 indicates that the row should be dropped.
1926  // If no selection exists then setup an identity map.
1927  if (row_compressed) {
1928  for (int row = 0; row < m_rows; ++row) {
1929  row_map[row] = -1;
1930  }
1931  for (int c_row = 0; c_row < m_num_rowsel; ++c_row) {
1932  row_map[m_rowsel[c_row]] = c_row;
1933  }
1934  }
1935  else {
1936  for (int row = 0; row < m_rows; ++row) {
1937  row_map[row] = row;
1938  }
1939  }
1940  #if defined(G_DEBUG_SPARSE_COMPRESSION)
1941  std::cout << " Row mapping ......:";
1942  for (int row = 0; row < m_rows; ++row) {
1943  std::cout << " " << row_map[row];
1944  }
1945  std::cout << std::endl;
1946  #endif
1947 
1948  // Setup column mapping array that maps original matrix column into
1949  // compressed matrix columns. An entry of -1 indicates that the
1950  // column should be dropped. If no selection exists then setup an
1951  // identity map.
1952  if (col_compressed) {
1953  for (int col = 0; col < m_cols; ++col) {
1954  col_map[col] = -1;
1955  }
1956  for (int c_col = 0; c_col < m_num_colsel; ++c_col) {
1957  col_map[m_colsel[c_col]] = c_col;
1958  }
1959  }
1960  else {
1961  for (int col = 0; col < m_cols; ++col) {
1962  col_map[col] = col;
1963  }
1964  }
1965  #if defined(G_DEBUG_SPARSE_COMPRESSION)
1966  std::cout << " Column mapping ...:";
1967  for (int col = 0; col < m_cols; ++col) {
1968  std::cout << " " << col_map[col];
1969  }
1970  std::cout << std::endl;
1971  #endif
1972 
1973  // Declare working vector
1974  GVector x(row_compressed ? m_num_rowsel : m_rows);
1975 
1976  // Compress input vector v -> c_v if required
1977  if (m_rowsel != NULL && m_num_rowsel < m_rows) {
1978  for (int c_row = 0; c_row < m_num_rowsel; ++c_row) {
1979  x[c_row] = vector[m_rowsel[c_row]];
1980  }
1981  }
1982  else {
1983  x = vector;
1984  }
1985  #if defined(G_DEBUG_SPARSE_COMPRESSION)
1986  std::cout << " Compressed vector : " << x << std::endl;
1987  #endif
1988 
1989  // Perform inverse permutation
1990  x = iperm(x, m_symbolic->m_pinv);
1991  #if defined(G_DEBUG_SPARSE_COMPRESSION)
1992  std::cout << " Permutated vector : " << x << std::endl;
1993  #endif
1994 
1995  // Inplace solve L\x=x. The column and row maps are just use to see
1996  // which columns or rows should be skipped in the calculations.
1997  for (int col = 0; col < m_cols; ++col) { // loop over columns
1998  int c_col = col_map[col];
1999  if (c_col >= 0) { // use only non-zero cols
2000  x[c_col] /= Lx[Lp[col]]; // divide by diag.
2001  for (int p = Lp[col]+1; p < Lp[col+1]; p++) { // loop over elements
2002  int c_row = row_map[Li[p]];
2003  if (c_row >= 0) { // use only non-zero rows
2004  x[c_row] -= Lx[p] * x[c_col];
2005  }
2006  }
2007  }
2008  }
2009  #if defined(G_DEBUG_SPARSE_COMPRESSION)
2010  std::cout << " Solve Lx=x .......: " << x << std::endl;
2011  #endif
2012 
2013  // Inplace solve L'\x=x. The column and row maps are just use to see
2014  // which columns or rows should be skipped in the calculations.
2015  for (int col = m_cols-1; col >= 0; --col) { // loop over columns
2016  int c_col = col_map[col];
2017  if (c_col >= 0) { // use only non-zero cols
2018  for (int p = Lp[col]+1; p < Lp[col+1]; p++) { // loop over elements
2019  int c_row = row_map[Li[p]];
2020  if (c_row >= 0) { // use only non-zero rows
2021  x[c_col] -= Lx[p] * x[c_row];
2022  }
2023  }
2024  x[c_col] /= Lx[Lp[col]];
2025  }
2026  }
2027  #if defined(G_DEBUG_SPARSE_COMPRESSION)
2028  std::cout << " Solve L'x=x ......: " << x << std::endl;
2029  #endif
2030 
2031  // Perform vector permutation
2032  x = perm(x, m_symbolic->m_pinv);
2033  #if defined(G_DEBUG_SPARSE_COMPRESSION)
2034  std::cout << " Permutated vector : " << x << std::endl;
2035  #endif
2036 
2037  // If column compression has been performed the expand the result
2038  // vector accordingly
2039  if (m_colsel != NULL && m_num_colsel < m_cols) {
2040  for (int c_col = 0; c_col < m_num_colsel; ++c_col) {
2041  result[m_colsel[c_col]] = x[c_col];
2042  }
2043  }
2044  else {
2045  result = x;
2046  }
2047  #if defined(G_DEBUG_SPARSE_COMPRESSION)
2048  std::cout << " Restored vector ..: " << result << std::endl;
2049  #endif
2050 
2051  // Free mapping arrays
2052  delete [] row_map;
2053  delete [] col_map;
2054 
2055  } // endelse: Case B
2056 
2057  } // endif: matrix was not empty
2058 
2059  // Return result vector
2060  return result;
2061 }
2062 
2063 
2064 /***********************************************************************//**
2065  * @brief Invert matrix using a Cholesky decomposition
2066  *
2067  * @param[in] compress Use zero-row/column compression.
2068  * @return Inverted matrix.
2069  *
2070  * Inverts the matrix using a Cholesky decomposition.
2071  ***************************************************************************/
2073 {
2074  // Generate Cholesky decomposition of matrix
2075  GMatrixSparse decomposition = cholesky_decompose(compress);
2076 
2077  // Allocate result matrix and unit vector
2078  GMatrixSparse matrix(m_rows, m_cols);
2079  GVector unit(m_rows);
2080 
2081  // Column-wise solving of the problem
2082  for (int col = 0; col < m_cols; ++col) {
2083 
2084  // Set unit vector
2085  unit[col] = 1.0;
2086 
2087  // Solve for column
2088  GVector x = decomposition.cholesky_solver(unit, compress);
2089 
2090  // Insert column in matrix
2091  matrix.column(col, x);
2092 
2093  // Clear unit vector for next round
2094  unit[col] = 0.0;
2095 
2096  }
2097 
2098  // Return matrix
2099  return matrix;
2100 }
2101 
2102 
2103 /***********************************************************************//**
2104  * @brief Print matrix
2105  *
2106  * @param[in] chatter Chattiness.
2107  * @return String containing matrix information
2108  ***************************************************************************/
2109 std::string GMatrixSparse::print(const GChatter& chatter) const
2110 {
2111  // Initialise result string
2112  std::string result;
2113 
2114  // Continue only if chatter is not silent
2115  if (chatter != SILENT) {
2116 
2117  // Determine number of elements
2118  int nonzero = 0;
2119  int elements = m_elements;
2120  if (m_cols > 0) {
2121  if (m_fill_val == 0.0) {
2122  nonzero = (m_colstart != NULL) ? m_colstart[m_cols] : 0;
2123  }
2124  else {
2125  nonzero = (m_colstart != NULL) ? m_colstart[m_cols]+1 : 0;
2126  elements++;
2127  }
2128  }
2129 
2130  // Append header
2131  result.append("=== GMatrixSparse ===");
2132 
2133  // Append information
2134  result.append("\n"+gammalib::parformat("Number of rows"));
2135  result.append(gammalib::str(m_rows));
2136  if (m_rowsel != NULL) {
2137  result.append(" (compressed "+gammalib::str(m_num_rowsel)+")");
2138  }
2139  result.append("\n"+gammalib::parformat("Number of columns"));
2140  result.append(gammalib::str(m_cols));
2141  if (m_colsel != NULL) {
2142  result.append(" (compressed "+gammalib::str(m_num_colsel)+")");
2143  }
2144  result.append("\n"+gammalib::parformat("Number of nonzero elements"));
2145  result.append(gammalib::str(nonzero));
2146  if (m_fill_val != 0.0) {
2147  result.append("\n"+gammalib::parformat("Pending element"));
2148  result.append("("+gammalib::str(m_fill_row)+
2149  ","+gammalib::str(m_fill_col)+")=");
2150  result.append(gammalib::str(m_fill_val));
2151  }
2152  result.append("\n"+gammalib::parformat("Number of allocated cells"));
2153  result.append(gammalib::str(m_alloc));
2154  result.append("\n"+gammalib::parformat("Memory block size"));
2155  result.append(gammalib::str(m_mem_block));
2156  result.append("\n"+gammalib::parformat("Sparse matrix fill"));
2157  result.append(gammalib::str(fill()));
2158  result.append("\n"+gammalib::parformat("Pending element"));
2159  result.append(gammalib::str(m_fill_val));
2160  result.append("\n"+gammalib::parformat("Fill stack size"));
2161  result.append(gammalib::str(m_stack_size));
2162  if (m_stack_data == NULL) {
2163  result.append(" (none)");
2164  }
2165 
2166  // Append elements and compression schemes
2167  result.append(print_elements(chatter));
2168  result.append(print_row_compression(chatter));
2169  result.append(print_col_compression(chatter));
2170 
2171  } // endif: chatter was not silent
2172 
2173  // Return result
2174  return result;
2175 }
2176 
2177 
2178 /***********************************************************************//**
2179  * @brief Initialises matrix filling stack
2180  *
2181  * @param[in] size Stack size.
2182  * @param[in] entries Maximum number of entries.
2183  *
2184  * The matrix filling stack is used to allow for a fast column-wise filling
2185  * of a sparse matrix. Columns are successively appended to a stack which is
2186  * regularily flushed when it is full. This reduces memory copies and
2187  * movements and increases filling speed.
2188  ***************************************************************************/
2189 void GMatrixSparse::stack_init(const int& size, const int& entries)
2190 {
2191  // Free exisiting stack
2193 
2194  // Initialise stack members
2196  m_stack_max_entries = (entries > 0) ? entries : m_cols;
2197  m_stack_size = (size > 0) ? size : G_SPARSE_MATRIX_DEFAULT_STACK_SIZE;
2198 
2199  // Allocate stack memory. Raise an exception if allocation fails
2201  m_stack_start = new int[m_stack_max_entries+1];
2202  m_stack_data = new double[m_stack_size];
2203  m_stack_rowinx = new int[m_stack_size];
2204  m_stack_work = new int[m_cols];
2205  m_stack_rows = new int[m_cols];
2206  m_stack_values = new double[m_cols];
2207 
2208  // Initialise next free stack location to the first stack element
2209  m_stack_start[0] = 0;
2210 
2211  // Return
2212  return;
2213 }
2214 
2215 
2216 /***********************************************************************//**
2217  * @brief Push a vector column on the matrix stack
2218  *
2219  * @param[in] vector Vector.
2220  * @param[in] col Matrix column [0,...,columns()[.
2221  *
2222  * @exception GException::out_of_range
2223  * Invalid matrix column specified.
2224  * @exception GException::invalid_argument
2225  * Matrix dimension mismatches the vector size.
2226  *
2227  * Adds the contect of a vector to the matrix stack. This method is
2228  * identical to the GMatrixSparse::stack_push_column method that uses a
2229  * compressed array, yet it takes a full vector. Internal working arrays
2230  * are used to convert the full column vector in a compressed array and to
2231  * hand it over to the compressed array version.
2232  ***************************************************************************/
2233 int GMatrixSparse::stack_push_column(const GVector& vector, const int& col)
2234 {
2235  // Raise an exception if the column is invalid
2236  #if defined(G_RANGE_CHECK)
2237  if (col < 0 || col >= m_cols) {
2238  throw GException::out_of_range(G_STACK_PUSH1, "Matrix column",
2239  col, m_cols);
2240  }
2241  #endif
2242 
2243  // Raise an exception if the matrix and vector dimensions are incompatible
2244  if (vector.size() != m_rows) {
2245  std::string msg = "Vector size "+gammalib::str(vector.size())+
2246  " does not match "+gammalib::str(m_rows)+" matrix "
2247  "rows. Please specify a vector of size "+
2248  gammalib::str(m_rows)+".";
2250  }
2251 
2252  // Initialise number of non-zero elements
2253  int non_zero = 0;
2254 
2255  // Compress vector in the buffer and set-up row index array
2256  for (int i = 0; i < vector.size(); ++i) {
2257  if (vector[i] != 0.0) {
2258  m_stack_values[non_zero] = vector[i];
2259  m_stack_rows[non_zero] = i;
2260  non_zero++;
2261  }
2262  }
2263 
2264  // Hand over to compressed array method
2265  int remaining = stack_push_column(m_stack_values, m_stack_rows, non_zero, col);
2266 
2267  // Return remaining number of elements
2268  return remaining;
2269 }
2270 
2271 
2272 /***********************************************************************//**
2273  * @brief Push a compressed array on the matrix stack
2274  *
2275  * @param[in] values Compressed array.
2276  * @param[in] rows Matrix rows of array elements.
2277  * @param[in] number Number of elements in array.
2278  * @param[in] col Matrix column [0,...,columns()[.
2279  *
2280  * @exception GException::out_of_range
2281  * Invalid matrix column specified.
2282  * @exception GException::invalid_argument
2283  * Matrix dimension mismatches the vector size.
2284  *
2285  * The method puts new data on the stack while assuring that each column
2286  * is present only once. If an already existing column is encountered
2287  * the data from the existing and new column are mixed and put into a new
2288  * entry one the stack; the old entry is signalled as obsolete.
2289  *
2290  * On return the method indicates the number of elements in the input
2291  * array that have not been put onto the stack (due to memory limits).
2292  * This value can be either '0' (i.e. all elements have been put on the
2293  * stack) or 'number' (i.e. none of the elements have been put on the
2294  * stack). Columns are not partially put onto the stack.
2295  ***************************************************************************/
2296 int GMatrixSparse::stack_push_column(const double* values, const int* rows,
2297  const int& number, const int& col)
2298 {
2299  // Raise an exception if the column is invalid
2300  #if defined(G_RANGE_CHECK)
2301  if (col < 0 || col >= m_cols) {
2302  throw GException::out_of_range(G_STACK_PUSH2, "Matrix column",
2303  col, m_cols);
2304  }
2305  #endif
2306 
2307  // Initialise return value
2308  int remaining = number;
2309 
2310  // Single loop for common exit point
2311  do {
2312 
2313  // If the array is empty there is nothing to do
2314  if (!values || !rows || (number < 1)) {
2315  continue;
2316  }
2317 
2318  // Raise an exception if the matrix and vector dimensions are
2319  // incompatible. This test needs to be done after the test on a
2320  // positive number!
2321  if (rows[number-1] >= m_rows) {
2322  std::string msg = "Row number "+gammalib::str(rows[number-1])+" for "
2323  "element "+gammalib::str(number)+" is not smaller "
2324  "than number "+gammalib::str(m_rows)+" of matrix "
2325  "rows. Please specify only row numbers smaller than "+
2326  gammalib::str(m_rows)+".";
2328  }
2329 
2330  // If there is no stack or the stack can not hold the requested
2331  // number of elements then report number of array elements to the
2332  // caller
2333  if (m_stack_data == NULL || number > m_stack_size) {
2334  continue;
2335  }
2336 
2337  // Debug header
2338  #if defined(G_DEBUG_SPARSE_STACK_PUSH)
2339  std::cout << "GMatrixSparse::stack_push_column(v, i, n, col=" << col << ")";
2340  std::cout << std::endl;
2341  std::cout << " Data to push on stack ...:";
2342  for (int i = 0; i < number; ++i) {
2343  std::cout << " " << values[i];
2344  }
2345  std::cout << std::endl;
2346  std::cout << " Row indices of data .....:";
2347  for (int i = 0; i < number; ++i) {
2348  std::cout << " " << rows[i];
2349  }
2350  std::cout << std::endl;
2351  #endif
2352 
2353  // If the stack is full then flush it before pushing the array onto it.
2354  // There are 2 conditions that may fill the stack: all entries are
2355  // occupied or there is not enough space to hold more elements
2357  (number >= (m_stack_size - m_stack_start[m_stack_entries]))) {
2358  stack_flush();
2359  }
2360 
2361  // If the specified column is already on the stack then mix new column
2362  // with old one and invalidate old column (by setting its column index
2363  // to -1).
2364  // Since we loop here over all stack entries we don't want to have too
2365  // many entries in the stack. So don't set the maximum number of stack
2366  // entries too large ...
2367  for (int entry = 0; entry < m_stack_entries; ++entry) {
2368 
2369  // Is the specified column already on the stack?
2370  if (col == m_stack_colinx[entry]) {
2371 
2372  // Get start and stop indices of existing column in stack
2373  int i_start = m_stack_start[entry];
2374  int i_stop = m_stack_start[entry+1];
2375 
2376  // Allocate variables to hold mixing results
2377  int num_1; // Number of elements that are only in stack column
2378  int num_2; // Number of elements that are only in new column
2379  int num_mix; // Number of elements that are in both columns
2380 
2381  // Determine how many elements are requested to hold the combined
2382  // column
2383  mix_column_prepare(&(m_stack_rowinx[i_start]), i_stop-i_start,
2384  rows, number, &num_1, &num_2, &num_mix);
2385  int num_request = num_1 + num_2 + num_mix;
2386 
2387  // If there is not enough space in the stack to hold the combined
2388  // column we flush the stack and exit the loop now. In this case
2389  // the new column is put as the first column in the empty (flushed
2390  // stack)
2391  if (num_request >= (m_stack_size - m_stack_start[m_stack_entries])) {
2392  stack_flush();
2393  break;
2394  }
2395 
2396  // There is enough space, so combine both columns into a fresh one
2397  int inx = m_stack_start[m_stack_entries];
2398  int num_new;
2399  mix_column(&(m_stack_data[i_start]), &(m_stack_rowinx[i_start]),
2400  i_stop-i_start,
2401  values, rows, number,
2402  &(m_stack_data[inx]), &(m_stack_rowinx[inx]), &num_new);
2403 
2404  // Store entry information and initialise start of next entry
2405  m_stack_colinx[m_stack_entries] = col; // Store column index for entry
2406  m_stack_entries++; // Increase the # of entries
2407  m_stack_start[m_stack_entries] = inx + num_new; // Set start pointer for next entry
2408 
2409  // Invalidate old column
2410  m_stack_colinx[entry] = -1;
2411 
2412  // Fall through to end
2413  remaining = 0;
2414  break;
2415 
2416  } // endif: we found an existing column
2417  } // endfor: looped over all columns
2418 
2419  // If column has already been inserted then fall through ...
2420  if (remaining == 0) {
2421  continue;
2422  }
2423 
2424  // Otherwise, push the array on the stack
2425  int inx = m_stack_start[m_stack_entries];
2426  for (int i = 0; i < number; ++i) {
2427  m_stack_data[inx] = values[i];
2428  m_stack_rowinx[inx] = rows[i];
2429  inx++;
2430  }
2431 
2432  // Store entry information and initialise start of next entry
2433  m_stack_colinx[m_stack_entries] = col; // Store column index for entry
2434  m_stack_entries++; // Increase the # of entries
2435  m_stack_start[m_stack_entries] = inx; // Set start pointer for next entry
2436 
2437  // Signal success
2438  remaining = 0;
2439 
2440  } while (0); // End of main do-loop
2441 
2442  // Debug: show stack information
2443  #if defined(G_DEBUG_SPARSE_STACK_PUSH)
2444  std::cout << " Number of stack entries .: " << m_stack_entries << std::endl;
2445  std::cout << " Stack entry columns .....:";
2446  for (int i = 0; i < m_stack_entries; ++i) {
2447  std::cout << " " << m_stack_colinx[i];
2448  }
2449  std::cout << std::endl;
2450  std::cout << " Stack entry starts ......:";
2451  for (int i = 0; i < m_stack_entries; ++i) {
2452  std::cout << " " << m_stack_start[i];
2453  }
2454  std::cout << " (next at " << m_stack_start[m_stack_entries] << ")" << std::endl;
2455  std::cout << " Stack data ..............:";
2456  for (int i = 0; i < m_stack_start[m_stack_entries]; ++i) {
2457  std::cout << " " << m_stack_data[i];
2458  }
2459  std::cout << std::endl;
2460  std::cout << " Stack rows ..............:";
2461  for (int i = 0; i < m_stack_start[m_stack_entries]; ++i) {
2462  std::cout << " " << m_stack_rowinx[i];
2463  }
2464  std::cout << std::endl;
2465  #endif
2466 
2467  // Return remaining number of elements
2468  return remaining;
2469 }
2470 
2471 
2472 /***********************************************************************//**
2473  * @brief Flush matrix stack
2474  *
2475  * Adds the content of the stack to the actual matrix. First, the total
2476  * number of matrix elements is determined. Then new memory is allocated to
2477  * hold all elements. Finally, all elements are filled into the new memory
2478  * that is then replacing the old matrix memory.
2479  *
2480  * The method uses the internal working array m_stack_work.
2481  *
2482  * NOTE: The present algorithm assumes that each column occurs only once
2483  * in the stack!
2484  ***************************************************************************/
2486 {
2487  // Do nothing if there is no stack, no stack entries or the matrix is
2488  // empty
2489  if (m_stack_data == NULL || m_stack_entries == 0 || is_empty()) {
2490  return;
2491  }
2492 
2493  // Fill pending value
2494  fill_pending();
2495 
2496  // Debug header
2497  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
2498  std::cout << "GMatrixSparse::stack_flush" << std::endl;
2499  std::cout << " Number of columns on stack : " << m_stack_entries << std::endl;
2500  std::cout << " Number of elements on stack: " << m_stack_start[m_stack_entries] << std::endl;
2501  std::cout << " Number of matrix elements .: " << m_elements << std::endl;
2502  std::cout << " Col.start at end of matrix : " << m_colstart[m_cols] << std::endl;
2503  #endif
2504 
2505  // Use stack working array to flag all columns that exist already in the
2506  // matrix by 1 and all non-existing columns by 0
2507  for (int col = 0; col < m_cols; ++col) {
2508  if (m_colstart[col] < m_colstart[col+1]) {
2509  m_stack_work[col] = 1;
2510  }
2511  else {
2512  m_stack_work[col] = 0;
2513  }
2514  }
2515 
2516  // Initialise some element counters
2517  int new_elements = 0; // Number of new elements to add to matrix
2518  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
2519  int num_columns = 0; // Number of valid columns in stack
2520  int num_matrix = 0; // Number of elements only in matrix
2521  int num_stack = 0; // Number of elements only on stack
2522  int num_both = 0; // Number of elements in matrix and on stack
2523  #endif
2524 
2525  // Loop over all stack entries and gather the number of elements that are
2526  // new with respect to the initial matrix. For each column set a flag with
2527  // the following meanings:
2528  // -(entry+2): Column exists and is also found in stack=entry
2529  // +(entry+2): Column exists in stack=entry only
2530  // 1: Column exists in matrix only
2531  // 0: Column does neither exist in matrix or stack
2532  for (int entry = 0; entry < m_stack_entries; ++entry) {
2533 
2534  // Get column for actual entry
2535  int col = m_stack_colinx[entry];
2536 
2537  // Consider only valid entries
2538  if (col >= 0) {
2539 
2540  // Debug option: update number of valid columns
2541  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
2542  num_columns++;
2543  #endif
2544 
2545  // If column exists already in matrix then determine total number
2546  // of additional elements
2547  if (m_stack_work[col] == 1) {
2548 
2549  // Flag that this column is a mixed column
2550  m_stack_work[col] = -(entry+2);
2551 
2552  // Setup index boundaries
2553  int i_start = m_colstart[col];
2554  int i_stop = m_colstart[col+1];
2555  int k_start = m_stack_start[entry];
2556  int k_stop = m_stack_start[entry+1];
2557 
2558  // Allocate output counters
2559  int num_1; // Number of elements only found in matrix column
2560  int num_2; // Number of elements only found in stack column
2561  int num_mix; // Number of elements found in both columns
2562 
2563  // Analyse column mixing
2564  mix_column_prepare(&(m_rowinx[i_start]), i_stop-i_start,
2565  &(m_stack_rowinx[k_start]), k_stop-k_start,
2566  &num_1, &num_2, &num_mix);
2567 
2568  // Update element counters
2569  new_elements += num_2;
2570  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
2571  num_matrix += num_1;
2572  num_stack += num_2;
2573  num_both += num_mix;
2574  #endif
2575 
2576  } // endif: column existed in the matrix
2577 
2578  // If column did not exist in the matrix then consider all
2579  // elements as new
2580  else {
2581  m_stack_work[col] = (entry+2);
2582  new_elements += (m_stack_start[entry+1] - m_stack_start[entry]);
2583  }
2584 
2585  } // endif: entry was valid
2586 
2587  } // endfor: looped over all entries
2588 
2589  // Debug option: Log number of valid columns and elements in new matrix
2590  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
2591  std::cout << " Valid columns on stack ....: " << num_columns << std::endl;
2592  std::cout << " Valid elements on stack ...: " << new_elements << std::endl;
2593  #endif
2594 
2595  // Compute total number of matrix elements
2596  int elements = m_elements + new_elements;
2597 
2598  // Allocate memory for new matrix (always keep some elbow room)
2599  m_alloc = elements + m_mem_block;
2600  double* new_data = new double[m_alloc];
2601  int* new_rowinx = new int[m_alloc];
2602 
2603  // Fill new matrix. For this purpose we loop over all matrix columns
2604  // and perform the operation that was identified in the previous scan
2605  int index = 0;
2606  for (int col = 0; col < m_cols; ++col) {
2607 
2608  // If column does not exist then skip. Note that the index is stored
2609  // in [col] instead of [col+1] to not disturb the existing matrix.
2610  // The column start indices are corrected later.
2611  if (m_stack_work[col] == 0) {
2612  m_colstart[col] = index;
2613  continue;
2614  }
2615 
2616  // If column exists in the matrix but not in the stack we copy the
2617  // existing column into the new matrix
2618  if (m_stack_work[col] == 1) {
2619  for (int i = m_colstart[col]; i < m_colstart[col+1]; ++i) {
2620  new_data[index] = m_data[i];
2621  new_rowinx[index] = m_rowinx[i];
2622  index++;
2623  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
2624  num_matrix++;
2625  #endif
2626  }
2627  }
2628 
2629  // If column is new (i.e. it did not exist in the matrix before)
2630  // we copy the stack into the new matrix
2631  else if (m_stack_work[col] > 1) {
2632  int entry = m_stack_work[col] - 2;
2633  for (int i = m_stack_start[entry]; i < m_stack_start[entry+1]; ++i) {
2634  new_data[index] = m_stack_data[i];
2635  new_rowinx[index] = m_stack_rowinx[i];
2636  index++;
2637  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
2638  num_stack++;
2639  #endif
2640  }
2641  }
2642 
2643  // If column exists in the matrix and in the stack then we have to
2644  // mix both
2645  else {
2646 
2647  // Get stack entry for column mix
2648  int entry = -(m_stack_work[col] + 2);
2649 
2650  // Setup index boundaries
2651  int i_start = m_colstart[col];
2652  int i_stop = m_colstart[col+1];
2653  int k_start = m_stack_start[entry];
2654  int k_stop = m_stack_start[entry+1];
2655 
2656  // Allocate output counter
2657  int num;
2658 
2659  // Perform column mixing
2660  mix_column(&(m_data[i_start]), &(m_rowinx[i_start]), i_stop-i_start,
2661  &(m_stack_data[k_start]), &(m_stack_rowinx[k_start]), k_stop-k_start,
2662  &(new_data[index]), &(new_rowinx[index]), &num);
2663 
2664  // Increment element index
2665  index += num;
2666 
2667  } // endelse: column mixing required
2668 
2669  // Store actual index in column start array. Note that the index is
2670  // stored in [col] instead of [col+1] to not disturb the existing
2671  // matrix. The column start indices are corrected later.
2672  m_colstart[col] = index;
2673 
2674  } // endfor: looped over all columns
2675 
2676  // Dump number of elements in new matrix after stack flushing
2677  #if defined(G_DEBUG_SPARSE_STACK_FLUSH)
2678  std::cout << " Added elements ............: " << index;
2679  std::cout << " (should be " << elements << ")" << std::endl;
2680  std::cout << " - Matrix only .............: " << num_matrix << std::endl;
2681  std::cout << " - Stack only ..............: " << num_stack << std::endl;
2682  std::cout << " - Matrix & Stack ..........: " << num_both << std::endl;
2683  #endif
2684 
2685  // Correct columns start array
2686  for (int col = m_cols; col > 0; --col) {
2687  m_colstart[col] = m_colstart[col-1];
2688  }
2689  m_colstart[0] = 0;
2690 
2691  // Delete old matrix memory
2692  if (m_data != NULL) delete [] m_data;
2693  if (m_rowinx != NULL) delete [] m_rowinx;
2694 
2695  // Update pointers to new memory and update element counter
2696  m_data = new_data;
2697  m_rowinx = new_rowinx;
2698  m_elements = elements;
2699 
2700  // Stack is empty now, so reset stack counters
2701  m_stack_entries = 0;
2702  m_stack_start[0] = 0;
2703 
2704  // Return
2705  return;
2706 }
2707 
2708 
2709 /***********************************************************************//**
2710  * @brief Destroy matrix stack
2711  *
2712  * Flush and destroy matrix stack
2713  ***************************************************************************/
2715 {
2716  // Flush stack first
2717  stack_flush();
2718 
2719  // Free stack members
2721 
2722  // Initialise stack (no entries)
2724 
2725  // Return
2726  return;
2727 }
2728 
2729 
2730 /*==========================================================================
2731  = =
2732  = Private methods =
2733  = =
2734  ==========================================================================*/
2735 
2736 /***********************************************************************//**
2737  * @brief Initialise class mambers
2738  ***************************************************************************/
2740 {
2741  // Initialise sparse matrix members
2742  m_rowinx = NULL;
2744  m_zero = 0.0;
2745  m_fill_val = 0.0;
2746  m_fill_row = 0;
2747  m_fill_col = 0;
2748  m_symbolic = NULL;
2749  m_numeric = NULL;
2750 
2751  // Initialise stack members
2753 
2754  // Return
2755  return;
2756 }
2757 
2758 
2759 /***********************************************************************//**
2760  * @brief Copy class members
2761  *
2762  * @param[in] matrix Matrix to be copied.
2763  ***************************************************************************/
2765 {
2766  // Copy GMatrixSparse members
2767  m_mem_block = matrix.m_mem_block;
2768  m_zero = matrix.m_zero;
2769  m_fill_val = matrix.m_fill_val;
2770  m_fill_row = matrix.m_fill_row;
2771  m_fill_col = matrix.m_fill_col;
2772 
2773  // Copy row indices if they exist
2774  if (matrix.m_rowinx != NULL && matrix.m_alloc > 0) {
2775  m_rowinx = new int[matrix.m_alloc];
2776  for (int i = 0; i < matrix.m_elements; ++i) {
2777  m_rowinx[i] = matrix.m_rowinx[i];
2778  }
2779  }
2780 
2781  // Clone symbolic decomposition if it exists
2782  if (matrix.m_symbolic != NULL) {
2783  m_symbolic = new GSparseSymbolic();
2784  *m_symbolic = *matrix.m_symbolic;
2785  }
2786 
2787  // Copy numeric decomposition if it exists
2788  if (matrix.m_numeric != NULL) {
2789  m_numeric = new GSparseNumeric();
2790  *m_numeric = *matrix.m_numeric;
2791  }
2792 
2793  // Return
2794  return;
2795 }
2796 
2797 
2798 /***********************************************************************//**
2799  * @brief Delete class members
2800  ***************************************************************************/
2802 {
2803  // De-allocate only if memory has indeed been allocated
2804  if (m_numeric != NULL) delete m_numeric;
2805  if (m_symbolic != NULL) delete m_symbolic;
2806  if (m_rowinx != NULL) delete [] m_rowinx;
2807 
2808  // Properly mark members as free
2809  m_rowinx = NULL;
2810  m_symbolic = NULL;
2811  m_numeric = NULL;
2812 
2813  // Free stack members
2815 
2816  // Return
2817  return;
2818 }
2819 
2820 
2821 /***********************************************************************//**
2822  * @brief Allocate matrix
2823  *
2824  * @param[in] rows Number of rows.
2825  * @param[in] columns Number of columns.
2826  * @param[in] elements Number of matrix elements to be physically allocated.
2827  *
2828  * This is the main constructor code that allocates and initialises memory
2829  * for matrix elements. The method only allocates elements if both @p rows
2830  * and @p columns are positive. Otherwise the method does nothing and will
2831  * set the m_rows and m_cols attributes to zero.
2832  ***************************************************************************/
2833 void GMatrixSparse::alloc_members(const int& rows, const int& columns,
2834  const int& elements)
2835 {
2836  // Continue only if rows and columns are valid
2837  if (rows > 0 && columns > 0) {
2838 
2839  // Free any existing memory
2840  if (m_colstart != NULL) delete [] m_colstart;
2841 
2842  // Allocate column start array. This is the only array that we can
2843  // allocate at this time. The other arrays can only be allocated
2844  // during filling of the matrix
2845  m_colstart = new int[columns+1];
2846 
2847  // Store (logical) matrix size
2848  m_rows = rows;
2849  m_cols = columns;
2850 
2851  // Initialise column start indices to 0
2852  for (int col = 0; col <= m_cols; ++col) {
2853  m_colstart[col] = 0;
2854  }
2855 
2856  // Optionally allocate memory for matrix elements
2857  if (elements > 0) {
2858  alloc_elements(0, elements);
2859  }
2860 
2861  } // endif: number of elements was positive
2862 
2863  // Return
2864  return;
2865 }
2866 
2867 
2868 /***********************************************************************//**
2869  * @brief Initialise fill stack
2870  *
2871  * The fill stack is used to fill the sparse matrix without any prior know-
2872  * ledge about the number of location of the non-zero matrix elements.
2873  * Values are first filled into the stack and flushed into the sparse
2874  * matrix once the stack is full.
2875  ***************************************************************************/
2877 {
2878  // Initialise stack members
2879  m_stack_max_entries = 0;
2880  m_stack_size = 0;
2881  m_stack_entries = 0;
2882  m_stack_colinx = NULL;
2883  m_stack_start = NULL;
2884  m_stack_data = NULL;
2885  m_stack_rowinx = NULL;
2886  m_stack_work = NULL;
2887  m_stack_rows = NULL;
2888  m_stack_values = NULL;
2889 
2890  // Return
2891  return;
2892 }
2893 
2894 
2895 /***********************************************************************//**
2896  * @brief Delete fill-stack class members
2897  *
2898  * Deletes all memory that has been allocated using the init_stack_members()
2899  * method.
2900  ***************************************************************************/
2902 {
2903  // Free stack members
2904  if (m_stack_colinx != NULL) delete [] m_stack_colinx;
2905  if (m_stack_start != NULL) delete [] m_stack_start;
2906  if (m_stack_data != NULL) delete [] m_stack_data;
2907  if (m_stack_rowinx != NULL) delete [] m_stack_rowinx;
2908  if (m_stack_work != NULL) delete [] m_stack_work;
2909  if (m_stack_rows != NULL) delete [] m_stack_rows;
2910  if (m_stack_values != NULL) delete [] m_stack_values;
2911 
2912  // Properly mark members as free
2913  m_stack_colinx = NULL;
2914  m_stack_start = NULL;
2915  m_stack_data = NULL;
2916  m_stack_rowinx = NULL;
2917  m_stack_work = NULL;
2918  m_stack_rows = NULL;
2919  m_stack_values = NULL;
2920 
2921  // Return
2922  return;
2923 }
2924 
2925 
2926 /***********************************************************************//**
2927  * @brief Determines element index for (row,column)
2928  *
2929  * @param[in] row Matrix row [0,...,m_rows[.
2930  * @param[in] column Matrix column [0,...,m_cols[.
2931  * @return Element index.
2932  *
2933  * @exception GException::out_of_range
2934  * Matrix row or column out of range
2935  *
2936  * Returns the index in the compressed array for (row,col). The following
2937  * special results exist:
2938  *
2939  * -1: Requested index does not exist in the matrix.
2940  * m_elements: Requested index is the pending element.
2941  ***************************************************************************/
2942 int GMatrixSparse::get_index(const int& row, const int& column) const
2943 {
2944  // Raise an exception if the row or column index is invalid
2945  #if defined(G_RANGE_CHECK)
2946  if (row < 0 || row >= m_rows) {
2947  throw GException::out_of_range(G_GET_INDEX, "Matrix row", row, m_rows);
2948  }
2949  if (column < 0 || column >= m_cols) {
2950  throw GException::out_of_range(G_GET_INDEX, "Matrix column", column, m_cols);
2951  }
2952  #endif
2953 
2954  // Initialise element to 'not found'
2955  int index = -1;
2956 
2957  // If we have a pending element and the pending element is requested
2958  // then return m_elements as index
2959  if ((m_fill_val != 0.0) && (row == m_fill_row && column == m_fill_col)) {
2960  index = m_elements;
2961  }
2962 
2963  // ... otherwise if there are elements in the matrix then search for
2964  // the requested row and column.
2965  else if (m_elements > 0) {
2966 
2967  // Get the start and stop of the elements
2968  int* ptr_colstart = m_colstart + column;
2969  int i_start = *ptr_colstart++;
2970  int i_stop = *ptr_colstart;
2971 
2972  // If column is not empty then search requested row by bisection
2973  if (i_stop > i_start) {
2974 
2975  // Find row by bisection
2976  while ((i_stop - i_start) > 1) {
2977  int mid = (i_start+i_stop) / 2;
2978  if (m_rowinx[mid] > row) {
2979  i_stop = mid;
2980  }
2981  else {
2982  i_start = mid;
2983  }
2984  }
2985 
2986  // If we found the row then store the index
2987  if (m_rowinx[i_start] == row) {
2988  index = i_start;
2989  }
2990 
2991  } // endif: column was not empty
2992 
2993  } // endif: matrix contained elements
2994 
2995  // Return index
2996  return index;
2997 }
2998 
2999 
3000 /***********************************************************************//**
3001  * @brief Fills pending matrix element
3002  *
3003  * If 'm_fill_val' is non-zero a pending matrix element exists that should
3004  * be filled into (row,col)=(m_fill_row,m_fill_col). This routine performs
3005  * the filling of the matrix with this element and resets 'm_fill_val' to
3006  * zero. This routine allows for element-by-element filling of a sparse
3007  * matrix. This is, of course, very time consuming and should in general
3008  * be avoided. However, it allows to design a sparse matrix class that
3009  * hides the matrix sparsity completely to the user.
3010  ***************************************************************************/
3012 {
3013  // If we have a pending element then fill it into the matrix
3014  if (m_fill_val != 0.0) {
3015 
3016  // Debugging
3017  #if defined(G_DEBUG_SPARSE_PENDING) || defined(G_DEBUG_SPARSE_INSERTION)
3018  std::cout << "GMatrixSparse::fill_pending(): pending value " <<
3019  m_fill_val << " will be filled in location (" <<
3020  m_fill_row << "," << m_fill_col << ")" << std::endl;
3021  #endif
3022 
3023  // If there are so far no elements in the matrix then append element ...
3024  int inx_insert;
3025  if (m_elements == 0) {
3026  inx_insert = 0;
3027  }
3028 
3029  // ... otherwise search for index to insert
3030  else {
3031  int* ptr_colstart = m_colstart + m_fill_col;
3032  int i_start = *ptr_colstart++;
3033  int i_stop = *ptr_colstart;
3034  int* ptr_rowinx = m_rowinx + i_start;
3035  for (inx_insert = i_start; inx_insert < i_stop; ++inx_insert) {
3036  int row_test = *ptr_rowinx++;
3037  if (row_test > m_fill_row) {
3038  break;
3039  }
3040  }
3041  }
3042 
3043  // Allocate memory for new element
3044  alloc_elements(inx_insert, 1);
3045 
3046  // Insert element
3047  m_data[inx_insert] = m_fill_val;
3048  m_rowinx[inx_insert] = m_fill_row;
3049 
3050  // Update column start indices
3051  for (int col = m_fill_col+1; col <= m_cols; ++col) {
3052  m_colstart[col] += 1;
3053  }
3054 
3055  // Reset fill value
3056  m_fill_val = 0.0;
3057  m_fill_row = 0;
3058  m_fill_col = 0;
3059 
3060  // Debugging: show sparse matrix after filling
3061  #if defined(G_DEBUG_SPARSE_PENDING) || defined(G_DEBUG_SPARSE_INSERTION)
3062  std::cout << " Data: ";
3063  for (int i = 0; i < m_elements; ++i) {
3064  std::cout << m_data[i] << " ";
3065  }
3066  std::cout << std::endl << " Row.: ";
3067  for (int i = 0; i < m_elements; ++i) {
3068  std::cout << m_rowinx[i] << " ";
3069  }
3070  std::cout << std::endl << " Col.: ";
3071  for (int i = 0; i < m_cols+1; ++i) {
3072  std::cout << m_colstart[i] << " ";
3073  }
3074  std::cout << std::endl;
3075  #endif
3076 
3077  } // endif: a pending matrix element was found
3078 
3079  // Return
3080  return;
3081 }
3082 
3083 
3084 /***********************************************************************//**
3085  * @brief Allocate memory for new matrix elements
3086  *
3087  * @param[in] start Index of first allocated element.
3088  * @param[in] num Number of elements to be allocated.
3089  *
3090  * Inserts a memory allocation for 'num' elements at the index 'start'.
3091  * The new elements are filled with 0.0 and the corresponding row indices
3092  * are set to 0.
3093  *
3094  * NOTE: This method does not take care of updating 'm_colstart'. This has
3095  * to be done by the client.
3096  ***************************************************************************/
3097 void GMatrixSparse::alloc_elements(int start, const int& num)
3098 {
3099  // Dump header
3100  #if defined(G_DEBUG_SPARSE_MALLOC)
3101  std::cout << "GMatrixSparse::alloc_elements(start=" << start << ", num=" <<
3102  num << ")" << std::endl;
3103  std::cout << " Before allocation : " << m_elements << " " << m_alloc << std::endl;
3104  #endif
3105 
3106  // Continue only if we need memory
3107  if (num > 0) {
3108 
3109  // If start is after the end then append memory
3110  if (start > m_elements) {
3111  start = m_elements;
3112  }
3113 
3114  // Determine the requested new logical size of the matrix
3115  int new_size = m_elements + num;
3116 
3117  // Case A: the requested memory is already available, so just move the
3118  // data to make space for new elements and initialise the new cells
3119  if (new_size <= m_alloc) {
3120 
3121  // Move up all elements after index to insert
3122  for (int i = m_elements - 1; i >= start; --i) {
3123  m_data[i+num] = m_data[i];
3124  m_rowinx[i+num] = m_rowinx[i];
3125  }
3126 
3127  // Clear new elements (zero content for row 0)
3128  for (int i = start; i < start+num; ++i) {
3129  m_data[i] = 0.0;
3130  m_rowinx[i] = 0;
3131  }
3132 
3133  // Update element counter
3134  m_elements += num;
3135 
3136  } // endif: Case A: memory already existed
3137 
3138  // Case B: more memory is needed, so allocate it, copy the existing
3139  // content, and initialise new cells
3140  else {
3141 
3142  // Make sure that enough memory is allocated
3143  int new_propose = m_alloc + m_mem_block;
3144  m_alloc = (new_size > new_propose) ? new_size : new_propose;
3145 
3146  // Allocate memory for new elements
3147  double* new_data = new double[m_alloc];
3148  int* new_rowinx = new int[m_alloc];
3149 
3150  // Copy all elements before index to insert
3151  for (int i = 0; i < start; ++i) {
3152  new_data[i] = m_data[i];
3153  new_rowinx[i] = m_rowinx[i];
3154  }
3155 
3156  // Clear new elements (zero content for row 0)
3157  for (int i = start; i < start+num; ++i) {
3158  new_data[i] = 0.0;
3159  new_rowinx[i] = 0;
3160  }
3161 
3162  // Copy all elements after index to insert
3163  for (int i = start; i < m_elements; ++i) {
3164  new_data[i+num] = m_data[i];
3165  new_rowinx[i+num] = m_rowinx[i];
3166  }
3167 
3168  // Delete old memory
3169  if (m_data != NULL) delete [] m_data;
3170  if (m_rowinx != NULL) delete [] m_rowinx;
3171 
3172  // Update pointers to new memory and update element counter
3173  m_data = new_data;
3174  m_rowinx = new_rowinx;
3175  m_elements += num;
3176 
3177  } // endelse: Case B: more memory was needed
3178 
3179  } // endif: needed new memory
3180 
3181  // Dump new memory size
3182  #if defined(G_DEBUG_SPARSE_MALLOC)
3183  std::cout << " After allocation .: " << m_elements << " " << m_alloc << std::endl;
3184  #endif
3185 
3186  // Return
3187  return;
3188 }
3189 
3190 
3191 /***********************************************************************//**
3192  * @brief Free memory for obsolete matrix elements
3193  *
3194  * @param[in] start Index of first freed element.
3195  * @param[in] num Number of elements to be freed.
3196  *
3197  * Free memory for 'num' elements starting from index 'start'.
3198  *
3199  * NOTE: This method does not take care of updating 'm_colstart'. This has
3200  * to be done by the client.
3201  ***************************************************************************/
3202 void GMatrixSparse::free_elements(const int& start, const int& num)
3203 {
3204  // Dump header
3205  #if defined(G_DEBUG_SPARSE_MALLOC)
3206  std::cout << "GMatrixSparse::free_elements(start=" << start << ", num=" <<
3207  num << ")" << std::endl;
3208  #endif
3209 
3210  // Continue only if we need to free memory and if start is within the
3211  // range
3212  if (num > 0 && start < m_elements) {
3213 
3214  // Determine the requested new logical size of the matrix
3215  int new_size = m_elements - num;
3216 
3217  // If there are no elements then simply delete all matrix elements ...
3218  if (new_size < 1) {
3219  if (m_data != NULL) delete [] m_data;
3220  if (m_rowinx != NULL) delete [] m_rowinx;
3221  m_data = NULL;
3222  m_rowinx = NULL;
3223  m_elements = 0;
3224  m_alloc = 0;
3225  }
3226 
3227  // ... otherwise shrink the array
3228  else {
3229 
3230  // Case A: If at least one entire memory block has been liberated then
3231  // physically shrink the matrix array
3232  if (m_alloc - new_size > m_mem_block) {
3233 
3234  // Shrink, but leave a memory block for possible future filling
3235  m_alloc = new_size + m_mem_block;
3236 
3237  // Allocate new memory
3238  double* new_data = new double[m_alloc];
3239  int* new_rowinx = new int[m_alloc];
3240 
3241  // Copy all elements before the starting index
3242  for (int i = 0; i < start; ++i) {
3243  new_data[i] = m_data[i];
3244  new_rowinx[i] = m_rowinx[i];
3245  }
3246 
3247  // Copy all elements after the starting index
3248  for (int i = start; i < new_size; ++i) {
3249  new_data[i] = m_data[i+num];
3250  new_rowinx[i] = m_rowinx[i+num];
3251  }
3252 
3253  // Delete old memory
3254  if (m_data != NULL) delete [] m_data;
3255  if (m_rowinx != NULL) delete [] m_rowinx;
3256 
3257  // Update pointers to new memory and update element counter
3258  m_data = new_data;
3259  m_rowinx = new_rowinx;
3260  m_elements = new_size;
3261 
3262  } // endif: Case A: memory shrinkage performed
3263 
3264  // Case B: we keep the memory and just move the elements
3265  else {
3266 
3267  // Move all elements after the starting index
3268  for (int i = start; i < new_size; ++i) {
3269  m_data[i] = m_data[i+num];
3270  m_rowinx[i] = m_rowinx[i+num];
3271  }
3272 
3273  // Update element counter
3274  m_elements = new_size;
3275 
3276  } // endelse: Case B
3277 
3278  } // endif: array shrinkage needed
3279  } // endif: needed new memory
3280 
3281  // Return
3282  return;
3283 }
3284 
3285 
3286 /***********************************************************************//**
3287  * @brief Remove rows and columns with zeros
3288  *
3289  * @exception GException::runtime_error
3290  * All matrix elements are zero.
3291  *
3292  * Remove all rows and columns that contain only zeros from matrix. This
3293  * function is needed for compressed matrix factorisation. The resulting
3294  * matrix has reduced size (number of rows and columns) and dimensions
3295  * (number of elements). Note that the physical memory is not reduced by
3296  * the method.
3297  ***************************************************************************/
3299 {
3300  // Dump header
3301  #if defined(G_DEBUG_SPARSE_COMPRESSION)
3302  std::cout << "GMatrixSparse::remove_zero_row_col" << std::endl;
3303  #endif
3304 
3305  // Fill pending value in matrix
3306  fill_pending();
3307 
3308  // Select non-zero rows and columns of matrix
3309  select_non_zero();
3310 
3311  // Stop if there is no compression
3312  if (m_rows == m_num_rowsel && m_cols == m_num_colsel) {
3313  return;
3314  }
3315 
3316  // Throw exception if all matrix elements are zero
3317  if (m_num_rowsel < 1 || m_num_colsel < 1) {
3318  std::string msg = "All matrix elements are zero.";
3320  }
3321 
3322  // Allocate row mapping array
3323  int* row_map = new int[m_rows];
3324 
3325  // Setup row mapping array that maps original matrix rows into compressed
3326  // matrix rows. An entry of -1 indicates that the row should be dropped
3327  for (int row = 0; row < m_rows; ++row) {
3328  row_map[row] = -1;
3329  }
3330  for (int c_row = 0; c_row < m_num_rowsel; ++c_row) {
3331  row_map[m_rowsel[c_row]] = c_row;
3332  }
3333 
3334  // Initialise pointers to compressed array
3335  double* d_data = m_data;
3336  int* d_rowinx = m_rowinx;
3337 
3338  // Initialise column start of first column to zero
3339  m_colstart[0] = 0;
3340 
3341  // Dump column start array
3342  #if defined(G_DEBUG_SPARSE_COMPRESSION)
3343  std::cout << " before compression (" << m_colstart[m_cols] << "):";
3344  for (int col = 0; col <= m_cols; ++col)
3345  std::cout << " " << m_colstart[col];
3346  std::cout << std::endl;
3347  #endif
3348 
3349  // Loop over all columns of compressed matrix
3350  for (int c_col = 0; c_col < m_num_colsel; ++c_col) {
3351 
3352  // Get index range for elements in original matrix
3353  int i_start = m_colstart[m_colsel[c_col]];
3354  int i_stop = m_colstart[m_colsel[c_col]+1];
3355 
3356  // Initialise number of element counter for the compressed column
3357  int num = 0;
3358 
3359  // Loop over all elements in original column
3360  for (int i = i_start; i < i_stop; ++i) {
3361  int c_row = row_map[m_rowinx[i]];
3362  if (c_row >= 0) {
3363  *d_data++ = m_data[i];
3364  *d_rowinx++ = c_row;
3365  num++;
3366  }
3367  }
3368 
3369  // Update column stop
3370  m_colstart[c_col+1] = m_colstart[c_col] + num;
3371 
3372  } // endfor: looped over all columns of compressed matrix
3373 
3374  // Dump column start array
3375  #if defined(G_DEBUG_SPARSE_COMPRESSION)
3376  std::cout << " after compression (" << m_colstart[m_num_colsel] << ") :";
3377  for (int c_col = 0; c_col <= m_num_colsel; ++c_col)
3378  std::cout << " " << m_colstart[c_col];
3379  std::cout << std::endl;
3380  #endif
3381 
3382  // Free row mapping array
3383  delete [] row_map;
3384 
3385  // Update matrix size
3386  m_rows = m_num_rowsel;
3387  m_cols = m_num_colsel;
3388 
3389  // Return
3390  return;
3391 }
3392 
3393 
3394 /***********************************************************************//**
3395  * @brief Insert zero rows and columns
3396  *
3397  * @param[in] rows Number of rows.
3398  * @param[in] cols Number of columns.
3399  *
3400  * Insert zero rows and columns into matrix. Since for a sparse matrix
3401  * this does not require any allocation of additional memory, the data are
3402  * not moved by this function, but the pointers are re-arranged.
3403  ***************************************************************************/
3404 void GMatrixSparse::insert_zero_row_col(const int& rows, const int& cols)
3405 {
3406  // Dump header
3407  #if defined(G_DEBUG_SPARSE_COMPRESSION)
3408  std::cout << "GMatrixSparse::insert_zero_row_col(" << rows << "," << cols <<
3409  ")" << std::endl;
3410  #endif
3411 
3412  // Fill pending value in matrix
3413  fill_pending();
3414 
3415  // Stop if there is no insertion
3416  if (m_rows == rows && m_cols == cols) {
3417  return;
3418  }
3419 
3420  // If row selection exists then restore row indices
3421  if (m_rowsel != NULL) {
3422  for (int i = 0; i < m_elements; ++i) {
3423  m_rowinx[i] = m_rowsel[m_rowinx[i]];
3424  }
3425  }
3426 
3427  // Dump column start array
3428  #if defined(G_DEBUG_SPARSE_COMPRESSION)
3429  std::cout << " before restoration (" << m_colstart[m_num_colsel] << "):";
3430  for (int c_col = 0; c_col <= m_num_colsel; ++c_col)
3431  std::cout << " " << m_colstart[c_col];
3432  std::cout << std::endl;
3433  #endif
3434 
3435  // If column selection exists then restore column counters
3436  if (m_colsel != NULL) {
3437 
3438  // Start insertion from the last original column
3439  int col_stop = cols - 1;
3440 
3441  // Loop over all compressed columns
3442  for (int c_col = m_num_colsel-1; c_col > 0; --c_col) {
3443  int col_start = m_colsel[c_col-1] + 1;
3444  int col_value = m_colstart[c_col];
3445  for (int col = col_start; col <= col_stop; ++col) {
3446  m_colstart[col] = col_value;
3447  }
3448  col_stop = col_start - 1;
3449  }
3450 
3451  // Set first columns if they are not yet set
3452  for (int col = 0; col <= col_stop; ++col) {
3453  m_colstart[col] = 0;
3454  }
3455 
3456  // Restore the number of elements
3457  m_colstart[cols] = m_elements;
3458  }
3459 
3460  // Dump column start array
3461  #if defined(G_DEBUG_SPARSE_COMPRESSION)
3462  std::cout << " after restoration (" << m_colstart[cols] << ") :";
3463  for (int col = 0; col <= cols; ++col) {
3464  std::cout << " " << m_colstart[col];
3465  }
3466  std::cout << std::endl;
3467  #endif
3468 
3469  // Update matrix size
3470  m_rows = rows;
3471  m_cols = cols;
3472 
3473  // Return
3474  return;
3475 }
3476 
3477 
3478 /***********************************************************************//**
3479  * @brief Prepare mix of sparse columns
3480  *
3481  * @param[in] src1_row Row index array [0...src1_num-1] of first column.
3482  * @param[in] src1_num Number of elements in first column.
3483  * @param[in] src2_row Row index array [0...src2_num-1] of second column.
3484  * @param[in] src2_num Number of elements in second column.
3485  * @param[out] num_1 Number of elements only found in column 1.
3486  * @param[out] num_2 Number of elements only found in column 2.
3487  * @param[out] num_mix Number of elements found in both columns.
3488  *
3489  * This method prepares the mix of two sparse matrix columns into a single
3490  * column. It counts the number of elements in two columns that have the
3491  * same row and puts this number into @p num_mix. The number of elements
3492  * that are only present in the first column will be put in @p num_1, the
3493  * number of elements that are only present in the second column will be put
3494  * in @p num_2.
3495  ***************************************************************************/
3496 void GMatrixSparse::mix_column_prepare(const int* src1_row, int src1_num,
3497  const int* src2_row, int src2_num,
3498  int* num_1, int* num_2, int* num_mix)
3499 {
3500  // Initialise element counters
3501  *num_1 = 0;
3502  *num_2 = 0;
3503  *num_mix = 0;
3504 
3505  // Initialise indices and row indices of both columns
3506  int inx_1 = 0; // Column 1 element index
3507  int inx_2 = 0; // Column 2 element index
3508  int row_1 = src1_row[inx_1]; // Column 1 first row index
3509  int row_2 = src2_row[inx_2]; // Column 2 first row index
3510 
3511  // Count number of elements with same row or elements that exist only in
3512  // either the first or the second column while both columns contain still
3513  // elements
3514  while (inx_1 < src1_num && inx_2 < src2_num) {
3515 
3516  // Case A: the element exist in both columns
3517  if (row_1 == row_2) {
3518  (*num_mix)++;
3519  inx_1++;
3520  inx_2++;
3521  if (inx_1 < src1_num) {
3522  row_1 = src1_row[inx_1];
3523  }
3524  if (inx_2 < src2_num) {
3525  row_2 = src2_row[inx_2];
3526  }
3527  }
3528 
3529  // Case B: the element exists only in first column
3530  else if (row_1 < row_2) {
3531  (*num_1)++;
3532  inx_1++;
3533  if (inx_1 < src1_num) {
3534  row_1 = src1_row[inx_1];
3535  }
3536  }
3537 
3538  // Case C: the element exists only in second column
3539  else {
3540  (*num_2)++;
3541  inx_2++;
3542  if (inx_2 < src2_num) {
3543  row_2 = src2_row[inx_2];
3544  }
3545  }
3546 
3547  } // endwhile: counted elements
3548 
3549  // At this point at least one column expired of elements. If there are
3550  // still elements remaining in the first column then count them now.
3551  if (inx_1 < src1_num) {
3552  *num_1 += (src1_num - inx_1);
3553  }
3554 
3555  // If there are still elements remaining in the second column then count
3556  // them now.
3557  if (inx_2 < src2_num) {
3558  *num_2 += (src2_num - inx_2);
3559  }
3560 
3561  // We're done
3562  return;
3563 }
3564 
3565 
3566 /***********************************************************************//**
3567  * @brief Mix of sparse columns
3568  *
3569  * @param[in] src1_data Data array [0...src1_num-1] of first column.
3570  * @param[in] src1_row Row index array [0...src1_num-1] of first column.
3571  * @param[in] src1_num Number of elements in first column.
3572  * @param[in] src2_data Data array [0...src2_num-1] of second column.
3573  * @param[in] src2_row Row index array [0...src2_num-1] of second column.
3574  * @param[in] src2_num Number of elements in second column.
3575  * @param[in] dst_data Data array [0...dst_num-1] of mixed column.
3576  * @param[in] dst_row Row index array [0...dst_num-1] of mixed column.
3577  * @param[out] dst_num Number of elements in mixed column.
3578  *
3579  * This method mixes two sparse matrix columns into a single column.
3580  ***************************************************************************/
3581 void GMatrixSparse::mix_column(const double* src1_data,
3582  const int* src1_row,
3583  int src1_num,
3584  const double* src2_data,
3585  const int* src2_row,
3586  int src2_num,
3587  double* dst_data,
3588  int* dst_row,
3589  int* dst_num)
3590 {
3591  // Initialise indices and row indices of both columns
3592  int inx = 0; // Mixed column element index
3593  int inx_1 = 0; // Column 1 element index
3594  int inx_2 = 0; // Column 2 element index
3595  int row_1 = src1_row[inx_1]; // Current row in column 1
3596  int row_2 = src2_row[inx_2]; // Current row in column 2
3597 
3598  // Mix elements of both columns while both contain still elements
3599  while (inx_1 < src1_num && inx_2 < src2_num) {
3600 
3601  // Case A: the element exists in both columns, so we add the values
3602  if (row_1 == row_2) {
3603  dst_data[inx] = src1_data[inx_1++] + src2_data[inx_2++];
3604  dst_row[inx] = row_1;
3605  if (inx_1 < src1_num) {
3606  row_1 = src1_row[inx_1];
3607  }
3608  if (inx_2 < src2_num) {
3609  row_2 = src2_row[inx_2];
3610  }
3611  }
3612 
3613  // Case B: the element exists only in first column, so we copy the
3614  // element from the first column
3615  else if (row_1 < row_2) {
3616  dst_data[inx] = src1_data[inx_1++];
3617  dst_row[inx] = row_1;
3618  if (inx_1 < src1_num) {
3619  row_1 = src1_row[inx_1];
3620  }
3621  }
3622 
3623  // Case C: the element exists only in second column, so we copy the
3624  // element from the second column
3625  else {
3626  dst_data[inx] = src2_data[inx_2++];
3627  dst_row[inx] = row_2;
3628  if (inx_2 < src2_num) {
3629  row_2 = src2_row[inx_2];
3630  }
3631  }
3632 
3633  // Update the destination index since we added a element
3634  inx++;
3635 
3636  } // endwhile: mixing
3637 
3638  // At this point either the first or the second column expired of elements
3639  // In the case that there are still elements remaining in the first column
3640  // we add them now ...
3641  for (int i = inx_1; i < src1_num; ++i, ++inx) {
3642  dst_data[inx] = src1_data[i];
3643  dst_row[inx] = src1_row[i];
3644  }
3645 
3646  // ... or in the case that there are still elements remaining in the second
3647  // column we add them now
3648  for (int i = inx_2; i < src2_num; ++i, ++inx) {
3649  dst_data[inx] = src2_data[i];
3650  dst_row[inx] = src2_row[i];
3651  }
3652 
3653  // Now store the number of elements in the mixed column
3654  *dst_num = inx;
3655 
3656  // We're done
3657  return;
3658 }
3659 
3660 
3661 /***********************************************************************//**
3662  * @brief Insert row in matrix
3663  *
3664  * @param[in] row Matrix row [0,...,row()[.
3665  * @param[in] vector Vector.
3666  * @param[in] add Add vector to existing elements?
3667  *
3668  * @exception GException::out_of_range
3669  * Invalid matrix row specified.
3670  * @exception GException::invalid_argument
3671  * Vector size incompatible with number of matrix columns.
3672  *
3673  * @todo Remove elements that are empty after addition
3674  ***************************************************************************/
3675 void GMatrixSparse::insert_row(const int& row,
3676  const GVector& vector,
3677  const bool& add)
3678 {
3679  // Raise an exception if the row index is invalid
3680  #if defined(G_RANGE_CHECK)
3681  if (row < 0 || row >= m_rows) {
3682  throw GException::out_of_range(G_SET_ROW, "Matrix row", row, m_rows);
3683  }
3684  #endif
3685 
3686  // Raise an exception if the vector length is invalid
3687  if (vector.size() != m_cols) {
3688  std::string msg = "Vector size "+gammalib::str(vector.size())+
3689  " does not match "+gammalib::str(m_cols)+" matrix "
3690  "columns. Please specify a vector of size "+
3691  gammalib::str(m_cols)+".";
3693  }
3694 
3695  // Fill pending element into matrix
3696  fill_pending();
3697 
3698  // Determine indices of new columns that need to be allocated and set
3699  // or add all vector elements that do not need a new memory allocation
3700  std::vector<int> new_columns;
3701  std::vector<int> k_insert;
3702  for (int i = m_cols-1; i >= 0; --i) {
3703  if (vector[i] != 0.0) {
3704  int k = m_colstart[i];
3705  int istop = m_colstart[i+1];
3706  bool add_column = true;
3707  for (; k < istop; ++k) {
3708  if (m_rowinx[k] == row) {
3709  if (add) {
3710  m_data[k] += vector[i];
3711  }
3712  else {
3713  m_data[k] = vector[i];
3714  }
3715  add_column = false;
3716  break;
3717  }
3718  else if (m_rowinx[k] > row) {
3719  break;
3720  }
3721  }
3722  if (add_column) {
3723  new_columns.push_back(i);
3724  k_insert.push_back(k);
3725  }
3726  }
3727  }
3728 
3729  // If new columns are needed then allocate new memory and insert
3730  // new columns
3731  if (!new_columns.empty()) {
3732 
3733  // Determine the requested new logical size of the matrix
3734  int new_size = m_elements + new_columns.size();
3735 
3736  // Initialise pointers to destination
3737  bool new_memory = false;
3738  double* new_data = m_data;
3739  int* new_rowinx = m_rowinx;
3740 
3741  // If requested size is smaller than allocated size then allocate
3742  // new memory
3743  if (new_size > m_alloc) {
3744 
3745  // Propose a new memory size
3746  int new_propose = m_alloc + m_mem_block;
3747 
3748  // Make sure that enough memory is allocated
3749  m_alloc = (new_size > new_propose) ? new_size : new_propose;
3750 
3751  // Allocate memory for new elements
3752  new_data = new double[m_alloc];
3753  new_rowinx = new int[m_alloc];
3754 
3755  // Signal that new memory was allocated
3756  new_memory = true;
3757 
3758  } // endif: memory allocation required
3759 
3760  // Loop over all new columns
3761  int n_insert = new_columns.size();
3762  int k_last = m_elements;
3763  for (int i = 0; i < new_columns.size(); ++i) {
3764 
3765  // Move data back by n_insert positions
3766  for (int k = k_insert[i]; k < k_last; ++k) {
3767  new_data[k+n_insert] = m_data[k];
3768  new_rowinx[k+n_insert] = m_rowinx[k];
3769  }
3770 
3771  // Insert vector element just before the block that was moved
3772  // back
3773  new_data[k_insert[i]+n_insert-1] = vector[new_columns[i]];
3774  new_rowinx[k_insert[i]+n_insert-1] = row;
3775 
3776  // Update number of positions to move back and end of block
3777  n_insert--;
3778  k_last = k_insert[i];
3779 
3780  } // endfor: looped over all new columns
3781 
3782  // Update column start and number of elements
3783  int i_insert = 0;
3784  n_insert = new_columns.size();
3785  for (int i = m_cols; i > 0; --i) {
3786  if (new_columns[i_insert] == i) {
3787  n_insert--;
3788  i_insert++;
3789  if (n_insert < 1) {
3790  break;
3791  }
3792  }
3793  m_colstart[i] += n_insert;
3794  }
3795 
3796  // Update number of elements
3797  m_elements = new_size;
3798 
3799  // If memory was allocated then free old memory and attach new
3800  // memory
3801  if (new_memory) {
3802 
3803  // Delete old memory
3804  if (m_data != NULL) delete [] m_data;
3805  if (m_rowinx != NULL) delete [] m_rowinx;
3806 
3807  // Update pointers to new memory and update element counter
3808  m_data = new_data;
3809  m_rowinx = new_rowinx;
3810 
3811  } // endif: memory was allocated
3812 
3813  } // endif: new columns were needed
3814 
3815  // Return
3816  return;
3817 }
3818 
3819 
3820 /*==========================================================================
3821  = =
3822  = Friend functions =
3823  = =
3824  ==========================================================================*/
3825 
3826 /***********************************************************************//**
3827  * @brief cs_symperm
3828  *
3829  * @param[in] matrix Matrix.
3830  * @param[in] pinv TBD.
3831  *
3832  * Returns matrix(p,p) where matrix and matrix(p,p) are symmetric the upper
3833  * part stored.
3834  ***************************************************************************/
3835 GMatrixSparse cs_symperm(const GMatrixSparse& matrix, const int* pinv)
3836 {
3837  // Declare loop variables
3838  int i, j, p, q, i2, j2;
3839 
3840  // Assign matrix attributes
3841  int n = matrix.m_cols;
3842  int* Ap = matrix.m_colstart;
3843  int* Ai = matrix.m_rowinx;
3844  double* Ax = matrix.m_data;
3845 
3846  // Allocate result matrix
3847  GMatrixSparse C(n, n, Ap[n]);
3848 
3849  // Allocate and initialise workspace
3850  int wrk_size = n;
3851  int* wrk_int = new int[wrk_size];
3852  for (i = 0; i < wrk_size; ++i) {
3853  wrk_int[i] = 0;
3854  }
3855 
3856  // Assign result matrix attributes
3857  int* Cp = C.m_colstart;
3858  int* Ci = C.m_rowinx;
3859  double* Cx = C.m_data;
3860 
3861  // Count entries in each column of C
3862  for (j = 0; j < n; j++) {
3863 
3864  // Column j of A is column j2 of C
3865  j2 = pinv ? pinv[j] : j;
3866 
3867  // Loop over entries in column j
3868  for (p = Ap[j]; p < Ap[j+1]; p++) {
3869  i = Ai [p];
3870  if (i > j) continue; // skip lower triangular part of A
3871  i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
3872  wrk_int[G_MAX(i2, j2)]++; // column count of C
3873  }
3874  }
3875 
3876  // Compute column pointers of C
3877  cs_cumsum(Cp, wrk_int, n);
3878 
3879  // Loop over all columns of A
3880  for (j = 0 ; j < n ; j++) {
3881 
3882  // Column j of A is column j2 of C
3883  j2 = pinv ? pinv[j] : j;
3884 
3885  // Loop over entries in column j
3886  for (p = Ap[j]; p < Ap[j+1]; p++) {
3887  i = Ai [p] ;
3888  if (i > j) continue; // skip lower triangular part of A
3889  i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
3890  Ci[q = wrk_int[G_MAX(i2,j2)]++] = G_MIN(i2,j2);
3891  if (Cx) Cx[q] = Ax[p];
3892  }
3893  }
3894 
3895  // Free workspace
3896  delete [] wrk_int;
3897 
3898  // Rectify the number of elements in matrix C
3899  C.free_elements(Cp[n], (C.m_elements-Cp[n]));
3900 
3901  // Return result
3902  return C;
3903 }
3904 
3905 
3906 /***************************************************************************
3907  * @brief Compute transpose matrix
3908  *
3909  * @param[in] matrix Matrix.
3910  * @param[in] values Flag that signals if values should be copied.
3911  *
3912  * The flag 'values' allows to avoid copying the actual data values. This
3913  * allows to perform a logical matrix transposition, as needed by the
3914  * symbolic matrix analysis class.
3915  *
3916  * Note that this method does not support pending elements (they have to be
3917  * filled before).
3918  *
3919  * The method does nothing on empty matrices.
3920  ***************************************************************************/
3921 GMatrixSparse cs_transpose(const GMatrixSparse& matrix, int values)
3922 {
3923  // Declare and allocate result matrix
3924  GMatrixSparse result(matrix.m_cols, matrix.m_rows, matrix.m_elements);
3925 
3926  // Transpose only if matrix is not empty
3927  if (!matrix.is_empty()) {
3928 
3929  // Allocate and initialise workspace
3930  int wrk_size = matrix.m_rows;
3931  int* wrk_int = new int[wrk_size];
3932  for (int i = 0; i < wrk_size; ++i) {
3933  wrk_int[i] = 0;
3934  }
3935 
3936  // Setup the number of non-zero elements in each row
3937  // for (p = 0 ; p < Ap [n] ; p++) w [Ai [p]]++ ;
3938  // Ap[n] = m.m_colstart[col]
3939  // n = m.m_cols
3940  // Ai[p] = m.m_rowinx[p]
3941  for (int p = 0; p < matrix.m_colstart[matrix.m_cols]; ++p) {
3942  wrk_int[matrix.m_rowinx[p]]++;
3943  }
3944 
3945  // Set row pointers. To use a GSparseSymbolic function we have to
3946  // allocate and object (but this does not take memory)
3947  cs_cumsum(result.m_colstart, wrk_int, matrix.m_rows);
3948 
3949  // Case A: Normal transponse, including assignment of values
3950  if (values) {
3951  for (int col = 0; col < matrix.m_cols; ++col) {
3952  for (int p = matrix.m_colstart[col]; p < matrix.m_colstart[col+1] ; ++p) {
3953  int i = wrk_int[matrix.m_rowinx[p]]++;
3954  result.m_rowinx[i] = col;
3955  result.m_data[i] = matrix.m_data[p] ;
3956  }
3957  }
3958  }
3959 
3960  // Case B: Logical transponse, no assignment of values is performed
3961  else {
3962  for (int col = 0; col < matrix.m_cols; ++col) {
3963  for (int p = matrix.m_colstart[col]; p < matrix.m_colstart[col+1] ; ++p) {
3964  result.m_rowinx[wrk_int[matrix.m_rowinx[p]]++] = col;
3965  }
3966  }
3967  }
3968 
3969  // Free workspace
3970  delete [] wrk_int;
3971 
3972  } // endif: matrix was not empty
3973 
3974  // Return transponse matrix
3975  return result;
3976 }
3977 
3978 
3979 /***********************************************************************//**
3980  * @brief cs_cumsum
3981  *
3982  * @param[out] p Integer array (n+1 elements).
3983  * @param[in] c Integer array (n elements).
3984  * @param[in] n Number of elements.
3985  *
3986  * Evaluate p[0..n] = cumulative sum of c[0..n-1]
3987  ***************************************************************************/
3988 double cs_cumsum(int* p, int* c, int n)
3989 {
3990  // Signal error if one of the input pointers is NULL
3991  if (!p || !c) return (-1);
3992 
3993  // Initialise sums (integer and double)
3994  int nz = 0;
3995  double nz2 = 0.0;
3996 
3997  // Evaluate p[0..n] = cumulative sum of c[0..n-1]
3998  for (int i = 0; i < n; ++i) {
3999  p[i] = nz ;
4000  nz += c[i];
4001  nz2 += c[i]; // also in double to avoid int overflow
4002  c[i] = p[i]; // also copy p[0..n-1] back into c[0..n-1]
4003  }
4004  p[n] = nz ;
4005 
4006  // Return cumulative sum of c[0..n-1]
4007  return nz2;
4008 }
virtual double & operator()(const int &row, const int &column)
Return reference to matrix element.
virtual GMatrixSparse & operator=(const GMatrixSparse &matrix)
Matrix assignment operator.
virtual double max(void) const
Return maximum matrix element.
GMatrixSparse cs_symperm(const GMatrixSparse &matrix, const int *pinv)
cs_symperm
Abstract matrix base class interface definition.
#define G_EXTRACT_COLUMN
GMatrixSparse * m_L
int get_index(const int &row, const int &column) const
Determines element index for (row,column)
std::string number(const std::string &noun, const int &number)
Convert singular noun into number noun.
Definition: GTools.cpp:1167
std::string print_row_compression(const GChatter &chatter=NORMAL) const
Print row compression scheme if it exists.
void cholesky_symbolic_analysis(int order, const GMatrixSparse &m)
int m_fill_col
Column into which element needs to be filled.
Sparse matrix numeric analysis class definition.
void stack_destroy(void)
Destroy matrix stack.
const int & size(void) const
Return number of matrix elements.
GMatrixSparse cholesky_invert(const bool &compress=true) const
Invert matrix using a Cholesky decomposition.
Sparse matrix class interface definition.
#define G_CONSTRUCTOR
virtual double fill(void) const
Returns fill of matrix.
void remove_zero_row_col(void)
Remove rows and columns with zeros.
#define G_GET_INDEX
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
void insert_zero_row_col(const int &rows, const int &cols)
Insert zero rows and columns.
Generic matrix class definition.
#define G_CHOL_SOLVE
GMatrixSparse cs_transpose(const GMatrixSparse &matrix, int values)
int m_rows
Number of rows.
Sparse matrix symbolic analysis class definition.
virtual GMatrixSparse * clone(void) const
Clone matrix.
virtual GMatrixSparse & operator*=(const GMatrixSparse &matrix)
Unary matrix multiplication operator.
virtual GVector column(const int &column) const
Extract column as vector from matrix.
Definition: GMatrix.cpp:718
Symmetric matrix class interface definition.
virtual ~GMatrixSparse(void)
Destructor.
Gammalib tools definition.
GMatrixSparse invert(void) const
Return inverted matrix.
virtual double & at(const int &row, const int &column)
Return reference to matrix element.
int m_num_colsel
Number of selected columns (for compressed decomposition)
virtual GVector column(const int &column) const
Extract column as vector from matrix.
GMatrixSparse(void)
Void matrix constructor.
virtual void clear(void)
Clear matrix.
int m_stack_max_entries
Maximum number of entries in the stack.
#define G_ADD_TO_COLUMN
#define G_STACK_PUSH2
#define G_SET_ROW
int * m_stack_work
Stack flush integer working array [m_cols].
GMatrixSparse transpose(void) const
Return transposed matrix.
virtual void add_to_column(const int &column, const GVector &vector)
Add vector column into matrix.
Sparse matrix symbolic analysis class.
#define G_OP_MUL_VEC
#define G_OP_SUB
void insert_row(const int &row, const GVector &vector, const bool &add)
Insert row in matrix.
GVector solve(const GVector &vector) const
Solves linear matrix equation.
int * m_stack_rowinx
Stack row indices [m_stack_size].
void stack_flush(void)
Flush matrix stack.
std::string print_elements(const GChatter &chatter=NORMAL, const int &num=15) const
Print all matrix elements.
GVector perm(const GVector &vector, const int *p)
Computes vector permutation.
Definition: GVector.cpp:1011
friend class GSparseNumeric
double * m_data
Matrix data.
const int & rows(void) const
Return number of matrix rows.
int non_zeros(void) const
Returns number of non-zero elements in vector.
Definition: GVector.cpp:559
int * m_stack_colinx
Column index for each entry [m_stack_entries].
GMatrixSparse abs(void) const
Return absolute of matrix.
int * m_pinv
Inverse row permutation for QR, fill reduce permutation for Cholesky.
virtual GMatrixSparse & operator-=(const GMatrixSparse &matrix)
Unary matrix subtraction operator.
void free_elements(const int &start, const int &num)
Free memory for obsolete matrix elements.
int m_cols
Number of columns.
virtual double sum(void) const
Sum matrix elements.
double * m_stack_data
Stack data [m_stack_size].
double cs_cumsum(int *p, int *c, int n)
cs_cumsum
void free_stack_members(void)
Delete fill-stack class members.
void init_stack_members(void)
Initialise fill stack.
GChatter
Definition: GTypemaps.hpp:33
int * m_rowinx
Row-indices of all elements.
#define G_MAX(a, b)
int m_alloc
Size of allocated matrix memory.
virtual void add_to_row(const int &row, const GVector &vector)
Add row to matrix elements.
Sparse matrix numeric analysis class.
int m_stack_entries
Number of entries in the stack.
std::string print_col_compression(const GChatter &chatter=NORMAL) const
Print column compression scheme if it exists.
int * m_stack_start
Start in stack for each entry [m_stack_entries+1].
#define G_SET_COLUMN
Vector class interface definition.
#define G_OP_ADD
void copy_members(const GMatrixSparse &m)
Copy class members.
void fill_pending(void)
Fills pending matrix element.
#define G_SPARSE_MATRIX_DEFAULT_STACK_SIZE
virtual GVector row(const int &row) const
Extract row as vector from matrix.
#define G_AT
GSparseSymbolic * m_symbolic
Holds GSparseSymbolic object after decomposition.
double m_fill_val
Element to be filled.
#define G_INSERT_ROW
#define G_REMOVE_ZERO
Symmetric matrix class definition.
virtual double min(void) const
Return minimum matrix element.
#define G_SET_COLUMN2
virtual GMatrixSparse operator-(void) const
Negate matrix elements.
#define G_ADD_TO_COLUMN2
int m_elements
Number of elements stored in matrix.
void mix_column_prepare(const int *src1_row, int src1_num, const int *src2_row, int src2_num, int *num_1, int *num_2, int *num_mix)
Prepare mix of sparse columns.
virtual GVector column(const int &column) const
Extract column as vector from matrix.
Sparse matrix class definition.
virtual bool operator==(const GMatrixSparse &matrix) const
Equalty operator.
GVector cholesky_solver(const GVector &vector, const bool &compress=true) const
Cholesky solver.
friend GMatrixSparse cs_transpose(const GMatrixSparse &matrix, int values)
int * m_stack_rows
Stack push integer working array [m_cols].
int m_fill_row
Row into which element needs to be filled.
void stack_init(const int &size=0, const int &entries=0)
Initialises matrix filling stack.
Exception handler interface definition.
#define G_MIN(a, b)
virtual GMatrixBase & operator=(const GMatrixBase &matrix)
Assignment operator.
void init_members(void)
Initialise class mambers.
const int & size(void) const
Return size of vector.
Definition: GVector.hpp:178
GVector iperm(const GVector &vector, const int *p)
Computes vector inverse permutation.
Definition: GVector.cpp:1038
Generic matrix class definition.
Definition: GMatrix.hpp:79
Vector class.
Definition: GVector.hpp:46
virtual std::string print(const GChatter &chatter=NORMAL) const
Print matrix.
virtual GVector operator*(const GVector &vector) const
Vector multiplication.
int m_num_rowsel
Number of selected rows (for compressed decomposition)
void alloc_members(const int &rows, const int &cols, const int &elements=0)
Allocate matrix.
GSparseNumeric * m_numeric
Holds GSparseNumeric object after decomposition.
virtual GMatrixSparse & operator+=(const GMatrixSparse &matrix)
Unary matrix addition operator.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
friend class GSparseSymbolic
void cholesky_numeric_analysis(const GMatrixSparse &m, const GSparseSymbolic &s)
int * m_colstart
Column start indices (m_cols+1)
int * m_rowsel
Row selection (for compressed decomposition)
#define G_EXTRACT_ROW
void select_non_zero(void)
Setup compressed matrix factorisation.
#define G_SPARSE_MATRIX_DEFAULT_MEM_BLOCK
int m_mem_block
Memory block to be allocated at once.
const int & columns(void) const
Return number of matrix columns.
#define G_STACK_PUSH1
GMatrixSparse cholesky_decompose(const bool &compress=true) const
Return Cholesky decomposition.
int m_stack_size
Maximum number of elements in the stack.
int * m_colsel
Column selection (for compressed decomposition)
bool is_empty(void) const
Check if matrix is empty.
virtual bool operator!=(const GMatrixSparse &matrix) const
Non-equality operator.
double m_zero
The zero element (needed for data access)
void alloc_elements(int start, const int &num)
Allocate memory for new matrix elements.
double * m_stack_values
Stack push double buffer [m_cols].
#define G_OP_MAT_MUL
void mix_column(const double *src1_data, const int *src1_row, int src1_num, const double *src2_data, const int *src2_row, int src2_num, double *dst_data, int *dst_row, int *dst_num)
Mix of sparse columns.
void free_members(void)
Delete class members.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
int stack_push_column(const GVector &vector, const int &col)
Push a vector column on the matrix stack.