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