Justine Devin / gammalib

Commit 316cb2d0411aa781b1a003554a12cc29dc52b3ef

Authored by knodlseder
1 parent bb0dce36

make config.h available

Inline Side-by-side
Showing 4 changed files with 218 additions and 190 deletions
include/GMatrix.hpp
  View file @316cb2d
... ... @@ -20,6 +20,9 @@
20 20 #define GMATRIX_HPP
21 21  
22 22 /* __ Includes ___________________________________________________________ */
  23 +#ifdef HAVE_CONFIG_H
  24 +#include <config.h>
  25 +#endif
23 26 #include "GMatrixBase.hpp"
24 27  
25 28  
... ...
include/GMatrixBase.hpp
  View file @316cb2d
... ... @@ -20,6 +20,9 @@
20 20 #define GMATRIXBASE_HPP
21 21  
22 22 /* __ Includes ___________________________________________________________ */
  23 +#ifdef HAVE_CONFIG_H
  24 +#include <config.h>
  25 +#endif
23 26 #include "GVector.hpp"
24 27  
25 28  
... ...
include/GSparseMatrix.hpp
  View file @316cb2d
1 1 /***************************************************************************
2 2 * GSparseMatrix.hpp - sparse matrix class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2006-2009 by Jurgen Knodlseder *
  4 + * copyright (C) 2006-2010 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software; you can redistribute it and/or modify *
... ... @@ -10,11 +10,23 @@
10 10 * (at your option) any later version. *
11 11 * *
12 12 ***************************************************************************/
  13 +/**
  14 + * @file GSparseMatrix.hpp
  15 + * @brief GSparseMatrix class definition.
  16 + * @author J. Knodlseder
  17 + */
13 18  
14 19 #ifndef GSPARSEMATRIX_HPP
15 20 #define GSPARSEMATRIX_HPP
16 21  
17 22 /* __ Includes ___________________________________________________________ */
  23 +#ifdef HAVE_CONFIG_H
  24 +#include <config.h>
  25 +#endif
  26 +#include <math.h>
  27 +#include <iostream>
  28 +#include "GException.hpp"
  29 +#include <stdlib.h>
18 30 #include "GMatrixBase.hpp"
19 31  
20 32 /* __ Definitions ________________________________________________________ */
... ... @@ -23,125 +35,129 @@
23 35 #define G_SPARSE_MATRIX_DEFAULT_STACK_SIZE 512 // Max # of stack elements
24 36  
25 37  
26   -/***************************************************************************
27   - * GSparseMatrix class definition *
28   - * ----------------------------------------------------------------------- *
29   - * GMatrix implements a sparse matrix storage class. It derives from the *
30   - * abstract base class GMatrixBase. *
  38 +/***********************************************************************//**
  39 + * @class GSparseMatrix
  40 + *
  41 + * @brief GSparseMatrix class interface defintion
  42 + *
  43 + * Implements a sparse matrix storage class. It derives from the abstract
  44 + * base class GMatrixBase.
31 45 ***************************************************************************/
32 46 class GSparseMatrix : public GMatrixBase {
33   - // Friend classes
34   - friend class GSparseSymbolic;
35   - friend class GSparseNumeric;
36   -
37   - // Binary operator friends
38   - friend GSparseMatrix operator* (const double& a, const GSparseMatrix& b);
39   - friend GSparseMatrix operator* (const GSparseMatrix& a, const double& b);
40   - friend GSparseMatrix operator/ (const GSparseMatrix& a, const double& b);
41   -
42   - // I/O friends
43   - friend std::ostream& operator<< (std::ostream& os, const GSparseMatrix& m);
44   -
45   - // Friend functions
46   - friend GSparseMatrix transpose(const GSparseMatrix& m);
47   - friend GSparseMatrix fabs(const GSparseMatrix& m);
48   - friend GSparseMatrix cholesky_decompose(const GSparseMatrix& m, int compress = 1);
49   - friend GSparseMatrix cholesky_invert(const GSparseMatrix& m, int compress = 1);
50   - //
51   - friend GSparseMatrix cs_symperm(const GSparseMatrix& m, const int* pinv);
52   - friend GSparseMatrix cs_transpose(const GSparseMatrix& m, int values);
  47 +
  48 + // Friend classes
  49 + friend class GSparseSymbolic;
  50 + friend class GSparseNumeric;
  51 +
  52 + // Binary operator friends
  53 + friend GSparseMatrix operator* (const double& a, const GSparseMatrix& b);
  54 + friend GSparseMatrix operator* (const GSparseMatrix& a, const double& b);
  55 + friend GSparseMatrix operator/ (const GSparseMatrix& a, const double& b);
  56 +
  57 + // I/O friends
  58 + friend std::ostream& operator<< (std::ostream& os, const GSparseMatrix& m);
  59 +
  60 + // Friend functions
  61 + friend GSparseMatrix transpose(const GSparseMatrix& m);
  62 + friend GSparseMatrix fabs(const GSparseMatrix& m);
  63 + friend GSparseMatrix cholesky_decompose(const GSparseMatrix& m,
  64 + int compress = 1);
  65 + friend GSparseMatrix cholesky_invert(const GSparseMatrix& m,
  66 + int compress = 1);
  67 + //
  68 + friend GSparseMatrix cs_symperm(const GSparseMatrix& m, const int* pinv);
  69 + friend GSparseMatrix cs_transpose(const GSparseMatrix& m, int values);
53 70  
54 71 public:
55   - // Constructors and destructors
56   - GSparseMatrix(int rows, int cols, int elements = 0);
57   -// GSparseMatrix(const GMatrix& m);
58   -// GSparseMatrix(const GSymMatrix& m);
59   - GSparseMatrix(const GSparseMatrix& m);
60   - virtual ~GSparseMatrix();
61   -
62   - // Operators
63   - virtual GSparseMatrix& operator= (const GSparseMatrix& m);
64   - virtual double& operator() (int row, int col);
65   - virtual const double& operator() (int row, int col) const;
66   - virtual GSparseMatrix operator+ (const GSparseMatrix& m) const;
67   - virtual GSparseMatrix operator- (const GSparseMatrix& m) const;
68   - virtual GSparseMatrix operator* (const GSparseMatrix& m) const;
69   - virtual GVector operator* (const GVector& v) const;
70   - virtual int operator== (const GSparseMatrix& m) const;
71   - virtual int operator!= (const GSparseMatrix& m) const;
72   - virtual GSparseMatrix operator- () const;
73   - virtual GSparseMatrix& operator+= (const GSparseMatrix& m);
74   - virtual GSparseMatrix& operator-= (const GSparseMatrix& m);
75   - virtual GSparseMatrix& operator*= (const GSparseMatrix& m);
76   - virtual GSparseMatrix& operator*= (const double& d);
77   - virtual GSparseMatrix& operator/= (const double& d);
78   -
79   - // Matrix methods
80   - virtual void add_col(const GVector& v, int col);
81   - virtual void add_col(const double* values, const int* rows, int number, int col);
82   - virtual void cholesky_decompose(int compress = 1);
83   - virtual GVector cholesky_solver(const GVector& v, int compress = 1);
84   - virtual void cholesky_invert(int compress = 1);
85   - virtual void clear();
86   - virtual GVector extract_row(int row) const;
87   - virtual GVector extract_col(int col) const;
88   - //TBD virtual GSparseMatrix extract_lower_triangle() const;
89   - //TBD virtual GSparseMatrix extract_upper_triangle() const;
90   - virtual double fill() const;
91   - virtual void insert_col(const GVector& v, int col);
92   - virtual void insert_col(const double* values, const int* rows, int number, int col);
93   - virtual double min() const;
94   - virtual double max() const;
95   - virtual void set_mem_block(int block);
96   - virtual void stack_init(int size = 0, int entries = 0);
97   - virtual int stack_push_column(const GVector& v, int col);
98   - virtual int stack_push_column(const double* values, const int* rows, int number, int col);
99   - virtual void stack_flush(void);
100   - virtual void stack_destroy(void);
101   - virtual double sum() const;
102   - virtual void transpose();
103   -
  72 + // Constructors and destructors
  73 + GSparseMatrix(int rows, int cols, int elements = 0);
  74 + GSparseMatrix(const GSparseMatrix& m);
  75 + virtual ~GSparseMatrix(void);
  76 +
  77 + // Operators
  78 + GSparseMatrix& operator= (const GSparseMatrix& m);
  79 + double& operator() (int row, int col);
  80 + const double& operator() (int row, int col) const;
  81 + GSparseMatrix operator+ (const GSparseMatrix& m) const;
  82 + GSparseMatrix operator- (const GSparseMatrix& m) const;
  83 + GSparseMatrix operator* (const GSparseMatrix& m) const;
  84 + GVector operator* (const GVector& v) const;
  85 + int operator== (const GSparseMatrix& m) const;
  86 + int operator!= (const GSparseMatrix& m) const;
  87 + GSparseMatrix operator- () const;
  88 + GSparseMatrix& operator+= (const GSparseMatrix& m);
  89 + GSparseMatrix& operator-= (const GSparseMatrix& m);
  90 + GSparseMatrix& operator*= (const GSparseMatrix& m);
  91 + GSparseMatrix& operator*= (const double& d);
  92 + GSparseMatrix& operator/= (const double& d);
  93 +
  94 + // Matrix methods
  95 + void clear(void);
  96 + void transpose(void);
  97 + void add_col(const GVector& v, int col);
  98 + void add_col(const double* values, const int* rows,
  99 + int number, int col);
  100 + void cholesky_decompose(int compress = 1);
  101 + GVector cholesky_solver(const GVector& v, int compress = 1);
  102 + void cholesky_invert(int compress = 1);
  103 + GVector extract_row(int row) const;
  104 + GVector extract_col(int col) const;
  105 + void insert_col(const GVector& v, int col);
  106 + void insert_col(const double* values, const int* rows,
  107 + int number, int col);
  108 + void set_mem_block(int block);
  109 + void stack_init(int size = 0, int entries = 0);
  110 + int stack_push_column(const GVector& v, int col);
  111 + int stack_push_column(const double* values, const int* rows,
  112 + int number, int col);
  113 + void stack_flush(void);
  114 + void stack_destroy(void);
  115 + double fill(void) const;
  116 + double min(void) const;
  117 + double max(void) const;
  118 + double sum(void) const;
  119 +
104 120 private:
105   - // Private methods
106   - void constructor(int rows, int cols, int elements = 0);
107   - void init_members(void);
108   - void init_stack_members(void);
109   - void copy_members(const GSparseMatrix& m);
110   - void free_members(void);
111   - void free_stack_members(void);
112   - int get_index(int row, int col) const;
113   - void fill_pending(void);
114   - void alloc_elements(int start, int num);
115   - void free_elements(int start, int num);
116   - void remove_zero_row_col(void);
117   - void insert_zero_row_col(int rows, int cols);
118   - void mix_column_prepare(const int* src1_row, int src1_num,
119   - const int* src2_row, int src2_num,
120   - int* num_1, int* num_2, int* num_mix);
121   - void mix_column(const double* src1_data, const int* src1_row, int src1_num,
122   - const double* src2_data, const int* src2_row, int src2_num,
123   - double* dst_data, int* dst_row, int* dst_num);
124   -
125   - // Private data members
126   - int* m_rowinx; // Row-indices of all elements
127   - double m_zero; // The zero element (needed for data access)
128   - double m_fill_val; // Element to be filled
129   - int m_fill_row; // Row into which element needs to be filled
130   - int m_fill_col; // Column into which element needs to be filled
131   - int m_mem_block; // Memory block to be allocated at once
132   - void* m_symbolic; // Holds GSparseSymbolic object after decomposition
133   - void* m_numeric; // Holds GSparseNumeric object after decomposition
134   -
135   - // Fill-stack
136   - int m_stack_max_entries; // Maximum number of entries in the stack
137   - int m_stack_size; // Maximum number of elements in the stack
138   - int m_stack_entries; // Number of entries in the stack
139   - int* m_stack_colinx; // Column index for each entry [m_stack_entries]
140   - int* m_stack_start; // Start in stack for each entry [m_stack_entries+1]
141   - double* m_stack_data; // Stack data [m_stack_size]
142   - int* m_stack_rowinx; // Stack row indices [m_stack_size]
143   - int* m_stack_work; // Stack integer working array [m_cols]
144   - double* m_stack_buffer; // Stack double buffer [m_cols]
  121 + // Private methods
  122 + void constructor(int rows, int cols, int elements = 0);
  123 + void init_members(void);
  124 + void init_stack_members(void);
  125 + void copy_members(const GSparseMatrix& m);
  126 + void free_members(void);
  127 + void free_stack_members(void);
  128 + int get_index(int row, int col) const;
  129 + void fill_pending(void);
  130 + void alloc_elements(int start, int num);
  131 + void free_elements(int start, int num);
  132 + void remove_zero_row_col(void);
  133 + void insert_zero_row_col(int rows, int cols);
  134 + void mix_column_prepare(const int* src1_row, int src1_num,
  135 + const int* src2_row, int src2_num,
  136 + int* num_1, int* num_2, int* num_mix);
  137 + void mix_column(const double* src1_data, const int* src1_row, int src1_num,
  138 + const double* src2_data, const int* src2_row, int src2_num,
  139 + double* dst_data, int* dst_row, int* dst_num);
  140 +
  141 + // Private data members
  142 + int* m_rowinx; //!< Row-indices of all elements
  143 + double m_zero; //!< The zero element (needed for data access)
  144 + double m_fill_val; //!< Element to be filled
  145 + int m_fill_row; //!< Row into which element needs to be filled
  146 + int m_fill_col; //!< Column into which element needs to be filled
  147 + int m_mem_block; //!< Memory block to be allocated at once
  148 + void* m_symbolic; //!< Holds GSparseSymbolic object after decomposition
  149 + void* m_numeric; //!< Holds GSparseNumeric object after decomposition
  150 +
  151 + // Fill-stack
  152 + int m_stack_max_entries; //!< Maximum number of entries in the stack
  153 + int m_stack_size; //!< Maximum number of elements in the stack
  154 + int m_stack_entries; //!< Number of entries in the stack
  155 + int* m_stack_colinx; //!< Column index for each entry [m_stack_entries]
  156 + int* m_stack_start; //!< Start in stack for each entry [m_stack_entries+1]
  157 + double* m_stack_data; //!< Stack data [m_stack_size]
  158 + int* m_stack_rowinx; //!< Stack row indices [m_stack_size]
  159 + int* m_stack_work; //!< Stack integer working array [m_cols]
  160 + double* m_stack_buffer; //!< Stack double buffer [m_cols]
145 161 };
146 162  
147 163  
... ... @@ -152,21 +168,22 @@ private:
152 168 inline
153 169 double& GSparseMatrix::operator() (int row, int col)
154 170 {
155   - #if defined(G_RANGE_CHECK)
156   - if (row >= m_rows || col >= m_cols)
157   - throw GException::out_of_range("GSparseMatrix::operator(int,int)", row, col, m_rows, m_cols);
158   - #endif
159   - fill_pending();
160   - int inx = get_index(row,col);
161   - double* value;
162   - if (inx < 0) {
163   - value = &m_fill_val;
164   - m_fill_row = row;
165   - m_fill_col = col;
166   - }
167   - else
168   - value = &(m_data[inx]);
169   - return *value;
  171 + #if defined(G_RANGE_CHECK)
  172 + if (row >= m_rows || col >= m_cols)
  173 + throw GException::out_of_range("GSparseMatrix::operator(int,int)",
  174 + row, col, m_rows, m_cols);
  175 + #endif
  176 + fill_pending();
  177 + int inx = get_index(row,col);
  178 + double* value;
  179 + if (inx < 0) {
  180 + value = &m_fill_val;
  181 + m_fill_row = row;
  182 + m_fill_col = col;
  183 + }
  184 + else
  185 + value = &(m_data[inx]);
  186 + return *value;
170 187 }
171 188  
172 189 // Matrix element access operator (const version)
... ... @@ -176,75 +193,76 @@ double&amp; GSparseMatrix::operator() (int row, int col)
176 193 inline
177 194 const double& GSparseMatrix::operator() (int row, int col) const
178 195 {
179   - #if defined(G_RANGE_CHECK)
180   - if (row >= m_rows || col >= m_cols)
181   - throw GException::out_of_range("GSparseMatrix::operator(int,int)", row, col, m_rows, m_cols);
182   - #endif
183   - int inx = get_index(row,col);
184   - double* value;
185   - if (inx < 0)
186   - value = (double*)&m_zero;
187   - else if (inx == m_elements)
188   - value = (double*)&m_fill_val;
189   - else
190   - value = (double*)&(m_data[inx]);
191   - return *value;
  196 + #if defined(G_RANGE_CHECK)
  197 + if (row >= m_rows || col >= m_cols)
  198 + throw GException::out_of_range("GSparseMatrix::operator(int,int)",
  199 + row, col, m_rows, m_cols);
  200 + #endif
  201 + int inx = get_index(row,col);
  202 + double* value;
  203 + if (inx < 0)
  204 + value = (double*)&m_zero;
  205 + else if (inx == m_elements)
  206 + value = (double*)&m_fill_val;
  207 + else
  208 + value = (double*)&(m_data[inx]);
  209 + return *value;
192 210 }
193 211  
194 212 // Binary matrix addition
195 213 inline
196 214 GSparseMatrix GSparseMatrix::operator+ (const GSparseMatrix& m) const
197 215 {
198   - GSparseMatrix result = *this;
199   - result += m;
200   - return result;
  216 + GSparseMatrix result = *this;
  217 + result += m;
  218 + return result;
201 219 }
202 220  
203 221 // Binary matrix subtraction
204 222 inline
205 223 GSparseMatrix GSparseMatrix::operator- (const GSparseMatrix& m) const
206 224 {
207   - GSparseMatrix result = *this;
208   - result -= m;
209   - return result;
  225 + GSparseMatrix result = *this;
  226 + result -= m;
  227 + return result;
210 228 }
211 229  
212 230 // Binary matrix multiplication
213 231 inline
214 232 GSparseMatrix GSparseMatrix::operator* (const GSparseMatrix& m) const
215 233 {
216   - GSparseMatrix result = *this;
217   - result *= m;
218   - return result;
  234 + GSparseMatrix result = *this;
  235 + result *= m;
  236 + return result;
219 237 }
220 238  
221 239 // Matrix scaling
222 240 inline
223 241 GSparseMatrix& GSparseMatrix::operator*= (const double& s)
224 242 {
225   - fill_pending();
226   - multiplication(s);
227   - return *this;
  243 + fill_pending();
  244 + multiplication(s);
  245 + return *this;
228 246 }
229 247  
230 248 // Matrix scalar division
231 249 inline
232 250 GSparseMatrix& GSparseMatrix::operator/= (const double& s)
233 251 {
234   - double inverse = 1.0/s;
235   - fill_pending();
236   - multiplication(inverse);
237   - return *this;
  252 + double inverse = 1.0/s;
  253 + fill_pending();
  254 + multiplication(inverse);
  255 + return *this;
238 256 }
239 257  
240 258 // Negation
241 259 inline
242 260 GSparseMatrix GSparseMatrix::operator- ( ) const
243 261 {
244   - GSparseMatrix result = *this;
245   - result.fill_pending();
246   - result.negation();
247   - return result;
  262 + GSparseMatrix result = *this;
  263 + result.fill_pending();
  264 + result.negation();
  265 + return result;
248 266 }
249 267  
250 268  
... ... @@ -255,8 +273,8 @@ GSparseMatrix GSparseMatrix::operator- ( ) const
255 273 inline
256 274 void GSparseMatrix::transpose()
257 275 {
258   - fill_pending();
259   - *this = cs_transpose(*this, 1);
  276 + fill_pending();
  277 + *this = cs_transpose(*this, 1);
260 278 }
261 279  
262 280  
... ... @@ -267,65 +285,65 @@ void GSparseMatrix::transpose()
267 285 inline
268 286 GSparseMatrix operator* (const GSparseMatrix& a, const double& b)
269 287 {
270   - GSparseMatrix result = a;
271   - result.fill_pending();
272   - result *= b;
273   - return result;
  288 + GSparseMatrix result = a;
  289 + result.fill_pending();
  290 + result *= b;
  291 + return result;
274 292 }
275 293  
276 294 // Binary matrix scaling (matrix is right operand)
277 295 inline
278 296 GSparseMatrix operator* (const double& a, const GSparseMatrix& b)
279 297 {
280   - GSparseMatrix result = b;
281   - result.fill_pending();
282   - result *= a;
283   - return result;
  298 + GSparseMatrix result = b;
  299 + result.fill_pending();
  300 + result *= a;
  301 + return result;
284 302 }
285 303  
286 304 // Binary matrix division (matrix is left operand)
287 305 inline
288 306 GSparseMatrix operator/ (const GSparseMatrix& a, const double& b)
289 307 {
290   - GSparseMatrix result = a;
291   - result.fill_pending();
292   - result /= b;
293   - return result;
  308 + GSparseMatrix result = a;
  309 + result.fill_pending();
  310 + result /= b;
  311 + return result;
294 312 }
295 313  
296 314 // Matrix transpose function
297 315 inline
298 316 GSparseMatrix transpose(const GSparseMatrix& m)
299 317 {
300   - GSparseMatrix result = m;
301   - result.transpose();
302   - return result;
  318 + GSparseMatrix result = m;
  319 + result.transpose();
  320 + return result;
303 321 }
304 322  
305 323 // Set memory block size
306 324 inline
307 325 void GSparseMatrix::set_mem_block(int block)
308 326 {
309   - m_mem_block = (block > 0) ? block : 1;
310   - return;
  327 + m_mem_block = (block > 0) ? block : 1;
  328 + return;
311 329 }
312 330  
313 331 // Cholesky decomposition
314 332 inline
315 333 GSparseMatrix cholesky_decompose(const GSparseMatrix& m, int compress)
316 334 {
317   - GSparseMatrix result = m;
318   - result.cholesky_decompose(compress);
319   - return result;
  335 + GSparseMatrix result = m;
  336 + result.cholesky_decompose(compress);
  337 + return result;
320 338 }
321 339  
322 340 // Matrix inversion using Cholesky decomposition
323 341 inline
324 342 GSparseMatrix cholesky_invert(const GSparseMatrix& m, int compress)
325 343 {
326   - GSparseMatrix result = m;
327   - result.cholesky_invert(compress);
328   - return result;
  344 + GSparseMatrix result = m;
  345 + result.cholesky_invert(compress);
  346 + return result;
329 347 }
330 348  
331 349  
... ...
include/GVector.hpp
  View file @316cb2d
... ... @@ -20,9 +20,13 @@
20 20 #define GVECTOR_HPP
21 21  
22 22 /* __ Includes ___________________________________________________________ */
  23 +#ifdef HAVE_CONFIG_H
  24 +#include <config.h>
  25 +#endif
  26 +#include <math.h>
  27 +#include <iostream>
  28 +#include <stdlib.h>
23 29 #include "GException.hpp"
24   -#include <math.h> // Mathematics functions
25   -#include <iostream> // std::ostream
26 30  
27 31  
28 32 /***********************************************************************//**
... ...