Justine Devin / gammalib

Commit 5b16109ee05eb1e6c4fccf7a231fd12688d66dde

Authored by Jürgen Knödlseder
1 parent ef62b5af

Rework matrix transformation methods. 

The methods transpose(), invert() and abs() now return a matrix instead of working on the actual matrix. The transpose(), inverse(), abs(), and negate() friend functions have been removed. The method solve() has been added that solves the linear matrix equation.

The Cholesky methods now also return a matrix instead of working on the actual object.
Inline Side-by-side
Showing 11 changed files with 485 additions and 622 deletions
include/GMatrix.hpp
  View file @5b16109
... ... @@ -34,6 +34,7 @@
34 34 /* __ Forward declarations _______________________________________________ */
35 35 class GMatrixSymmetric;
36 36 class GMatrixSparse;
  37 +class GVector;
37 38  
38 39  
39 40 /***********************************************************************//**
... ... @@ -115,10 +116,6 @@ public:
115 116 virtual void column(const int& column, const GVector& vector);
116 117 virtual void add_to_row(const int& row, const GVector& vector);
117 118 virtual void add_to_column(const int& column, const GVector& vector);
118   - virtual void transpose(void);
119   - virtual void invert(void);
120   - virtual void negate(void);
121   - virtual void abs(void);
122 119 virtual double fill(void) const;
123 120 virtual double min(void) const;
124 121 virtual double max(void) const;
... ... @@ -126,11 +123,16 @@ public:
126 123 virtual std::string print(const GChatter& chatter = NORMAL) const;
127 124  
128 125 // Other methods
129   - virtual GMatrix extract_lower_triangle(void) const;
130   - virtual GMatrix extract_upper_triangle(void) const;
131   - virtual void eulerx(const double& angle);
132   - virtual void eulery(const double& angle);
133   - virtual void eulerz(const double& angle);
  126 + GMatrix transpose(void) const;
  127 + GMatrix invert(void) const;
  128 + GVector solve(const GVector& vector) const;
  129 + GMatrix negate(void) const;
  130 + GMatrix abs(void) const;
  131 + GMatrix extract_lower_triangle(void) const;
  132 + GMatrix extract_upper_triangle(void) const;
  133 + void eulerx(const double& angle);
  134 + void eulery(const double& angle);
  135 + void eulerz(const double& angle);
134 136  
135 137 private:
136 138 // Private methods
... ... @@ -262,23 +264,6 @@ GMatrix& GMatrix::operator/=(const double& scalar)
262 264  
263 265  
264 266 /***********************************************************************//**
265   - * @brief Negate matrix elements
266   - *
267   - * @return Matrix with negated elements.
268   - *
269   - * Returns a matrix where each element has been replaced by its negative
270   - * element.
271   - ***************************************************************************/
272   -inline
273   -GMatrix GMatrix::operator-(void) const
274   -{
275   - GMatrix result = *this;
276   - result.negate();
277   - return result;
278   -}
279   -
280   -
281   -/***********************************************************************//**
282 267 * @brief Return minimum matrix element
283 268 *
284 269 * @return Minimum element in matrix.
... ... @@ -373,57 +358,4 @@ GMatrix operator/(const GMatrix& matrix, const double& scalar)
373 358 return result;
374 359 }
375 360  
376   -
377   -/***********************************************************************//**
378   - * @brief Return transpose of matrix
379   - *
380   - * @param[in] matrix Matrix.
381   - * @return Transpose of matrix.
382   - *
383   - * Returns the transpose of a matrix.
384   - ***************************************************************************/
385   -inline
386   -GMatrix transpose(const GMatrix& matrix)
387   -{
388   - GMatrix result = matrix;
389   - result.transpose();
390   - return result;
391   -}
392   -
393   -
394   -/***********************************************************************//**
395   - * @brief Return inverse of matrix
396   - *
397   - * @param[in] matrix Matrix.
398   - * @return Inverse of matrix.
399   - *
400   - * Returns the inverse of a matrix.
401   - ***************************************************************************/
402   -inline
403   -GMatrix invert(const GMatrix& matrix)
404   -{
405   - GMatrix result = matrix;
406   - result.invert();
407   - return result;
408   -}
409   -
410   -
411   -/***********************************************************************//**
412   - * @brief Return matrix with absolute values of all elements
413   - *
414   - * @param[in] matrix Matrix.
415   - * @return Matrix with elements being the absolute elements of the input
416   - * matrix.
417   - *
418   - * Return a matrix where each elements has been replaced by its absolute
419   - * value.
420   - ***************************************************************************/
421   -inline
422   -GMatrix abs(const GMatrix& matrix)
423   -{
424   - GMatrix result = matrix;
425   - result.abs();
426   - return result;
427   -}
428   -
429 361 #endif /* GMATRIX_HPP */
... ...
include/GMatrixBase.hpp
  View file @5b16109
... ... @@ -118,20 +118,12 @@
118 118 * where the size of @p vector needs to be identical to the number of
119 119 * columns in the matrix.
120 120 *
121   - * Furthermore, the following methods exist for matrix operations:
  121 + * Furthermore, the following methods exist for matrix transformations:
122 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);
  123 + * new_matrix = matrix.transpose(); // Transposes matrix
  124 + * new_matrix = matrix.invert(); // Inverts matrix
  125 + * new_vector = matrix.solve(vector); // Solves linear matrix equation
  126 + * new_matrix = matrix.abs(); // Returns matrix with all elements replaced by their absolute values
135 127 *
136 128 * The following methods allow to access matrix attributes:
137 129 *
... ... @@ -172,10 +164,6 @@ public:
172 164 virtual void column(const int& column, const GVector& vector) = 0;
173 165 virtual void add_to_row(const int& row, const GVector& vector) = 0;
174 166 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 167 virtual double fill(void) const = 0;
180 168 virtual double min(void) const = 0;
181 169 virtual double max(void) const = 0;
... ...
include/GMatrixSparse.hpp
  View file @5b16109
... ... @@ -40,6 +40,7 @@
40 40 class GMatrix;
41 41 class GMatrixSymmetric;
42 42 class GMatrixSparse;
  43 +class GVector;
43 44 class GSparseSymbolic;
44 45 class GSparseNumeric;
45 46  
... ... @@ -240,10 +241,6 @@ public:
240 241 virtual void column(const int& column, const GVector& vector);
241 242 virtual void add_to_row(const int& row, const GVector& vector);
242 243 virtual void add_to_column(const int& column, const GVector& vector);
243   - virtual void transpose(void);
244   - virtual void invert(void);
245   - virtual void negate(void);
246   - virtual void abs(void);
247 244 virtual double fill(void) const;
248 245 virtual double min(void) const;
249 246 virtual double max(void) const;
... ... @@ -251,20 +248,24 @@ public:
251 248 virtual std::string print(const GChatter& chatter = NORMAL) const;
252 249  
253 250 // Other methods
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);
258   - void cholesky_decompose(bool compress = true);
259   - GVector cholesky_solver(const GVector& vector, bool compress = true);
260   - void cholesky_invert(bool compress = true);
261   - void set_mem_block(const int& block);
262   - void stack_init(const int& size = 0, const int& entries = 0);
263   - int stack_push_column(const GVector& vector, const int& col);
264   - int stack_push_column(const double* values, const int* rows,
265   - const int& number, const int& col);
266   - void stack_flush(void);
267   - void stack_destroy(void);
  251 + void column(const int& column, const double* values,
  252 + const int* rows, int number);
  253 + void add_to_column(const int& column, const double* values,
  254 + const int* rows, int number);
  255 + GMatrixSparse transpose(void) const;
  256 + GMatrixSparse invert(void) const;
  257 + GVector solve(const GVector& vector) const;
  258 + GMatrixSparse abs(void) const;
  259 + GMatrixSparse cholesky_decompose(bool compress = true) const;
  260 + GVector cholesky_solver(const GVector& vector, bool compress = true) const;
  261 + GMatrixSparse cholesky_invert(bool compress = true) const;
  262 + void set_mem_block(const int& block);
  263 + void stack_init(const int& size = 0, const int& entries = 0);
  264 + int stack_push_column(const GVector& vector, const int& col);
  265 + int stack_push_column(const double* values, const int* rows,
  266 + const int& number, const int& col);
  267 + void stack_flush(void);
  268 + void stack_destroy(void);
268 269  
269 270 private:
270 271 // Private methods
... ... @@ -416,20 +417,6 @@ GMatrixSparse& GMatrixSparse::operator/=(const double& scalar)
416 417  
417 418  
418 419 /***********************************************************************//**
419   - * @brief Negate matrix elements
420   - *
421   - * @return Matrix with negated elements.
422   - ***************************************************************************/
423   -inline
424   -GMatrixSparse GMatrixSparse::operator-(void) const
425   -{
426   - GMatrixSparse result = *this;
427   - result.negate();
428   - return result;
429   -}
430   -
431   -
432   -/***********************************************************************//**
433 420 * @brief Multiply matrix by scalar
434 421 *
435 422 * @param[in] matrix Matrix.
... ... @@ -479,84 +466,4 @@ GMatrixSparse operator/(const GMatrixSparse& matrix, const double& scalar)
479 466 return result;
480 467 }
481 468  
482   -
483   -/***********************************************************************//**
484   - * @brief Return transpose of matrix
485   - *
486   - * @param[in] matrix Matrix.
487   - * @return Transpose of matrix.
488   - ***************************************************************************/
489   -inline
490   -GMatrixSparse transpose(const GMatrixSparse& matrix)
491   -{
492   - GMatrixSparse result = matrix;
493   - result.transpose();
494   - return result;
495   -}
496   -
497   -
498   -/***********************************************************************//**
499   - * @brief Return inverse of matrix
500   - *
501   - * @param[in] matrix Matrix.
502   - * @return Inverse of matrix.
503   - ***************************************************************************/
504   -inline
505   -GMatrixSparse invert(const GMatrixSparse& matrix)
506   -{
507   - GMatrixSparse result = matrix;
508   - result.invert();
509   - return result;
510   -}
511   -
512   -
513   -/***********************************************************************//**
514   - * @brief Return matrix with absolute values of all elements
515   - *
516   - * @param[in] matrix Matrix.
517   - * @return Matrix with elements being the absolute elements of the input
518   - * matrix.
519   - ***************************************************************************/
520   -inline
521   -GMatrixSparse abs(const GMatrixSparse& matrix)
522   -{
523   - GMatrixSparse result = matrix;
524   - result.abs();
525   - return result;
526   -}
527   -
528   -
529   -/***********************************************************************//**
530   - * @brief Return Cholesky decomposition of matrix
531   - *
532   - * @param[in] matrix Matrix.
533   - * @param[in] compress Use matrix compression (defaults to true).
534   - * @return Cholesky decomposition of matrix.
535   - ***************************************************************************/
536   -inline
537   -GMatrixSparse cholesky_decompose(const GMatrixSparse& matrix,
538   - bool compress = true)
539   -{
540   - GMatrixSparse result = matrix;
541   - result.cholesky_decompose(compress);
542   - return result;
543   -}
544   -
545   -
546   -/***********************************************************************//**
547   - * @brief Return inverse matrix using Cholesky decomposition
548   - *
549   - * @param[in] matrix Matrix.
550   - * @param[in] compress Use matrix compression (defaults to true).
551   - * @return Inverse of matrix.
552   - ***************************************************************************/
553   -inline
554   -GMatrixSparse cholesky_invert(const GMatrixSparse& matrix,
555   - bool compress = true)
556   -{
557   - GMatrixSparse result = matrix;
558   - result.cholesky_invert(compress);
559   - return result;
560   -}
561   -
562 469 #endif /* GMATRIXSPARSE_HPP */
... ...
include/GMatrixSymmetric.hpp
  View file @5b16109
... ... @@ -34,6 +34,7 @@
34 34 /* __ Forward declarations _______________________________________________ */
35 35 class GMatrix;
36 36 class GMatrixSparse;
  37 +class GVector;
37 38  
38 39  
39 40 /***********************************************************************//**
... ... @@ -106,10 +107,6 @@ public:
106 107 virtual void column(const int& column, const GVector& vector);
107 108 virtual void add_to_row(const int& row, const GVector& vector);
108 109 virtual void add_to_column(const int& column, const GVector& vector);
109   - virtual void transpose(void);
110   - virtual void invert(void);
111   - virtual void negate(void);
112   - virtual void abs(void);
113 110 virtual double fill(void) const;
114 111 virtual double min(void) const;
115 112 virtual double max(void) const;
... ... @@ -117,11 +114,15 @@ public:
117 114 virtual std::string print(const GChatter& chatter = NORMAL) const;
118 115  
119 116 // Other methods
120   - virtual GMatrix extract_lower_triangle(void) const;
121   - virtual GMatrix extract_upper_triangle(void) const;
122   - virtual void cholesky_decompose(bool compress = true);
123   - virtual GVector cholesky_solver(const GVector& vector, bool compress = true);
124   - virtual void cholesky_invert(bool compress = true);
  117 + GMatrixSymmetric transpose(void) const;
  118 + GMatrixSymmetric invert(void) const;
  119 + GVector solve(const GVector& vector) const;
  120 + GMatrixSymmetric abs(void) const;
  121 + GMatrix extract_lower_triangle(void) const;
  122 + GMatrix extract_upper_triangle(void) const;
  123 + GMatrixSymmetric cholesky_decompose(bool compress = true) const;
  124 + GVector cholesky_solver(const GVector& vector, bool compress = true) const;
  125 + GMatrixSymmetric cholesky_invert(bool compress = true) const;
125 126  
126 127 private:
127 128 // Private methods
... ... @@ -138,16 +139,18 @@ private:
138 139  
139 140  
140 141 /***********************************************************************//**
141   - * @brief Transpose matrix
  142 + * @brief Return transposed matrix
142 143 *
143   - * Transpose matrix. As the transposed matrix of a symmetric matrix is
144   - * identical to the original matrix, this method simply returns without
145   - * taking any action.
  144 + * @return Transposed matrix.
  145 + *
  146 + * Returns transposed matrix of the matrix. As the transposed matrix of a
  147 + * symmetric matrix is identical to the original matrix, this method simply
  148 + * returns the actual matrix.
146 149 ***************************************************************************/
147 150 inline
148   -void GMatrixSymmetric::transpose(void)
  151 +GMatrixSymmetric GMatrixSymmetric::transpose(void) const
149 152 {
150   - return;
  153 + return (*this);
151 154 }
152 155  
153 156  
... ... @@ -250,22 +253,6 @@ GMatrixSymmetric& GMatrixSymmetric::operator/=(const double& scalar)
250 253  
251 254  
252 255 /***********************************************************************//**
253   - * @brief Negate matrix elements
254   - *
255   - * @return Matrix with negated elements.
256   - *
257   - * Returns a matrix where each element is replaced by its negative element.
258   - ***************************************************************************/
259   -inline
260   -GMatrixSymmetric GMatrixSymmetric::operator-(void) const
261   -{
262   - GMatrixSymmetric result = *this;
263   - result.negate();
264   - return result;
265   -}
266   -
267   -
268   -/***********************************************************************//**
269 256 * @brief Multiply matrix by scalar
270 257 *
271 258 * @param[in] matrix Matrix.
... ... @@ -318,93 +305,4 @@ GMatrixSymmetric operator/(const GMatrixSymmetric& matrix, const double& scalar)
318 305 return result;
319 306 }
320 307  
321   -
322   -/***********************************************************************//**
323   - * @brief Return transpose of matrix
324   - *
325   - * @param[in] matrix Matrix.
326   - * @return Transpose of matrix.
327   - *
328   - * Returns the transpose of the matrix. This is trivial for a symmetric
329   - * matrix, as its transpose is the original matrix.
330   - ***************************************************************************/
331   -inline
332   -GMatrixSymmetric transpose(const GMatrixSymmetric& matrix)
333   -{
334   - return matrix;
335   -}
336   -
337   -
338   -/***********************************************************************//**
339   - * @brief Return inverse of matrix
340   - *
341   - * @param[in] matrix Matrix.
342   - * @return Inverse of matrix.
343   - *
344   - * Returns the inverse of the matrix.
345   - ***************************************************************************/
346   -inline
347   -GMatrixSymmetric invert(const GMatrixSymmetric& matrix)
348   -{
349   - GMatrixSymmetric result = matrix;
350   - result.invert();
351   - return result;
352   -}
353   -
354   -
355   -/***********************************************************************//**
356   - * @brief Return matrix with absolute values of all elements
357   - *
358   - * @param[in] matrix Matrix.
359   - * @return Matrix with elements being the absolute elements of the input
360   - * matrix.
361   - *
362   - * Returns a matrix where each element is replaced by its absolute value.
363   - ***************************************************************************/
364   -inline
365   -GMatrixSymmetric abs(const GMatrixSymmetric& matrix)
366   -{
367   - GMatrixSymmetric result = matrix;
368   - result.abs();
369   - return result;
370   -}
371   -
372   -
373   -/***********************************************************************//**
374   - * @brief Return Cholesky decomposition of matrix
375   - *
376   - * @param[in] matrix Matrix.
377   - * @param[in] compress Use matrix compression (defaults to true).
378   - * @return Cholesky decomposition of matrix.
379   - *
380   - * Returns a Cholesky decomposition of a matrix.
381   - ***************************************************************************/
382   -inline
383   -GMatrixSymmetric cholesky_decompose(const GMatrixSymmetric& matrix,
384   - bool compress = true)
385   -{
386   - GMatrixSymmetric result = matrix;
387   - result.cholesky_decompose(compress);
388   - return result;
389   -}
390   -
391   -
392   -/***********************************************************************//**
393   - * @brief Return inverse matrix using Cholesky decomposition
394   - *
395   - * @param[in] matrix Matrix.
396   - * @param[in] compress Use matrix compression (defaults to true).
397   - * @return Inverse of matrix.
398   - *
399   - * Returns the inverse of the matrix by performing a Cholesky decomposition.
400   - ***************************************************************************/
401   -inline
402   -GMatrixSymmetric cholesky_invert(const GMatrixSymmetric& matrix,
403   - bool compress = true)
404   -{
405   - GMatrixSymmetric result = matrix;
406   - result.cholesky_invert(compress);
407   - return result;
408   -}
409   -
410 308 #endif /* GMATRIXSYMMETRIC_HPP */
... ...
pyext/GMatrix.i
  View file @5b16109
... ... @@ -72,21 +72,20 @@ public:
72 72 virtual void column(const int& column, const GVector& vector);
73 73 virtual void add_to_row(const int& row, const GVector& vector);
74 74 virtual void add_to_column(const int& column, const GVector& vector);
75   - virtual void transpose(void);
76   - virtual void invert(void);
77   - virtual void negate(void);
78   - virtual void abs(void);
79 75 virtual double fill(void) const;
80 76 virtual double min(void) const;
81 77 virtual double max(void) const;
82 78 virtual double sum(void) const;
83 79  
84 80 // Other methods
85   - virtual GMatrix extract_lower_triangle(void) const;
86   - virtual GMatrix extract_upper_triangle(void) const;
87   - virtual void eulerx(const double& angle);
88   - virtual void eulery(const double& angle);
89   - virtual void eulerz(const double& angle);
  81 + GMatrix transpose(void) const;
  82 + GMatrix invert(void) const;
  83 + GMatrix abs(void) const;
  84 + GMatrix extract_lower_triangle(void) const;
  85 + GMatrix extract_upper_triangle(void) const;
  86 + void eulerx(const double& angle);
  87 + void eulery(const double& angle);
  88 + void eulerz(const double& angle);
90 89 };
91 90  
92 91  
... ... @@ -120,11 +119,3 @@ public:
120 119 return (*self);
121 120 }
122 121 };
123   -
124   -
125   -/***********************************************************************//**
126   - * @brief GMatrix friends
127   - ***************************************************************************/
128   -GMatrix transpose(const GMatrix& matrix);
129   -GMatrix invert(const GMatrix& matrix);
130   -GMatrix abs(const GMatrix& matrix);
... ...
pyext/GMatrixBase.i
  View file @5b16109
... ... @@ -55,10 +55,6 @@ public:
55 55 virtual void column(const int& column, const GVector& vector) = 0;
56 56 virtual void add_to_row(const int& row, const GVector& vector) = 0;
57 57 virtual void add_to_column(const int& column, const GVector& vector) = 0;
58   - virtual void transpose(void) = 0;
59   - virtual void invert(void) = 0;
60   - virtual void negate(void) = 0;
61   - virtual void abs(void) = 0;
62 58 virtual double fill(void) const = 0;
63 59 virtual double min(void) const = 0;
64 60 virtual double max(void) const = 0;
... ...
pyext/GMatrixSparse.i
  View file @5b16109
... ... @@ -76,30 +76,29 @@ public:
76 76 virtual void column(const int& column, const GVector& vector);
77 77 virtual void add_to_row(const int& row, const GVector& vector);
78 78 virtual void add_to_column(const int& column, const GVector& vector);
79   - virtual void transpose(void);
80   - virtual void invert(void);
81   - virtual void negate(void);
82   - virtual void abs(void);
83 79 virtual double fill(void) const;
84 80 virtual double min(void) const;
85 81 virtual double max(void) const;
86 82 virtual double sum(void) const;
87 83  
88 84 // Other methods
89   - void column(const int& column, const double* values,
90   - const int* rows, int number);
91   - void add_to_column(const int& column, const double* values,
92   - const int* rows, int number);
93   - void cholesky_decompose(bool compress = true);
94   - GVector cholesky_solver(const GVector& vector, bool compress = true);
95   - void cholesky_invert(bool compress = true);
96   - void set_mem_block(const int& block);
97   - void stack_init(const int& size = 0, const int& entries = 0);
98   - int stack_push_column(const GVector& vector, const int& col);
99   - int stack_push_column(const double* values, const int* rows,
100   - const int& number, const int& col);
101   - void stack_flush(void);
102   - void stack_destroy(void);
  85 + void column(const int& column, const double* values,
  86 + const int* rows, int number);
  87 + void add_to_column(const int& column, const double* values,
  88 + const int* rows, int number);
  89 + GMatrixSparse transpose(void) const;
  90 + GMatrixSparse invert(void) const;
  91 + GMatrixSparse abs(void) const;
  92 + GMatrixSparse cholesky_decompose(bool compress = true);
  93 + GVector cholesky_solver(const GVector& vector, bool compress = true);
  94 + GMatrixSparse cholesky_invert(bool compress = true);
  95 + void set_mem_block(const int& block);
  96 + void stack_init(const int& size = 0, const int& entries = 0);
  97 + int stack_push_column(const GVector& vector, const int& col);
  98 + int stack_push_column(const double* values, const int* rows,
  99 + const int& number, const int& col);
  100 + void stack_flush(void);
  101 + void stack_destroy(void);
103 102 };
104 103  
105 104  
... ... @@ -133,11 +132,3 @@ public:
133 132 return (*self);
134 133 }
135 134 };
136   -
137   -
138   -/***********************************************************************//**
139   - * @brief GMatrixSparse friends
140   - ***************************************************************************/
141   -GMatrixSparse transpose(const GMatrixSparse& matrix);
142   -GMatrixSparse invert(const GMatrixSparse& matrix);
143   -GMatrixSparse abs(const GMatrixSparse& matrix);
... ...
pyext/GMatrixSymmetric.i
  View file @5b16109
... ... @@ -70,21 +70,21 @@ public:
70 70 virtual void column(const int& column, const GVector& vector);
71 71 virtual void add_to_row(const int& row, const GVector& vector);
72 72 virtual void add_to_column(const int& column, const GVector& vector);
73   - virtual void transpose(void);
74   - virtual void invert(void);
75   - virtual void negate(void);
76   - virtual void abs(void);
77 73 virtual double fill(void) const;
78 74 virtual double min(void) const;
79 75 virtual double max(void) const;
80 76 virtual double sum(void) const;
81 77  
82 78 // Other methods
83   - virtual GMatrix extract_lower_triangle(void) const;
84   - virtual GMatrix extract_upper_triangle(void) const;
85   - virtual void cholesky_decompose(bool compress = true);
86   - virtual GVector cholesky_solver(const GVector& vector, bool compress = true);
87   - virtual void cholesky_invert(bool compress = true);
  79 + GMatrixSymmetric transpose(void) const;
  80 + GMatrixSymmetric invert(void) const;
  81 + GVector solve(const GVector& vector) const;
  82 + GMatrixSymmetric abs(void) const;
  83 + GMatrix extract_lower_triangle(void) const;
  84 + GMatrix extract_upper_triangle(void) const;
  85 + GMatrixSymmetric cholesky_decompose(bool compress = true) const;
  86 + GVector cholesky_solver(const GVector& vector, bool compress = true) const;
  87 + GMatrixSymmetric cholesky_invert(bool compress = true) const;
88 88 };
89 89  
90 90  
... ... @@ -118,11 +118,3 @@ public:
118 118 return (*self);
119 119 }
120 120 };
121   -
122   -
123   -/***********************************************************************//**
124   - * @brief GMatrixSymmetric friends
125   - ***************************************************************************/
126   -GMatrixSymmetric transpose(const GMatrixSymmetric& matrix);
127   -GMatrixSymmetric invert(const GMatrixSymmetric& matrix);
128   -GMatrixSymmetric abs(const GMatrixSymmetric& matrix);
... ...
src/linalg/GMatrix.cpp
  View file @5b16109
... ... @@ -51,6 +51,7 @@
51 51 #define G_ADD_TO_ROW "GMatrix::add_to_row(int&, GVector&)"
52 52 #define G_ADD_TO_COLUMN "GMatrix::add_to_column(int&, GVector&)"
53 53 #define G_INVERT "GMatrix::invert()"
  54 +#define G_SOLVE "GMatrix::solve(GVector&)"
54 55 #define G_EXTRACT_LOWER "GMatrix::extract_lower_triangle()"
55 56 #define G_EXTRACT_UPPER "GMatrix::extract_upper_triangle()"
56 57  
... ... @@ -305,6 +306,30 @@ GVector GMatrix::operator*(const GVector& vector) const
305 306  
306 307  
307 308 /***********************************************************************//**
  309 + * @brief Negate matrix elements
  310 + *
  311 + * @return Matrix with negated elements.
  312 + *
  313 + * Returns a matrix where each element has been replaced by its negative
  314 + * element.
  315 + ***************************************************************************/
  316 +GMatrix GMatrix::operator-(void) const
  317 +{
  318 + // Copy matrix
  319 + GMatrix matrix = *this;
  320 +
  321 + // Take the absolute value of all matrix elements
  322 + double* ptr = matrix.m_data;
  323 + for (int i = 0; i < matrix.m_elements; ++i, ++ptr) {
  324 + *ptr = -(*ptr);
  325 + }
  326 +
  327 + // Return matrix
  328 + return matrix;
  329 +}
  330 +
  331 +
  332 +/***********************************************************************//**
308 333 * @brief Unary matrix addition operator
309 334 *
310 335 * @param[in] matrix Matrix.
... ... @@ -773,98 +798,106 @@ void GMatrix::add_to_column(const int&amp; column, const GVector&amp; vector)
773 798  
774 799  
775 800 /***********************************************************************//**
776   - * @brief Transpose matrix
  801 + * @brief Return transposed matrix
  802 + *
  803 + * @return Transposed matrix.
777 804 *
778   - * The transpose operation exchanges the number of rows against the number
779   - * of columns. For a square matrix the exchange is done inplace. Otherwise
780   - * a copy of the matrix is made.
  805 + * Returns transposed matrix of the matrix.
781 806 ***************************************************************************/
782   -void GMatrix::transpose(void)
  807 +GMatrix GMatrix::transpose(void) const
783 808 {
784   - // Case A: Matrix is square then simply swap the elements
785   - if (m_rows == m_cols) {
786   - double swap;
787   - double* ptr_dst;
788   - double* ptr_src;
789   - for (int row = 0; row < m_rows; ++row) {
790   - for (int col = row; col < m_cols; ++col) {
791   - ptr_dst = m_data + m_cols*row + col;
792   - ptr_src = m_data + m_rows*col + row;
793   - swap = *ptr_dst;
794   - *ptr_dst = *ptr_src;
795   - *ptr_src = swap;
796   - }
797   - }
798   - }
  809 + // Allocate result matrix
  810 + GMatrix matrix(m_cols, m_rows);
799 811  
800   - // Case B: Non-rectangular transpose. No code optimization has been
801   - // done so far.
802   - else {
803   - GMatrix result(m_cols, m_rows);
804   - for (int row = 0; row < m_rows; ++row) {
805   - for (int col = 0; col < m_cols; ++col) {
806   - result(col, row) = (*this)(row, col);
807   - }
  812 + // Transpose matrix
  813 + for (int row = 0; row < m_rows; ++row) {
  814 + for (int col = 0; col < m_cols; ++col) {
  815 + matrix(col, row) = (*this)(row, col);
808 816 }
809   - *this = result;
810 817 }
811 818  
812   - // Return
813   - return;
  819 + // Return matrix
  820 + return matrix;
814 821 }
815 822  
816 823  
817 824 /***********************************************************************//**
818   - * @brief Invert matrix
  825 + * @brief Return inverted matrix
  826 + *
  827 + * @return Inverted matrix.
819 828 *
820 829 * @exception GException::feature_not_implemented
821 830 * Feature not yet implemented.
822 831 *
  832 + * Returns inverse of matrix.
  833 + *
823 834 * @todo Needs to be implemented.
824 835 ***************************************************************************/
825   -void GMatrix::invert(void)
  836 +GMatrix GMatrix::invert(void) const
826 837 {
  838 + // Allocate result matrix
  839 + GMatrix matrix(m_cols, m_rows);
  840 +
827 841 // Throw exception
828 842 throw GException::feature_not_implemented(G_INVERT);
829 843  
830   - // Return
831   - return;
  844 + // Return matrix
  845 + return matrix;
832 846 }
833 847  
834 848  
835 849 /***********************************************************************//**
836   - * @brief Negate all matrix elements
  850 + * @brief Solves linear matrix equation
  851 + *
  852 + * @param[in] vector Solution vector.
  853 + *
  854 + * @exception GException::feature_not_implemented
  855 + * Feature not yet implemented.
  856 + *
  857 + * Solves the linear equation
  858 + *
  859 + * \f[M \times {\tt solution} = {\tt vector} \f]
  860 + *
  861 + * where \f$M\f$ is the matrix, \f${\tt vector}\f$ is the result, and
  862 + * \f${\tt solution}\f$ is the solution.
837 863 *
838   - * Negates all matrix elements.
  864 + * @todo Needs to be implemented.
839 865 ***************************************************************************/
840   -void GMatrix::negate(void)
  866 +GVector GMatrix::solve(const GVector& vector) const
841 867 {
842   - // Negate all matrix elements
843   - double* ptr = m_data;
844   - for (int i = 0; i < m_elements; ++i, ++ptr) {
845   - *ptr = -(*ptr);
846   - }
  868 + // Allocate result vector
  869 + GVector result;
847 870  
848   - // Return
849   - return;
  871 + // Throw exception
  872 + throw GException::feature_not_implemented(G_SOLVE);
  873 +
  874 + // Return vector
  875 + return result;
850 876 }
851 877  
852 878  
853 879 /***********************************************************************//**
854   - * @brief Take absolute value of matrix elements
  880 + * @brief Return absolute of matrix
855 881 *
856   - * Replaces all elements of the matrix by their absolute values.
  882 + * @return Absolute of matrix
  883 + *
  884 + * Returns matrix where all elements of the matrix have been replaced by
  885 + * their absolute values.
857 886 ***************************************************************************/
858   -void GMatrix::abs(void)
  887 +GMatrix GMatrix::abs(void) const
859 888 {
  889 + // Allocate result matrix
  890 + GMatrix matrix(m_rows, m_cols);
  891 +
860 892 // Take the absolute value of all matrix elements
861   - double* ptr = m_data;
862   - for (int i = 0; i < m_elements; ++i, ++ptr) {
863   - *ptr = std::abs(*ptr);
  893 + double* src = m_data;
  894 + double* dst = matrix.m_data;
  895 + for (int i = 0; i < m_elements; ++i) {
  896 + *dst++ = std::abs(*src++);
864 897 }
865 898  
866   - // Return
867   - return;
  899 + // Return matrix
  900 + return matrix;
868 901 }
869 902  
870 903  
... ...
src/linalg/GMatrixSparse.cpp
  View file @5b16109
... ... @@ -453,6 +453,30 @@ GVector GMatrixSparse::operator*(const GVector&amp; vector) const
453 453  
454 454  
455 455 /***********************************************************************//**
  456 + * @brief Negate matrix elements
  457 + *
  458 + * @return Matrix with negated elements.
  459 + ***************************************************************************/
  460 +GMatrixSparse GMatrixSparse::operator-(void) const
  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 +/***********************************************************************//**
456 480 * @brief Equalty operator
457 481 *
458 482 * @param[in] matrix Matrix.
... ... @@ -1380,81 +1404,103 @@ void GMatrixSparse::add_to_column(const int&amp; column, const double* values,
1380 1404  
1381 1405  
1382 1406 /***********************************************************************//**
1383   - * @brief Transpose matrix
  1407 + * @brief Return transposed matrix
1384 1408 *
1385   - * The transpose operation exchanges the number of rows against the number
1386   - * of columns.
  1409 + * @return Transposed matrix.
  1410 + *
  1411 + * Returns transposed matrix of the matrix.
1387 1412 ***************************************************************************/
1388   -void GMatrixSparse::transpose(void)
  1413 +GMatrixSparse GMatrixSparse::transpose(void) const
1389 1414 {
  1415 + // Copy matrix
  1416 + GMatrixSparse matrix = *this;
  1417 +
1390 1418 // Fill pending element
1391   - fill_pending();
  1419 + matrix.fill_pending();
1392 1420  
1393 1421 // Compute the transpose
1394   - *this = cs_transpose(*this, 1);
  1422 + matrix = cs_transpose(matrix, 1);
1395 1423  
1396   - // Return
1397   - return;
  1424 + // Return matrix
  1425 + return matrix;
1398 1426 }
1399 1427  
1400 1428  
1401 1429 /***********************************************************************//**
1402   - * @brief Invert matrix
  1430 + * @brief Return inverted matrix
  1431 + *
  1432 + * @return Inverted matrix.
1403 1433 *
1404   - * Inverts the matrix using the Cholesky decomposition. This does not work
1405   - * on any kind of matrix.
  1434 + * Returns inverse of matrix. Inversion is done for the moment using Cholesky
  1435 + * decomposition. This does not work on any kind of matrix.
1406 1436 *
1407 1437 * @todo Specify in documentation for which kind of matrix the method works.
1408 1438 ***************************************************************************/
1409   -void GMatrixSparse::invert(void)
  1439 +GMatrixSparse GMatrixSparse::invert(void) const
1410 1440 {
  1441 + // Allocate result matrix
  1442 + GMatrixSparse matrix(m_cols, m_rows);
  1443 +
1411 1444 // Invert matrix
1412   - cholesky_invert(true);
  1445 + matrix.cholesky_invert(true);
1413 1446  
1414   - // Return
1415   - return;
  1447 + // Return matrix
  1448 + return matrix;
1416 1449 }
1417 1450  
1418 1451  
1419 1452 /***********************************************************************//**
1420   - * @brief Negate all matrix elements
  1453 + * @brief Solves linear matrix equation
1421 1454 *
1422   - * Negates all matrix elements.
  1455 + * @param[in] vector Solution vector.
  1456 + *
  1457 + * Solves the linear equation
  1458 + *
  1459 + * \f[M \times {\tt solution} = {\tt vector} \f]
  1460 + *
  1461 + * where \f$M\f$ is the matrix, \f${\tt vector}\f$ is the result, and
  1462 + * \f${\tt solution}\f$ is the solution. Solving is done using Cholesky
  1463 + * decomposition. This does not work on any kind of matrix.
  1464 + *
  1465 + * @todo Specify in documentation for which kind of matrix the method works.
1423 1466 ***************************************************************************/
1424   -void GMatrixSparse::negate(void)
  1467 +GVector GMatrixSparse::solve(const GVector& vector) const
1425 1468 {
1426   - // Fill pending element
1427   - fill_pending();
  1469 + // Get Cholesky decomposition of matrix
  1470 + GMatrixSparse decomposition = cholesky_decompose(true);
1428 1471  
1429   - // Negate all matrix elements
1430   - double* ptr = m_data;
1431   - for (int i = 0; i < m_elements; ++i, ++ptr) {
1432   - *ptr = -(*ptr);
1433   - }
1434   -
1435   - // Return
1436   - return;
  1472 + // Solve linear equation
  1473 + GVector result = decomposition.cholesky_solver(vector);
  1474 +
  1475 + // Return result
  1476 + return result;
1437 1477 }
1438 1478  
1439 1479  
1440 1480 /***********************************************************************//**
1441   - * @brief Take absolute value of matrix elements
  1481 + * @brief Return absolute of matrix
  1482 + *
  1483 + * @return Absolute of matrix
1442 1484 *
1443   - * Replaces all elements of the matrix by their absolute values.
  1485 + * Returns matrix where all elements of the matrix have been replaced by
  1486 + * their absolute values.
1444 1487 ***************************************************************************/
1445   -void GMatrixSparse::abs(void)
  1488 +GMatrixSparse GMatrixSparse::abs(void) const
1446 1489 {
  1490 + // Copy matrix
  1491 + GMatrixSparse matrix = *this;
  1492 +
1447 1493 // Fill pending element
1448   - fill_pending();
  1494 + matrix.fill_pending();
1449 1495  
1450 1496 // Take the absolute value of all matrix elements
1451   - double* ptr = m_data;
1452   - for (int i = 0; i < m_elements; ++i, ++ptr) {
  1497 + double* ptr = matrix.m_data;
  1498 + for (int i = 0; i < matrix.m_elements; ++i, ++ptr) {
1453 1499 *ptr = std::abs(*ptr);
1454 1500 }
1455   -
1456   - // Return
1457   - return;
  1501 +
  1502 + // Return matrix
  1503 + return matrix;
1458 1504 }
1459 1505  
1460 1506  
... ... @@ -1561,25 +1607,29 @@ double GMatrixSparse::sum(void) const
1561 1607  
1562 1608  
1563 1609 /***********************************************************************//**
1564   - * @brief Perform a Cholesky decomposition
  1610 + * @brief Return Cholesky decomposition
1565 1611 *
1566   - * @param[in] compress Use zero-row/column compression (default: true).
  1612 + * @param[in] compress Use zero-row/column compression (defaults to true).
  1613 + * @return Cholesky decomposition of matrix
1567 1614 *
1568   - * Cholesky decomposition of a sparse matrix. The decomposition is stored
1569   - * within the GMatrixSparse without destroying the original matrix.
  1615 + * Returns the Cholesky decomposition of a sparse matrix. The decomposition
  1616 + * is stored within a GMatrixSparse object.
1570 1617 ***************************************************************************/
1571   -void GMatrixSparse::cholesky_decompose(bool compress)
  1618 +GMatrixSparse GMatrixSparse::cholesky_decompose(bool compress) const
1572 1619 {
  1620 + // Create copy of matrix
  1621 + GMatrixSparse matrix = *this;
  1622 +
1573 1623 // Save original matrix size
1574   - int matrix_rows = m_rows;
1575   - int matrix_cols = m_cols;
  1624 + int matrix_rows = matrix.m_rows;
  1625 + int matrix_cols = matrix.m_cols;
1576 1626  
1577 1627 // Delete any existing symbolic and numeric analysis object and reset
1578 1628 // pointers
1579   - if (m_symbolic != NULL) delete m_symbolic;
1580   - if (m_numeric != NULL) delete m_numeric;
1581   - m_symbolic = NULL;
1582   - m_numeric = NULL;
  1629 + if (matrix.m_symbolic != NULL) delete matrix.m_symbolic;
  1630 + if (matrix.m_numeric != NULL) delete matrix.m_numeric;
  1631 + matrix.m_symbolic = NULL;
  1632 + matrix.m_numeric = NULL;
1583 1633  
1584 1634 // Allocate symbolic analysis object
1585 1635 GSparseSymbolic* symbolic = new GSparseSymbolic();
... ... @@ -1590,42 +1640,42 @@ void GMatrixSparse::cholesky_decompose(bool compress)
1590 1640 GSparseNumeric numeric;
1591 1641  
1592 1642 // Fill pending element into matrix
1593   - fill_pending();
  1643 + matrix.fill_pending();
1594 1644  
1595 1645 // Remove rows and columns containing only zeros if matrix compression
1596 1646 // has been selected
1597 1647 if (compress) {
1598   - remove_zero_row_col();
  1648 + matrix.remove_zero_row_col();
1599 1649 }
1600 1650  
1601 1651 // Ordering an symbolic analysis of matrix. This sets up an array 'pinv'
1602 1652 // which contains the fill-in reducing permutations
1603   - symbolic->cholesky_symbolic_analysis(1, *this);
  1653 + symbolic->cholesky_symbolic_analysis(1, matrix);
1604 1654  
1605 1655 // Store symbolic pointer in sparse matrix object
1606   - m_symbolic = symbolic;
  1656 + matrix.m_symbolic = symbolic;
1607 1657  
1608 1658 // Perform numeric Cholesky decomposition
1609   - numeric.cholesky_numeric_analysis(*this, *symbolic);
  1659 + numeric.cholesky_numeric_analysis(matrix, *symbolic);
1610 1660  
1611 1661 // Copy L matrix into this object
1612   - free_elements(0, m_elements);
1613   - alloc_elements(0, numeric.m_L->m_elements);
1614   - for (int i = 0; i < m_elements; ++i) {
1615   - m_data[i] = numeric.m_L->m_data[i];
1616   - m_rowinx[i] = numeric.m_L->m_rowinx[i];
  1662 + matrix.free_elements(0, matrix.m_elements);
  1663 + matrix.alloc_elements(0, numeric.m_L->m_elements);
  1664 + for (int i = 0; i < matrix.m_elements; ++i) {
  1665 + matrix.m_data[i] = numeric.m_L->m_data[i];
  1666 + matrix.m_rowinx[i] = numeric.m_L->m_rowinx[i];
1617 1667 }
1618   - for (int col = 0; col <= m_cols; ++col) {
1619   - m_colstart[col] = numeric.m_L->m_colstart[col];
  1668 + for (int col = 0; col <= matrix.m_cols; ++col) {
  1669 + matrix.m_colstart[col] = numeric.m_L->m_colstart[col];
1620 1670 }
1621 1671  
1622 1672 // Insert zero rows and columns if they have been removed previously.
1623 1673 if (compress) {
1624   - insert_zero_row_col(matrix_rows, matrix_cols);
  1674 + matrix.insert_zero_row_col(matrix_rows, matrix_cols);
1625 1675 }
1626 1676  
1627   - // Return
1628   - return;
  1677 + // Return matrix
  1678 + return matrix;
1629 1679 }
1630 1680  
1631 1681  
... ... @@ -1644,7 +1694,8 @@ void GMatrixSparse::cholesky_decompose(bool compress)
1644 1694 * This function is to be applied on a GMatrixSparse matrix for which a
1645 1695 * Choleksy factorization has been produced using 'cholesky_decompose'.
1646 1696 ***************************************************************************/
1647   -GVector GMatrixSparse::cholesky_solver(const GVector& vector, bool compress)
  1697 +GVector GMatrixSparse::cholesky_solver(const GVector& vector,
  1698 + bool compress) const
1648 1699 {
1649 1700 // Dump header
1650 1701 #if defined(G_DEBUG_SPARSE_COMPRESSION)
... ... @@ -1862,22 +1913,20 @@ GVector GMatrixSparse::cholesky_solver(const GVector&amp; vector, bool compress)
1862 1913  
1863 1914  
1864 1915 /***********************************************************************//**
1865   - * @brief Cholesky invert
  1916 + * @brief Invert matrix using a Cholesky decomposition
1866 1917 *
1867   - * @param[in] compress Use zero-row/column compression (default: true).
  1918 + * @param[in] compress Use zero-row/column compression (defaults to true).
  1919 + * @return Inverted matrix.
1868 1920 *
1869   - * @exception GException::matrix_zero
1870   - * All matrix elements are zero.
1871   - *
1872   - * Inplace matrix inversion using the Cholesky decomposition.
  1921 + * Inverts the matrix using a Cholesky decomposition.
1873 1922 ***************************************************************************/
1874   -void GMatrixSparse::cholesky_invert(bool compress)
  1923 +GMatrixSparse GMatrixSparse::cholesky_invert(bool compress) const
1875 1924 {
1876 1925 // Generate Cholesky decomposition of matrix
1877   - cholesky_decompose(compress);
  1926 + GMatrixSparse decomposition = cholesky_decompose(compress);
1878 1927  
1879 1928 // Allocate result matrix and unit vector
1880   - GMatrixSparse result(m_rows, m_cols);
  1929 + GMatrixSparse matrix(m_rows, m_cols);
1881 1930 GVector unit(m_rows);
1882 1931  
1883 1932 // Column-wise solving of the problem
... ... @@ -1887,27 +1936,25 @@ void GMatrixSparse::cholesky_invert(bool compress)
1887 1936 unit[col] = 1.0;
1888 1937  
1889 1938 // Solve for column
1890   - GVector x = cholesky_solver(unit, compress);
  1939 + GVector x = decomposition.cholesky_solver(unit, compress);
1891 1940  
1892 1941 // Insert column in matrix
1893   - result.column(col, x);
  1942 + matrix.column(col, x);
1894 1943  
1895 1944 // Clear unit vector for next round
1896 1945 unit[col] = 0.0;
1897 1946  
1898 1947 }
1899 1948  
1900   - // Assign result
1901   - *this = result;
1902   -
1903   - // Return
1904   - return;
  1949 + // Return matrix
  1950 + return matrix;
1905 1951 }
1906 1952  
1907 1953  
1908 1954 /***********************************************************************//**
1909 1955 * @brief Print matrix
1910 1956 *
  1957 + * @param[in] chatter Chattiness (defaults to NORMAL).
1911 1958 * @return String containing matrix information
1912 1959 ***************************************************************************/
1913 1960 std::string GMatrixSparse::print(const GChatter& chatter) const
... ...
src/linalg/GMatrixSymmetric.cpp
  View file @5b16109
... ... @@ -356,6 +356,30 @@ GVector GMatrixSymmetric::operator*(const GVector&amp; vector) const
356 356  
357 357  
358 358 /***********************************************************************//**
  359 + * @brief Negate matrix elements
  360 + *
  361 + * @return Matrix with negated elements.
  362 + *
  363 + * Returns a matrix where each element has been replaced by its negative
  364 + * element.
  365 + ***************************************************************************/
  366 +GMatrixSymmetric GMatrixSymmetric::operator-(void) const
  367 +{
  368 + // Copy matrix
  369 + GMatrixSymmetric matrix = *this;
  370 +
  371 + // Take the absolute value of all matrix elements
  372 + double* ptr = matrix.m_data;
  373 + for (int i = 0; i < matrix.m_elements; ++i, ++ptr) {
  374 + *ptr = -(*ptr);
  375 + }
  376 +
  377 + // Return matrix
  378 + return matrix;
  379 +}
  380 +
  381 +
  382 +/***********************************************************************//**
359 383 * @brief Unary matrix addition operator
360 384 *
361 385 * @param[in] matrix Matrix.
... ... @@ -687,56 +711,78 @@ void GMatrixSymmetric::add_to_column(const int&amp; column, const GVector&amp; vector)
687 711  
688 712  
689 713 /***********************************************************************//**
690   - * @brief Invert matrix
  714 + * @brief Return inverted matrix
691 715 *
692   - * Inverts the matrix using the Cholesky decomposition. This does not work
693   - * on any kind of matrix.
  716 + * @return Inverted matrix.
  717 + *
  718 + * Returns inverse of matrix. Inversion is done for the moment using Cholesky
  719 + * decomposition. This does not work on any kind of matrix.
694 720 *
695 721 * @todo Specify in documentation for which kind of matrix the method works.
696 722 ***************************************************************************/
697   -void GMatrixSymmetric::invert(void)
  723 +GMatrixSymmetric GMatrixSymmetric::invert(void) const
698 724 {
  725 + // Allocate result matrix
  726 + GMatrixSymmetric matrix(m_cols, m_rows);
  727 +
699 728 // Invert matrix
700   - cholesky_invert(true);
  729 + matrix.cholesky_invert(true);
701 730  
702   - // Return
703   - return;
  731 + // Return matrix
  732 + return matrix;
704 733 }
705 734  
706 735  
707 736 /***********************************************************************//**
708   - * @brief Negate all matrix elements
  737 + * @brief Solves linear matrix equation
  738 + *
  739 + * @param[in] vector Solution vector.
  740 + *
  741 + * Solves the linear equation
709 742 *
710   - * Negates all matrix elements.
  743 + * \f[M \times {\tt solution} = {\tt vector} \f]
  744 + *
  745 + * where \f$M\f$ is the matrix, \f${\tt vector}\f$ is the result, and
  746 + * \f${\tt solution}\f$ is the solution. Solving is done using Cholesky
  747 + * decomposition. This does not work on any kind of matrix.
  748 + *
  749 + * @todo Specify in documentation for which kind of matrix the method works.
711 750 ***************************************************************************/
712   -void GMatrixSymmetric::negate(void)
  751 +GVector GMatrixSymmetric::solve(const GVector& vector) const
713 752 {
714   - // Negate all matrix elements
715   - double* ptr = m_data;
716   - for (int i = 0; i < m_elements; ++i, ++ptr) {
717   - *ptr = -(*ptr);
718   - }
719   -
720   - // Return
721   - return;
  753 + // Get Cholesky decomposition of matrix
  754 + GMatrixSymmetric decomposition = cholesky_decompose(true);
  755 +
  756 + // Solve linear equation
  757 + GVector result = decomposition.cholesky_solver(vector);
  758 +
  759 + // Return result
  760 + return result;
722 761 }
723 762  
724 763  
725 764 /***********************************************************************//**
726   - * @brief Take absolute value of matrix elements
  765 + * @brief Return absolute of matrix
  766 + *
  767 + * @return Absolute of matrix
727 768 *
728   - * Replaces all elements of the matrix by their absolute values.
  769 + * Returns matrix where all elements of the matrix have been replaced by
  770 + * their absolute values.
729 771 ***************************************************************************/
730   -void GMatrixSymmetric::abs(void)
  772 +GMatrixSymmetric GMatrixSymmetric::abs(void) const
731 773 {
  774 + // Allocate result matrix
  775 + GMatrixSymmetric matrix(m_rows, m_cols);
  776 +
732 777 // Take the absolute value of all matrix elements
733   - double* ptr = m_data;
734   - for (int i = 0; i < m_elements; ++i, ++ptr) {
735   - *ptr = std::abs(*ptr);
  778 + double* src = m_data;
  779 + double* dst = matrix.m_data;
  780 + for (int i = 0; i < m_elements; ++i) {
  781 + *dst++ = std::abs(*src++);
736 782 }
737   -
738   - // Return
739   - return;
  783 +
  784 + // Return matrix
  785 + return matrix;
740 786 }
741 787  
742 788  
... ... @@ -848,60 +894,69 @@ GMatrix GMatrixSymmetric::extract_upper_triangle(void) const
848 894  
849 895  
850 896 /***********************************************************************//**
851   - * @brief Perform a Cholesky decomposition
  897 + * @brief Return Cholesky decomposition
852 898 *
853   - * @param[in] compress Use zero-row/column compression (default: true).
  899 + * @param[in] compress Use zero-row/column compression (defaults to true).
  900 + * @return Cholesky decomposition of matrix
854 901 *
855 902 * @exception GException::matrix_not_pos_definite
856 903 * Matrix is not positive definite.
857 904 * @exception GException::matrix_zero
858 905 * All matrix elements are zero.
859 906 *
860   - * Inplace Cholesky decomposition inspired by Numerical Recipes algorithm.
  907 + * Returns the Cholesky decomposition of a sparse matrix. The decomposition
  908 + * is stored within a GMatrixSymmetric object.
  909 + *
  910 + * The method is inspired by the algorithm found in Numerical Recipes.
861 911 * The decomposition, which is a matrix occupying only the lower triange,
862 912 * is stored in the elements of the symmetric matrix. To visualise the
863 913 * matrix one has to use 'lower_triangle()' to extract the relevant part.
864 914 * Case A operates on a full matrix, Case B operates on a (logically)
865 915 * compressed matrix where zero rows/columns have been removed.
866 916 ***************************************************************************/
867   -void GMatrixSymmetric::cholesky_decompose(bool compress)
  917 +GMatrixSymmetric GMatrixSymmetric::cholesky_decompose(bool compress) const
868 918 {
  919 + // Create copy of matrix
  920 + GMatrixSymmetric matrix = *this;
  921 +
869 922 // Set-up incides of non zero rows if matrix compression is requested
870 923 if (compress) {
871   - set_inx();
  924 + matrix.set_inx();
872 925 }
873 926  
874 927 // Check if zero-row/col compression is needed
875   - int no_zeros = ((compress && (m_num_inx == m_rows)) || !compress);
  928 + int no_zeros = ((compress && (matrix.m_num_inx == matrix.m_rows)) || !compress);
876 929  
877 930 // Case A: no zero-row/col compression needed
878 931 if (no_zeros) {
879 932  
880 933 // Loop over upper triangle (col >= row)
881 934 double diag = 0.0;
882   - for (int row = 0; row < m_rows; ++row) {
883   - double* ptr = m_data + m_colstart[row];
884   - for (int col = row; col < m_cols; ++col, ++ptr) {
885   - double sum = *ptr; // sum = M(row,col)
  935 + for (int row = 0; row < matrix.m_rows; ++row) {
  936 + double* ptr = matrix.m_data + matrix.m_colstart[row];
  937 + for (int col = row; col < matrix.m_cols; ++col, ++ptr) {
  938 + // sum = M(row,col)
  939 + double sum = *ptr;
886 940 for (int k = 0; k < row; ++k) {
887   - int offset = m_colstart[k] - k; // is always positive
888   - sum -= m_data[offset+row] * m_data[offset+col]; // sum -= M(row,k)*M(col,k)
  941 + int offset = matrix.m_colstart[k] - k; // is always positive
  942 + sum -= matrix.m_data[offset+row] * matrix.m_data[offset+col]; // sum -= M(row,k)*M(col,k)
889 943 }
890 944 if (row == col) {
891 945 if (sum <= 0.0) {
892 946 throw GException::matrix_not_pos_definite(G_CHOL_DECOMP, row, sum);
893 947 }
894   - *ptr = sqrt(sum); // M(row,row) = sqrt(sum)
  948 + *ptr = std::sqrt(sum); // M(row,row) = sqrt(sum)
895 949 diag = 1.0/(*ptr);
896 950 }
897   - else
898   - *ptr = sum*diag; // M(row,col) = sum/M(row,row)
  951 + else {
  952 + *ptr = sum*diag; // M(row,col) = sum/M(row,row)
  953 + }
899 954 }
900 955 }
901 956 } // endif: there were no zero rows/cols in matrix
902 957  
903 958 // Case B: zero-row/col compression needed
904   - else if (m_num_inx > 0) {
  959 + else if (matrix.m_num_inx > 0) {
905 960  
906 961 // Allocate loop variables and pointers
907 962 int row;
... ... @@ -913,25 +968,27 @@ void GMatrixSymmetric::cholesky_decompose(bool compress)
913 968  
914 969 // Loop over upper triangle (col >= row)
915 970 double diag = 0.0;
916   - for (row = 0, row_ptr = m_inx; row < m_num_inx; ++row, ++row_ptr) {
917   - double* ptr_0 = m_data + m_colstart[*row_ptr] - *row_ptr;
918   - for (col = row, col_ptr = m_inx + row; col < m_num_inx; ++col, ++col_ptr) {
  971 + for (row = 0, row_ptr = matrix.m_inx; row < matrix.m_num_inx; ++row, ++row_ptr) {
  972 + double* ptr_0 = matrix.m_data + matrix.m_colstart[*row_ptr] - *row_ptr;
  973 + for (col = row, col_ptr = matrix.m_inx + row; col < matrix.m_num_inx; ++col, ++col_ptr) {
919 974 double* ptr = ptr_0 + *col_ptr;
920 975 double sum = *ptr; // sum = M(row,col)
921   - for (k = 0, k_ptr = m_inx; k < row; ++k, ++k_ptr) {
922   - int offset = m_colstart[*k_ptr] - *k_ptr; // is always positive
923   - sum -= m_data[offset+*row_ptr] * m_data[offset+*col_ptr];
  976 + for (k = 0, k_ptr = matrix.m_inx; k < row; ++k, ++k_ptr) {
  977 + int offset = matrix.m_colstart[*k_ptr] - *k_ptr; // is always positive
  978 + sum -= matrix.m_data[offset+*row_ptr] *
  979 + matrix.m_data[offset+*col_ptr];
924 980 // sum -= M(row,k)*M(col,k)
925 981 }
926 982 if (*row_ptr == *col_ptr) {
927 983 if (sum <= 0.0) {
928 984 throw GException::matrix_not_pos_definite(G_CHOL_DECOMP, *row_ptr, sum);
929 985 }
930   - *ptr = sqrt(sum); // M(row,row) = sqrt(sum)
  986 + *ptr = std::sqrt(sum); // M(row,row) = sqrt(sum)
931 987 diag = 1.0/(*ptr);
932 988 }
933   - else
934   - *ptr = sum*diag; // M(row,col) = sum/M(row,row)
  989 + else {
  990 + *ptr = sum*diag; // M(row,col) = sum/M(row,row)
  991 + }
935 992 }
936 993 }
937 994 } // endelse: zero-row/col compression needed
... ... @@ -941,8 +998,8 @@ void GMatrixSymmetric::cholesky_decompose(bool compress)
941 998 throw GException::matrix_zero(G_CHOL_DECOMP);
942 999 }
943 1000  
944   - // Return
945   - return;
  1001 + // Return matrix
  1002 + return matrix;
946 1003 }
947 1004  
948 1005  
... ... @@ -962,7 +1019,8 @@ void GMatrixSymmetric::cholesky_decompose(bool compress)
962 1019 * that is produced by 'cholesky_decompose'. Case A operates on a full
963 1020 * matrix, Case B on a zero rows/columns (logically) compressed matrix.
964 1021 ***************************************************************************/
965   -GVector GMatrixSymmetric::cholesky_solver(const GVector& vector, bool compress)
  1022 +GVector GMatrixSymmetric::cholesky_solver(const GVector& vector,
  1023 + bool compress) const
966 1024 {
967 1025 // Raise an exception if the matrix and vector dimensions are not compatible
968 1026 if (m_rows != vector.size()) {
... ... @@ -1041,60 +1099,72 @@ GVector GMatrixSymmetric::cholesky_solver(const GVector&amp; vector, bool compress)
1041 1099  
1042 1100  
1043 1101 /***********************************************************************//**
1044   - * @brief Cholesky invert
  1102 + * @brief Invert matrix using a Cholesky decomposition
1045 1103 *
1046   - * @param[in] compress Use zero-row/column compression (default: true).
  1104 + * @param[in] compress Use zero-row/column compression (defaults to true).
  1105 + * @return Inverted matrix.
1047 1106 *
1048 1107 * @exception GException::matrix_zero
1049 1108 * All matrix elements are zero.
1050 1109 *
1051   - * Inplace matrix inversion using the Cholesky decomposition. Case A
1052   - * operates on a full matrix while Case B operates on a (logically)
1053   - * compressed matrix where all zero rows/columns are skipped.
  1110 + * Inverts the matrix using a Cholesky decomposition.
  1111 + *
  1112 + * The method distinguish two cases. Case A operates on a full matrix while
  1113 + * Case B operates on a (logically) compressed matrix where all zero
  1114 + * rows/columns are skipped.
1054 1115 ***************************************************************************/
1055   -void GMatrixSymmetric::cholesky_invert(bool compress)
  1116 +GMatrixSymmetric GMatrixSymmetric::cholesky_invert(bool compress) const
1056 1117 {
1057 1118 // Generate Cholesky decomposition of matrix
1058   - this->cholesky_decompose(compress);
  1119 + GMatrixSymmetric matrix = cholesky_decompose(compress);
1059 1120  
1060 1121 // Check if zero-row/col compression is needed
1061   - int no_zeros = ((compress && (m_num_inx == m_rows)) || !compress);
  1122 + int no_zeros = ((compress && (matrix.m_num_inx == matrix.m_rows)) || !compress);
1062 1123  
1063 1124 // Case A: no zero-row/col compression needed
1064 1125 if (no_zeros) {
1065 1126  
1066 1127 // Generate inverse of Cholesky decomposition (col>row)
1067   - for (int row = 0; row < m_rows; ++row) {
1068   - double* ptr = m_data + m_colstart[row];
1069   - *ptr = 1.0/(*ptr); // M(row,row) = 1/M(row,row)
1070   - for (int col = row+1; col < m_cols; ++col) {
  1128 + for (int row = 0; row < matrix.m_rows; ++row) {
  1129 +
  1130 + // M(row,row) = 1/M(row,row)
  1131 + double* ptr = matrix.m_data + matrix.m_colstart[row];
  1132 + *ptr = 1.0/(*ptr);
  1133 +
  1134 + for (int col = row+1; col < matrix.m_cols; ++col) {
  1135 +
  1136 + // sum -= M(col,k)*M(k,row)
1071 1137 double sum = 0.0;
1072   - double* ptr1 = m_data + col - row;
  1138 + double* ptr1 = matrix.m_data + col - row;
1073 1139 double* ptr2 = ptr;
1074 1140 for (int k = row; k < col; ++k) {
1075   - sum -= *(ptr1-- + m_colstart[k]) * *ptr2++; // sum -= M(col,k)*M(k,row)
  1141 + sum -= *(ptr1-- + matrix.m_colstart[k]) * *ptr2++;
1076 1142 }
1077   - *(ptr+col-row) = sum/m_data[m_colstart[col]]; // M(col,row) = sum/M(col,col)
  1143 +
  1144 + // M(col,row) = sum/M(col,col)
  1145 + *(ptr+col-row) = sum/matrix.m_data[matrix.m_colstart[col]];
1078 1146 }
1079 1147 }
1080 1148  
1081 1149 // Matrix multiplication (col>=row)
1082   - for (int row = 0; row < m_rows; ++row) {
1083   - double* ptr = m_data + m_colstart[row];
1084   - for (int col = row; col < m_cols; ++col) {
  1150 + for (int row = 0; row < matrix.m_rows; ++row) {
  1151 + double* ptr = matrix.m_data + matrix.m_colstart[row];
  1152 + for (int col = row; col < matrix.m_cols; ++col) {
  1153 + // sum += M(row,k)*M(k,col)
1085 1154 double sum = 0.0;
1086 1155 double* ptr1 = ptr + col - row;
1087   - double* ptr2 = m_data + m_colstart[col];
1088   - for (int k = col; k < m_cols; ++k) {
1089   - sum += *ptr1++ * *ptr2++; // sum += M(row,k)*M(k,col)
  1156 + double* ptr2 = matrix.m_data + matrix.m_colstart[col];
  1157 + for (int k = col; k < matrix.m_cols; ++k) {
  1158 + sum += *ptr1++ * *ptr2++;
1090 1159 }
1091   - *(ptr+col-row) = sum; // M(row,col) = sum
  1160 + // M(row,col) = sum
  1161 + *(ptr+col-row) = sum;
1092 1162 }
1093 1163 }
1094 1164 } // endif: no zero-row/col compression needed
1095 1165  
1096 1166 // Case B: zero-row/col compression needed
1097   - else if (m_num_inx > 0) {
  1167 + else if (matrix.m_num_inx > 0) {
1098 1168  
1099 1169 // Allocate loop variables and pointers
1100 1170 int row;
... ... @@ -1105,33 +1175,51 @@ void GMatrixSymmetric::cholesky_invert(bool compress)
1105 1175 int* k_ptr;
1106 1176  
1107 1177 // Generate inverse of Cholesky decomposition (col>row)
1108   - for (row = 0, row_ptr = m_inx; row < m_num_inx; ++row, ++row_ptr) {
1109   - double* ptr_diag = m_data + m_colstart[*row_ptr];
  1178 + for (row = 0, row_ptr = matrix.m_inx;
  1179 + row < matrix.m_num_inx; ++row, ++row_ptr) {
  1180 +
  1181 + // M(row,row) = 1/M(row,row)
  1182 + double* ptr_diag = matrix.m_data + matrix.m_colstart[*row_ptr];
1110 1183 double* ptr_2 = ptr_diag - *row_ptr;
1111   - *ptr_diag = 1.0/(*ptr_diag); // M(row,row) = 1/M(row,row)
1112   - for (col = row+1, col_ptr = m_inx+row+1; col < m_num_inx; ++col, ++col_ptr) {
  1184 + *ptr_diag = 1.0/(*ptr_diag);
  1185 +
  1186 + for (col = row+1, col_ptr = matrix.m_inx+row+1;
  1187 + col < matrix.m_num_inx; ++col, ++col_ptr) {
  1188 +
  1189 + // sum -= M(col,k)*M(k,row)
1113 1190 double sum = 0.0;
1114   - double* ptr_1 = m_data + *col_ptr;
1115   - for (k = row, k_ptr = m_inx+row; k < col; ++k, ++k_ptr) {
1116   - sum -= *(ptr_1 + m_colstart[*k_ptr] - *k_ptr) *
1117   - *(ptr_2 + *k_ptr); // sum -= M(col,k)*M(k,row)
  1191 + double* ptr_1 = matrix.m_data + *col_ptr;
  1192 + for (k = row, k_ptr = matrix.m_inx+row;
  1193 + k < col; ++k, ++k_ptr) {
  1194 + sum -= *(ptr_1 + matrix.m_colstart[*k_ptr] - *k_ptr) *
  1195 + *(ptr_2 + *k_ptr);
1118 1196 }
1119   - *(ptr_2 + *col_ptr) = sum/m_data[m_colstart[*col_ptr]];
1120   - // M(col,row) = sum/M(col,col)
  1197 +
  1198 + // M(col,row) = sum/M(col,col)
  1199 + *(ptr_2 + *col_ptr) =
  1200 + sum/matrix.m_data[matrix.m_colstart[*col_ptr]];
1121 1201 }
1122 1202 }
1123 1203  
1124 1204 // Matrix multiplication (col>=row)
1125   - for (row = 0, row_ptr = m_inx; row < m_num_inx; ++row, ++row_ptr) {
1126   - double* ptr_diag = m_data + m_colstart[*row_ptr];
  1205 + for (row = 0, row_ptr = matrix.m_inx;
  1206 + row < matrix.m_num_inx; ++row, ++row_ptr) {
  1207 + double* ptr_diag = matrix.m_data + matrix.m_colstart[*row_ptr];
1127 1208 double* ptr_1 = ptr_diag - *row_ptr;
1128   - for (col = row, col_ptr = m_inx+row; col < m_num_inx; ++col, ++col_ptr) {
  1209 +
  1210 + for (col = row, col_ptr = matrix.m_inx+row;
  1211 + col < matrix.m_num_inx; ++col, ++col_ptr) {
  1212 +
  1213 + // sum += M(row,k)*M(k,col)
1129 1214 double sum = 0.0;
1130   - double* ptr_2 = m_data + m_colstart[*col_ptr] - *col_ptr;
1131   - for (k = col, k_ptr = m_inx+col; k < m_num_inx; ++k, ++k_ptr) {
1132   - sum += *(ptr_1 + *k_ptr) * *(ptr_2 + *k_ptr); // sum += M(row,k)*M(k,col)
  1215 + double* ptr_2 = matrix.m_data + matrix.m_colstart[*col_ptr] - *col_ptr;
  1216 + for (k = col, k_ptr = matrix.m_inx+col;
  1217 + k < matrix.m_num_inx; ++k, ++k_ptr) {
  1218 + sum += *(ptr_1 + *k_ptr) * *(ptr_2 + *k_ptr);
1133 1219 }
1134   - *(ptr_1 + *col_ptr) = sum; // M(row,col) = sum
  1220 +
  1221 + // M(row,col) = sum
  1222 + *(ptr_1 + *col_ptr) = sum;
1135 1223 }
1136 1224 }
1137 1225 } // endelse: zero-row/col compression needed
... ... @@ -1141,8 +1229,8 @@ void GMatrixSymmetric::cholesky_invert(bool compress)
1141 1229 throw GException::matrix_zero(G_CHOL_INVERT);
1142 1230 }
1143 1231  
1144   - // Return
1145   - return;
  1232 + // Return matrix
  1233 + return matrix;
1146 1234 }
1147 1235  
1148 1236  
... ...