Justine Devin / gammalib

Commit 305067486e2cf831fc2edfa05960a4ed28936f07

Authored by Jürgen Knödlseder
1 parent 92d08c9c

Improve Doxygen documentation, add at() methods and value assignment operator.

Inline Side-by-side
Showing 5 changed files with 389 additions and 121 deletions
include/GMatrixBase.hpp
  View file @3050674
... ... @@ -38,9 +38,110 @@
38 38 *
39 39 * @brief Abstract matrix base class interface defintion
40 40 *
41   - * This class is an abstract base class from which all matrix classes will
42   - * be derived. It defines the common interface for all matrix classes and
43   - * implements common methods and data members.
  41 + * This class provides the abstract interface for all matrix classes. A
  42 + * GammaLib matrix is a two dimensional array of double precision values
  43 + * that is commonly used for linear algebra computations.
  44 + *
  45 + * The matrix classes were designed in a manner so that matrices can be used
  46 + * in a natural syntax that is comparable to the syntax used for scalar
  47 + * values. In particular, the following matrix operations are available:
  48 + *
  49 + * assigned_matrix = matrix;
  50 + * assigned_matrix = value;
  51 + * added_matrix = matrix1 + matrix2;
  52 + * added_matrix += matrix;
  53 + * subtracted_matrix = matrix1 - matrix2;
  54 + * subtracted_matrix -= matrix;
  55 + * multiplied_matrix = matrix1 * matrix2;
  56 + * multiplied_matrix *= matrix;
  57 + * scaled_matrix = matrix * scale;
  58 + * scaled_matrix = scale * matrix;
  59 + * scaled_matrix *= scale;
  60 + * divided_matrix = matrix / factor;
  61 + * divided_matrix /= factor;
  62 + * negated_matrix = -matrix;
  63 + *
  64 + * where @p matrix, @p matrix1 and @p matrix2 are matrices, and
  65 + * @p value, @p scale and @p factor are double precision values. For the
  66 + * value assignment operator, each matrix element is assigned the specified
  67 + * value. For the scaling and division operators, every matrix element is
  68 + * multiplied or divided by the specified value.
  69 + *
  70 + * Matrices can also be compared using
  71 + *
  72 + * if (matrix == another_matrix) ...
  73 + * if (matrix != another_matrix) ...
  74 + *
  75 + * where identify is defined as having the same number of rows and columns
  76 + * and every element having the same value.
  77 + *
  78 + * Matrix elements are accessed using the access operators or the at() methods
  79 + *
  80 + * double element = matrix(row, column);
  81 + * double element = matrix.at(row, column);
  82 + *
  83 + * the difference being that the access operators do NOT but the at()
  84 + * methods do perform validity checking of the @p row and @p column
  85 + * arguments. Access operators and at() methods exist for non-const and
  86 + * const references.
  87 + *
  88 + * For more efficient matrix access, methods operating on matrix rows and
  89 + * columns are available. The
  90 + *
  91 + * GVector row = matrix.row(index);
  92 + * GVector column = matrix.column(index);
  93 + *
  94 + * methods return a row or a column in a GVector object, while the
  95 + *
  96 + * GVector row;
  97 + * GVector column;
  98 + * ...
  99 + * matrix.row(index, row);
  100 + * matrix.column(index, column);
  101 + *
  102 + * store a row or a column provided by a GVector object into the matrix.
  103 + * Obviously, the vector needs to be of the same size as the number of rows
  104 + * or columns available in the matrix. Adding of matrix elements row-by-row
  105 + * or column-by-column is done using the
  106 + *
  107 + * matrix.add_to_row(index, row);
  108 + * matrix.add_to_column(index, column);
  109 + *
  110 + * methods.
  111 + *
  112 + * Matrix multiplication with a vector is done using
  113 + *
  114 + * GVector vector;
  115 + * ...
  116 + * GVector product = matrix * vector;
  117 + *
  118 + * where the size of @p vector needs to be identical to the number of
  119 + * columns in the matrix.
  120 + *
  121 + * Furthermore, the following methods exist for matrix operations:
  122 + *
  123 + * matrix.transpose(); // Transposes matrix
  124 + * matrix.invert(); // Inverts matrix
  125 + * matrix.negate(); // Negates all matrix elements
  126 + * matrix.abs(); // Sets all matrix elements to their absolute values
  127 + *
  128 + * Friend functions are implemented to perform the same operations also in a
  129 + * more natural syntax that will generate a new matrix
  130 + *
  131 + * transposed_matrix = transpose(matrix);
  132 + * inverted_matrix = invert(matrix);
  133 + * negated_matrix = -matrix;
  134 + * absolute_matrix = abs(matrix);
  135 + *
  136 + * The following methods allow to access matrix attributes:
  137 + *
  138 + * matrix.fill(); // Percentage of non-zero elements in matrix [0,1]
  139 + * matrix.min(); // Minimum matrix element
  140 + * matrix.max(); // Maximum matrix element
  141 + * matrix.sum(); // Sum of all matrix elements
  142 + * matrix.rows(); // Number of rows in matrix
  143 + * matrix.columns(); // Number of columns in matrix
  144 + * matrix.size(); // Number of elements in matrix
44 145 ***************************************************************************/
45 146 class GMatrixBase : public GBase {
46 147  
... ... @@ -61,23 +162,25 @@ public:
61 162 virtual bool operator!=(const GMatrixBase& matrix) const;
62 163  
63 164 // Pure virtual methods
64   - virtual void clear(void) = 0;
65   - virtual GMatrixBase* clone(void) const = 0;
66   - virtual GVector row(const int& row) const = 0;
67   - virtual void row(const int& row, const GVector& vector) = 0;
68   - virtual GVector column(const int& column) const = 0;
69   - virtual void column(const int& column, const GVector& vector) = 0;
70   - virtual void add_to_row(const int& row, const GVector& vector) = 0;
71   - virtual void add_to_column(const int& column, const GVector& vector) = 0;
72   - virtual void transpose(void) = 0;
73   - virtual void invert(void) = 0;
74   - virtual void negate(void) = 0;
75   - virtual void abs(void) = 0;
76   - virtual double fill(void) const = 0;
77   - virtual double min(void) const = 0;
78   - virtual double max(void) const = 0;
79   - virtual double sum(void) const = 0;
80   - virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
  165 + virtual void clear(void) = 0;
  166 + virtual GMatrixBase* clone(void) const = 0;
  167 + virtual double& at(const int& row, const int& column) = 0;
  168 + virtual const double& at(const int& row, const int& column) const = 0;
  169 + virtual GVector row(const int& row) const = 0;
  170 + virtual void row(const int& row, const GVector& vector) = 0;
  171 + virtual GVector column(const int& column) const = 0;
  172 + virtual void column(const int& column, const GVector& vector) = 0;
  173 + virtual void add_to_row(const int& row, const GVector& vector) = 0;
  174 + virtual void add_to_column(const int& column, const GVector& vector) = 0;
  175 + virtual void transpose(void) = 0;
  176 + virtual void invert(void) = 0;
  177 + virtual void negate(void) = 0;
  178 + virtual void abs(void) = 0;
  179 + virtual double fill(void) const = 0;
  180 + virtual double min(void) const = 0;
  181 + virtual double max(void) const = 0;
  182 + virtual double sum(void) const = 0;
  183 + virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
81 184  
82 185 // Base class methods
83 186 const int& size(void) const;
... ...
include/GMatrixSparse.hpp
  View file @3050674
... ... @@ -49,16 +49,20 @@ class GSparseNumeric;
49 49 *
50 50 * @brief Sparse matrix class interface defintion
51 51 *
52   - * This class implements a sparse matrix class. The class only stores
53   - * non-zero elements, which can considerably reduce the memory requirements
54   - * for large systems that are sparsly filled.
  52 + * This class implements a sparse matrix. The class only stores non-zero
  53 + * elements, which can considerably reduce the memory requirements for large
  54 + * systems that are sparsely filled. Also the computations will be faster
  55 + * because only non-zero elements will be used in the operations.
  56 + *
  57 + * For a description of the common matrix methods, please refer to the
  58 + * GMatrixBase class.
55 59 *
56 60 * The class has been inspired from the code CSparse that can be downloaded
57   - * from http://www.cise.ufl.edu/research/sparse/CSparse/
58   - * CSparse has been written by Timothy A. Davis to accompany his book
59   - * "Direct Methods for Sparse Linear Systems". In the CSparse code, a sparse
60   - * matrix is implemented by the structure cs_sparse. The members of this
61   - * structure map as follows to the member of the GMatrixSparse class:
  61 + * from http://www.cise.ufl.edu/research/sparse/CSparse/ which has been
  62 + * written by Timothy A. Davis to accompany his book "Direct Methods for
  63 + * Sparse Linear Systems". In the CSparse code, a sparse matrix is
  64 + * implemented by the structure cs_sparse. The members of this
  65 + * structure map as follows to the members of the GMatrixSparse class:
62 66 *
63 67 * <table border>
64 68 * <tr>
... ... @@ -106,6 +110,76 @@ class GSparseNumeric;
106 110 * </tr>
107 111 * </table>
108 112 *
  113 + * Matrix allocation is done using the constructors
  114 + *
  115 + * GMatrixSparse sparsematrix(rows, columns);
  116 + * GMatrixSparse sparsematrix(rows, columns, elements);
  117 + * GMatrixSparse sparsematrix(matrix);
  118 + * GMatrixSparse sparsematrix(sparsematrix);
  119 + * GMatrixSparse sparsematrix(symmetricmatrix);
  120 + *
  121 + * where @p rows and @p columns specify the number of rows and columns of
  122 + * the matrix and @p elements is the size of the memory that should be
  123 + * pre-allocated for matrix storage (optional). Storage conversion
  124 + * constructors exist that allow allocating a sparse matrix by copying from
  125 + * a general matrix of type GMatrix and a symmetrix matrix of type
  126 + * GMatrixSymmetric.
  127 + *
  128 + * Because the location of a specific matrix element is not necessarily known
  129 + * in advance, the class implements a so called pending element for matrix
  130 + * element setting. A pending element is an element that is scheduled for
  131 + * insertion in the matrix, but which has not yet been inserted. The matrix
  132 + * element access operators () and at() return the reference to this pending
  133 + * element if a memory location for a given row and column does not yet
  134 + * exist. The protected method fill_pending() is then called before all
  135 + * operations that are done on the matrix to insert the element before the
  136 + * operation is executed.
  137 + *
  138 + * Another mechanism that has been implemented to speed up the fill of the
  139 + * sparse matrix is a "fill stack" that is used for insertion or addition
  140 + * of matrix columns by the
  141 + *
  142 + * matrix.column(index, vector);
  143 + * matrix.column(index, values, rows, number);
  144 + * matrix.add_to_column(index, vector);
  145 + * matrix.add_to_column(index, values, rows, number);
  146 + *
  147 + * methods. The fill stack is a buffer that implements a queue for columns
  148 + * that are to be inserted into the matrix. This is more efficient than
  149 + * insertion of elements one by one, as every insertion of an element may
  150 + * require to shift all elements back by one memory location.
  151 + *
  152 + * The fill stack is used as follows:
  153 + *
  154 + * matrix.stack_init(size, entries);
  155 + * ...
  156 + * matrix.column(index, vector); // (or one of the other column methods)
  157 + * ...
  158 + * matrix.stack_flush();
  159 + * ...
  160 + * matrix.stack_destroy();
  161 + *
  162 + * The stack_init(size, entries) method initialises the fill stack, where
  163 + * @p size is the size of the allocated memory buffer and @p entries is the
  164 + * maximum number of entries (i.e. columns) that will be held by the buffer.
  165 + * If @p size is set to 0 (the default value), a default @p size value of
  166 + * 512 is used. If @p entries is set to 0 (the default value), the number of
  167 + * matrix columns is taken as default @p entries value. Note that a too large
  168 + * number of elements will produce some overhead due to fill stack
  169 + * management, hence @p entries should not exceeed a value of the order of
  170 + * 10-100.
  171 + *
  172 + * The stack_flush() method flushes the stack, which is mandatory before any
  173 + * usage of the matrix. Note that the fill stack IS NOT INSERTED AUTOMATICALLY
  174 + * before any matrix operation, hence manual stack flushing is needed to make
  175 + * all filled matrix elements available for usage. The stack_destroy() method
  176 + * will flush the stack and free all stack elements. This method
  177 + * should be called once no filling is required anymore. If stack_destroy()
  178 + * is called immediately after filling, no call to stack_flush() is needed as
  179 + * the stack_destroy() method flushes the stack before destroying it. The
  180 + * matrix stack is also destroyed by the sparse matrix destructor, hence
  181 + * manual stack destruction is not mandatory.
  182 + *
109 183 * Except for *m_rowinx which is implemented on the level of GMatrixSparse,
110 184 * all other members are implemented by the base class GMatrixBase.
111 185 ***************************************************************************/
... ... @@ -144,6 +218,7 @@ public:
144 218  
145 219 // Other operators
146 220 virtual GMatrixSparse& operator=(const GMatrixSparse& matrix);
  221 + virtual GMatrixSparse& operator=(const double& value);
147 222 virtual GMatrixSparse operator+(const GMatrixSparse& matrix) const;
148 223 virtual GMatrixSparse operator-(const GMatrixSparse& matrix) const;
149 224 virtual GMatrixSparse operator*(const GMatrixSparse& matrix) const;
... ... @@ -157,6 +232,8 @@ public:
157 232 // Implemented pure virtual base class methods
158 233 virtual void clear(void);
159 234 virtual GMatrixSparse* clone(void) const;
  235 + virtual double& at(const int& row, const int& column);
  236 + virtual const double& at(const int& row, const int& column) const;
160 237 virtual GVector row(const int& row) const;
161 238 virtual void row(const int& row, const GVector& vector);
162 239 virtual GVector column(const int& column) const;
... ... @@ -174,10 +251,10 @@ public:
174 251 virtual std::string print(const GChatter& chatter = NORMAL) const;
175 252  
176 253 // Other methods
177   - void column(const double* values, const int* rows,
178   - int number, const int& col);
179   - void add_to_column(const double* values, const int* rows,
180   - int number, const int& col);
  254 + void column(const int& column, const double* values,
  255 + const int* rows, int number);
  256 + void add_to_column(const int& column, const double* values,
  257 + const int* rows, int number);
181 258 void cholesky_decompose(bool compress = true);
182 259 GVector cholesky_solver(const GVector& vector, bool compress = true);
183 260 void cholesky_invert(bool compress = true);
... ...
pyext/GMatrixBase.i
  View file @3050674
... ... @@ -42,29 +42,31 @@ public:
42 42 virtual ~GMatrixBase(void);
43 43  
44 44 // Pure virtual operators
45   - virtual GVector operator*(const GVector& vector) const = 0;
  45 + virtual GVector operator*(const GVector& vector) const = 0;
46 46  
47 47 // Implemented base class operators
48   - virtual bool operator==(const GMatrixBase& matrix) const;
49   - virtual bool operator!=(const GMatrixBase& matrix) const;
  48 + virtual bool operator==(const GMatrixBase& matrix) const;
  49 + virtual bool operator!=(const GMatrixBase& matrix) const;
50 50  
51 51 // Pure virtual methods
52   - virtual void clear(void) = 0;
53   - virtual GMatrixBase* clone(void) const = 0;
54   - virtual GVector row(const int& row) const = 0;
55   - virtual void row(const int& row, const GVector& vector) = 0;
56   - virtual GVector column(const int& column) const = 0;
57   - virtual void column(const int& column, const GVector& vector) = 0;
58   - virtual void add_to_row(const int& row, const GVector& vector) = 0;
59   - virtual void add_to_column(const int& column, const GVector& vector) = 0;
60   - virtual void transpose(void) = 0;
61   - virtual void invert(void) = 0;
62   - virtual void negate(void) = 0;
63   - virtual void abs(void) = 0;
64   - virtual double fill(void) const = 0;
65   - virtual double min(void) const = 0;
66   - virtual double max(void) const = 0;
67   - virtual double sum(void) const = 0;
  52 + virtual void clear(void) = 0;
  53 + virtual GMatrixBase* clone(void) const = 0;
  54 + virtual double& at(const int& row, const int& column) = 0;
  55 + virtual const double& at(const int& row, const int& column) const = 0;
  56 + virtual GVector row(const int& row) const = 0;
  57 + virtual void row(const int& row, const GVector& vector) = 0;
  58 + virtual GVector column(const int& column) const = 0;
  59 + virtual void column(const int& column, const GVector& vector) = 0;
  60 + virtual void add_to_row(const int& row, const GVector& vector) = 0;
  61 + virtual void add_to_column(const int& column, const GVector& vector) = 0;
  62 + virtual void transpose(void) = 0;
  63 + virtual void invert(void) = 0;
  64 + virtual void negate(void) = 0;
  65 + virtual void abs(void) = 0;
  66 + virtual double fill(void) const = 0;
  67 + virtual double min(void) const = 0;
  68 + virtual double max(void) const = 0;
  69 + virtual double sum(void) const = 0;
68 70  
69 71 // Base class methods
70 72 const int& size(void) const;
... ...
pyext/GMatrixSparse.i
  View file @3050674
... ... @@ -73,6 +73,8 @@ public:
73 73 // Implemented pure virtual base class methods
74 74 virtual void clear(void);
75 75 virtual GMatrixSparse* clone(void) const;
  76 + virtual double& at(const int& row, const int& column);
  77 + virtual const double& at(const int& row, const int& column) const;
76 78 virtual GVector row(const int& row) const;
77 79 virtual void row(const int& row, const GVector& vector);
78 80 virtual GVector column(const int& column) const;
... ... @@ -89,10 +91,10 @@ public:
89 91 virtual double sum(void) const;
90 92  
91 93 // Other methods
92   - void add_to_column(const double* values, const int* rows,
93   - int number, const int& column);
94   - void column(const double* values, const int* rows,
95   - int number, const int& column);
  94 + void column(const int& column, const double* values,
  95 + const int* rows, int number);
  96 + void add_to_column(const int& column, const double* values,
  97 + const int* rows, int number);
96 98 void cholesky_decompose(bool compress = true);
97 99 GVector cholesky_solver(const GVector& vector, bool compress = true);
98 100 void cholesky_invert(bool compress = true);
... ... @@ -125,4 +127,8 @@ public:
125 127 GMatrixSparse copy() {
126 128 return (*self);
127 129 }
  130 + GMatrixSparse set(const double& value) {
  131 + (*self) = value;
  132 + return (*self);
  133 + }
128 134 };
... ...
src/linalg/GMatrixSparse.cpp
  View file @3050674
... ... @@ -17,19 +17,6 @@
17 17 * You should have received a copy of the GNU General Public License *
18 18 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19 19 * *
20   - * ----------------------------------------------------------------------- *
21   - * This class implements the compressed sparse column format. The *
22   - * following arrays are allocated: *
23   - * *
24   - * m_rows (n) Number of rows *
25   - * m_cols (m) Number of columns *
26   - * m_data (Ax) Holds all 'elements' non-zero values, in column order *
27   - * m_colstart (Ap) Holds the index of the first element of each column *
28   - * m_rowinx (Ai) Holds the row indices for all elements *
29   - * *
30   - * Column 'col' covers therefore [m_colstart[j], ..., m_colstart[j+1]-1] *
31   - * Note that 'm_colstart' has m_cols+1 elements. *
32   - * *
33 20 ***************************************************************************/
34 21 /**
35 22 * @file GMatrixSparse.cpp
... ... @@ -53,20 +40,20 @@
53 40  
54 41 /* __ Method name definitions ____________________________________________ */
55 42 #define G_CONSTRUCTOR "GMatrixSparse::GMatrixSparse(int&, int&, int&)"
56   -#define G_ACCESS "GMatrixSparse::operator(int&, int&)"
57 43 #define G_OP_MUL_VEC "GMatrixSparse::operator*(GVector&)"
58 44 #define G_OP_ADD "GMatrixSparse::operator+=(GMatrixSparse&)"
59 45 #define G_OP_SUB "GMatrixSparse::operator-=(GMatrixSparse&)"
60 46 #define G_OP_MAT_MUL "GMatrixSparse::operator*=(GMatrixSparse&)"
  47 +#define G_AT "GMatrixSparse::at(int&, int&)"
61 48 #define G_EXTRACT_ROW "GMatrixSymmetric::row(int&)"
62 49 #define G_SET_ROW "GMatrixSymmetric::row(int&, GVector&)"
63 50 #define G_EXTRACT_COLUMN "GMatrixSymmetric::column(int&)"
64 51 #define G_SET_COLUMN "GMatrixSymmetric::column(int&, GVector&)"
65   -#define G_SET_COLUMN2 "GMatrixSymmetric::column(double*, int*, int, int&)"
  52 +#define G_SET_COLUMN2 "GMatrixSymmetric::column(int&, double*, int*, int)"
66 53 #define G_ADD_TO_ROW "GMatrixSymmetric::add_to_row(int&, GVector&)"
67 54 #define G_ADD_TO_COLUMN "GMatrixSymmetric::add_to_column(int&, GVector&)"
68   -#define G_ADD_TO_COLUMN2 "GMatrixSymmetric::add_to_column(double*, int*,"\
69   - " int, int&)"
  55 +#define G_ADD_TO_COLUMN2 "GMatrixSymmetric::add_to_column(int&, double*,"
  56 + " int*, int)"
70 57 #define G_INVERT "GMatrixSparse::invert(void)"
71 58 #define G_CHOL_DECOMP "GMatrixSparse::cholesky_decompose(bool)"
72 59 #define G_CHOL_SOLVE "GMatrixSparse::cholesky_solver(GVector&, bool)"
... ... @@ -137,6 +124,9 @@ GMatrixSparse::GMatrixSparse(void) : GMatrixBase()
137 124 * of rows and columns. The elements parameter allows to specify how much
138 125 * physical memory should be allocated initially. By default, no memory
139 126 * will be allocated.
  127 + *
  128 + * @todo Is there a real need to throw an empty exception? The class should
  129 + * be able to operate on empty matrices.
140 130 ***************************************************************************/
141 131 GMatrixSparse::GMatrixSparse(const int& rows,
142 132 const int& columns,
... ... @@ -257,9 +247,12 @@ GMatrixSparse::~GMatrixSparse(void)
257 247 ==========================================================================*/
258 248  
259 249 /***********************************************************************//**
260   - * @brief Assignment operator
  250 + * @brief Matrix assignment operator
261 251 *
262 252 * @param[in] matrix Matrix.
  253 + * @return Matrix.
  254 + *
  255 + * Assigns the content of another matrix to the actual matrix instance.
263 256 ***************************************************************************/
264 257 GMatrixSparse& GMatrixSparse::operator=(const GMatrixSparse& matrix)
265 258 {
... ... @@ -288,25 +281,65 @@ GMatrixSparse&amp; GMatrixSparse::operator=(const GMatrixSparse&amp; matrix)
288 281  
289 282  
290 283 /***********************************************************************//**
291   - * @brief Access operator
  284 + * @brief Value assignment operator
  285 + *
  286 + * @param[in] value Value.
  287 + * @return Matrix.
  288 + *
  289 + * Assigns the specified @p value to all elements of the matrix.
  290 + ***************************************************************************/
  291 +GMatrixSparse& GMatrixSparse::operator=(const double& value)
  292 +{
  293 + // Fill any pending element to have a non-pending state
  294 + fill_pending();
  295 +
  296 + // If value is 0 then simply reinitialize column start indices
  297 + if (value == 0) {
  298 +
  299 + // Initialise column start indices to 0
  300 + for (int col = 0; col <= m_cols; ++col) {
  301 + m_colstart[col] = 0;
  302 + }
  303 +
  304 + }
  305 +
  306 + // ... otherwise fill column-wise
  307 + else {
  308 +
  309 + // Set column vector
  310 + GVector column(m_rows);
  311 + column = value;
  312 +
  313 + // Column-wise setting
  314 + for (int col = 0; col < m_cols; ++col) {
  315 + this->column(col, column);
  316 + }
  317 +
  318 + }
  319 +
  320 + // Return this object
  321 + return *this;
  322 +}
  323 +
  324 +
  325 +/***********************************************************************//**
  326 + * @brief Return reference to matrix element
292 327 *
293 328 * @param[in] row Matrix row [0,...,rows()-1].
294 329 * @param[in] column Matrix column [0,...,columns()-1].
  330 + * @return Reference to matrix element.
295 331 *
296   - * @exception GException::out_of_range
297   - * Row or column index out of range.
  332 + * Returns the reference to the matrix element at @p row and @p column. If
  333 + * the matrix element does not yet exist, a reference to the pending element
  334 + * with a value of 0.0 is returned.
298 335 ***************************************************************************/
299 336 double& GMatrixSparse::operator()(const int& row, const int& column)
300 337 {
301   - // Compile option: perform range check
302   - #if defined(G_RANGE_CHECK)
303   - if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) {
304   - throw GException::out_of_range(G_ACCESS, row, column, m_rows, m_cols);
305   - }
306   - #endif
  338 + // Fill pending element. This will set the value of the pending element
  339 + // to 0.0.
  340 + fill_pending();
307 341  
308 342 // Get element
309   - fill_pending();
310 343 int inx = get_index(row, column);
311 344 double* value;
312 345 if (inx < 0) {
... ... @@ -324,29 +357,24 @@ double&amp; GMatrixSparse::operator()(const int&amp; row, const int&amp; column)
324 357  
325 358  
326 359 /***********************************************************************//**
327   - * @brief Access operator (const version)
  360 + * @brief Return reference to matrix element (const version)
328 361 *
329 362 * @param[in] row Matrix row [0,...,rows()-1].
330 363 * @param[in] column Matrix column [0,...,columns()-1].
  364 + * @return Const reference to matrix element.
331 365 *
332   - * @exception GException::out_of_range
333   - * Row or column index out of range.
334   - *
335   - * We need here the zero element to return also a pointer for 0.0 entry that
336   - * is not stored. Since we have the const version we don't have to care about
337   - * modification of this zero value.
  366 + * Returns a const reference to the matrix element at @p row and @p column.
  367 + * If the matrix element does not yet exist, a reference to the zero element
  368 + * is returned. If the matrix element corresponds to the pending element,
  369 + * a reference to the pending element is returned. Otherwise, a reference
  370 + * to the matrix elements is returned.
338 371 ***************************************************************************/
339 372 const double& GMatrixSparse::operator()(const int& row,
340 373 const int& column) const
341 374 {
342   - // Compile option: perform range check
343   - #if defined(G_RANGE_CHECK)
344   - if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) {
345   - throw GException::out_of_range(G_ACCESS, row, column, m_rows, m_cols);
346   - }
347   - #endif
348   -
349   - // Get element
  375 + // Get element. We need here the zero element to return also a pointer
  376 + // for 0.0 entry that is not stored. Since we have the const version we
  377 + // don't have to care about modification of this zero value.
350 378 int inx = get_index(row, column);
351 379 double* value;
352 380 if (inx < 0) {
... ... @@ -635,6 +663,50 @@ GMatrixSparse* GMatrixSparse::clone(void) const
635 663  
636 664  
637 665 /***********************************************************************//**
  666 + * @brief Return reference to matrix element
  667 + *
  668 + * @param[in] row Matrix row [0,...,rows()-1].
  669 + * @param[in] column Matrix column [0,...,columns()-1].
  670 + * @return Reference to matrix element.
  671 + *
  672 + * @exception GException::out_of_range
  673 + * Row or column index out of range.
  674 + ***************************************************************************/
  675 +double& GMatrixSparse::at(const int& row, const int& column)
  676 +{
  677 + // Raise exception if row or column index is out of range
  678 + if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) {
  679 + throw GException::out_of_range(G_AT, row, column, m_rows, m_cols);
  680 + }
  681 +
  682 + // Return element
  683 + return ((*this)(row, column));
  684 +}
  685 +
  686 +
  687 +/***********************************************************************//**
  688 + * @brief Return reference to matrix element (const version)
  689 + *
  690 + * @param[in] row Matrix row [0,...,rows()-1].
  691 + * @param[in] column Matrix column [0,...,columns()-1].
  692 + * @return Reference to matrix element.
  693 + *
  694 + * @exception GException::out_of_range
  695 + * Row or column index out of range.
  696 + ***************************************************************************/
  697 +const double& GMatrixSparse::at(const int& row, const int& column) const
  698 +{
  699 + // Raise exception if row or column index is out of range
  700 + if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) {
  701 + throw GException::out_of_range(G_AT, row, column, m_rows, m_cols);
  702 + }
  703 +
  704 + // Return element
  705 + return ((*this)(row, column));
  706 +}
  707 +
  708 +
  709 +/***********************************************************************//**
638 710 * @brief Extract row as vector from matrix
639 711 *
640 712 * @param[in] row Row to be extracted (starting from 0).
... ... @@ -888,10 +960,10 @@ void GMatrixSparse::column(const int&amp; column, const GVector&amp; vector)
888 960 /***********************************************************************//**
889 961 * @brief Insert compressed array into matrix column
890 962 *
  963 + * @param[in] column Column index [0,...,columns()-1].
891 964 * @param[in] values Compressed array.
892 965 * @param[in] rows Row indices of array.
893 966 * @param[in] number Number of elements in array.
894   - * @param[in] column Column index [0,...,columns()-1].
895 967 *
896 968 * @exception GException::out_of_range
897 969 * Invalid column index specified.
... ... @@ -904,12 +976,13 @@ void GMatrixSparse::column(const int&amp; column, const GVector&amp; vector)
904 976 * This is the main driver routine to insert data into a matrix. Note that
905 977 * there is another instance of this function that takes a vector.
906 978 ***************************************************************************/
907   -void GMatrixSparse::column(const double* values, const int* rows,
908   - int number, const int& column)
  979 +void GMatrixSparse::column(const int& column, const double* values,
  980 + const int* rows, int number)
909 981 {
910 982 // Debug header
911 983 #if defined(G_DEBUG_SPARSE_INSERTION)
912   - std::cout << "GMatrixSparse::insert_col(v, i, n, " << column << "):" << std::endl;
  984 + std::cout << "GMatrixSparse::column(";
  985 + std::cout << column << ", values, rows, " << number << "):" << std::endl;
913 986 std::cout << " Matrix Data : ";
914 987 for (int i = 0; i < m_elements; ++i) {
915 988 std::cout << m_data[i] << " ";
... ... @@ -1166,10 +1239,10 @@ void GMatrixSparse::add_to_column(const int&amp; column, const GVector&amp; vector)
1166 1239 /***********************************************************************//**
1167 1240 * @brief Add compressed array into matrix column
1168 1241 *
  1242 + * @param[in] column Column index [0,...,columns()-1].
1169 1243 * @param[in] values Compressed array.
1170 1244 * @param[in] rows Row indices of array.
1171 1245 * @param[in] number Number of elements in array.
1172   - * @param[in] col Column index (starting from 0).
1173 1246 *
1174 1247 * @exception GException::out_of_range
1175 1248 * Invalid column index specified.
... ... @@ -1182,12 +1255,13 @@ void GMatrixSparse::add_to_column(const int&amp; column, const GVector&amp; vector)
1182 1255 * normal and stack-based filled. Note that there is another instance of this
1183 1256 * method that takes a vector.
1184 1257 ***************************************************************************/
1185   -void GMatrixSparse::add_to_column(const double* values, const int* rows,
1186   - int number, const int& col)
  1258 +void GMatrixSparse::add_to_column(const int& column, const double* values,
  1259 + const int* rows, int number)
1187 1260 {
1188 1261 // Debug header
1189 1262 #if defined(G_DEBUG_SPARSE_ADDITION)
1190   - std::cout << "GMatrixSparse::add_col(v, i, n, " << col << "):" << std::endl;
  1263 + std::cout << "GMatrixSparse::add_col(";
  1264 + std::cout << column << ", values, rows, " << number << "):" << std::endl;
1191 1265 std::cout << " Matrix Data : ";
1192 1266 for (int i = 0; i < m_elements; ++i) {
1193 1267 std::cout << m_data[i] << " ";
... ... @@ -1216,7 +1290,7 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows,
1216 1290 // stack_push_column does its own argument verifications, so to avoid
1217 1291 // double checking we don't do anything before this call ...
1218 1292 if (m_stack_data != NULL) {
1219   - number = stack_push_column(values, rows, number, col);
  1293 + number = stack_push_column(values, rows, number, column);
1220 1294 if (number == 0) {
1221 1295 return;
1222 1296 }
... ... @@ -1231,8 +1305,8 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows,
1231 1305  
1232 1306 // Raise an exception if the column index is invalid
1233 1307 #if defined(G_RANGE_CHECK)
1234   - if (col < 0 || col >= m_cols) {
1235   - throw GException::out_of_range(G_ADD_TO_COLUMN2, col, 0, m_cols-1);
  1308 + if (column < 0 || column >= m_cols) {
  1309 + throw GException::out_of_range(G_ADD_TO_COLUMN2, column, 0, m_cols-1);
1236 1310 }
1237 1311 #endif
1238 1312  
... ... @@ -1245,8 +1319,8 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows,
1245 1319 } // endelse: there was no stack
1246 1320  
1247 1321 // Get indices of column in matrix
1248   - int i_start = m_colstart[col];
1249   - int i_stop = m_colstart[col+1];
  1322 + int i_start = m_colstart[column];
  1323 + int i_stop = m_colstart[column+1];
1250 1324  
1251 1325 // Case A: the column exists in the matrix, so mix new elements with existing
1252 1326 // data
... ... @@ -1267,7 +1341,7 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows,
1267 1341 wrk_double, wrk_int, &num_mix);
1268 1342  
1269 1343 // Insert mixed column
1270   - this->column(wrk_double, wrk_int, num_mix, col);
  1344 + this->column(column, wrk_double, wrk_int, num_mix);
1271 1345  
1272 1346 // Free workspace
1273 1347 delete [] wrk_int;
... ... @@ -1277,20 +1351,23 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows,
1277 1351  
1278 1352 // Case B: the column does not yet exist in the matrix, so just insert it
1279 1353 else {
1280   - this->column(values, rows, number, col);
  1354 + this->column(column, values, rows, number);
1281 1355 }
1282 1356  
1283 1357 // Debugging: show sparse matrix after insertion
1284 1358 #if defined(G_DEBUG_SPARSE_ADDITION)
1285 1359 std::cout << " Out Data: ";
1286   - for (int i = 0; i < m_elements; ++i)
  1360 + for (int i = 0; i < m_elements; ++i) {
1287 1361 std::cout << m_data[i] << " ";
  1362 + }
1288 1363 std::cout << std::endl << " Out Row : ";
1289   - for (int i = 0; i < m_elements; ++i)
  1364 + for (int i = 0; i < m_elements; ++i) {
1290 1365 std::cout << m_rowinx[i] << " ";
  1366 + }
1291 1367 std::cout << std::endl << " Out Col : ";
1292   - for (int i = 0; i < m_cols+1; ++i)
  1368 + for (int i = 0; i < m_cols+1; ++i) {
1293 1369 std::cout << m_colstart[i] << " ";
  1370 + }
1294 1371 std::cout << std::endl;
1295 1372 #endif
1296 1373  
... ... @@ -1321,15 +1398,18 @@ void GMatrixSparse::transpose(void)
1321 1398 /***********************************************************************//**
1322 1399 * @brief Invert matrix
1323 1400 *
1324   - * @exception GException::feature_not_implemented
1325   - * Feature not yet implemented.
  1401 + * Inverts the matrix using the Cholesky decomposition. This does not work
  1402 + * on any kind of matrix.
1326 1403 *
1327   - * @todo Needs to be implemented.
  1404 + * @todo Specify in documentation for which kind of matrix the method works.
1328 1405 ***************************************************************************/
1329 1406 void GMatrixSparse::invert(void)
1330 1407 {
1331   - // Throw exception
1332   - throw GException::feature_not_implemented(G_INVERT);
  1408 + // Perform Chloesky decomposition
  1409 + cholesky_decompose(true);
  1410 +
  1411 + // Invert matrix
  1412 + cholesky_invert(true);
1333 1413  
1334 1414 // Return
1335 1415 return;
... ...