Francois Brun / gammalib

Commit 40789616c3747e153e2decefa41ae04e3a055e6a

Authored by knodlseder
1 parent 40bd93ae

implement void constructor and allow for zero-element vector

Inline Side-by-side
Showing 2 changed files with 507 additions and 374 deletions
include/GVector.hpp
  View file @4078961
1 1 /***************************************************************************
2 2 * GVector.hpp - vector class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright : (C) 2006 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,6 +10,11 @@
10 10 * (at your option) any later version. *
11 11 * *
12 12 ***************************************************************************/
  13 +/**
  14 + * @file GVector.hpp
  15 + * @brief GVector class definition.
  16 + * @author J. Knodlseder
  17 + */
13 18  
14 19 #ifndef GVECTOR_HPP
15 20 #define GVECTOR_HPP
... ... @@ -17,91 +22,101 @@
17 22 /* __ Includes ___________________________________________________________ */
18 23 #include "GException.hpp"
19 24 #include <math.h> // Mathematics functions
20   -#include <stddef.h> // NULL
21 25 #include <iostream> // std::ostream
22 26  
23 27  
24   -/***************************************************************************
25   - * GVector class definition *
  28 +/***********************************************************************//**
  29 + * @class GVector
  30 + *
  31 + * @brief GVector class interface defintion
26 32 ***************************************************************************/
27 33 class GVector {
28   - // Friend classes
29   - friend class GMatrix;
30   - friend class GSymMatrix;
31   - friend class GSparseMatrix;
32   -
33   - // Operator friends
34   - friend GVector operator+ (const GVector &a, const GVector &b);
35   - friend GVector operator+ (const GVector &a, const double &b);
36   - friend GVector operator+ (const double &a, const GVector &b);
37   - friend GVector operator- (const GVector &a, const GVector &b);
38   - friend GVector operator- (const GVector &a, const double &b);
39   - friend GVector operator- (const double &a, const GVector &b);
40   - friend double operator* (const GVector &a, const GVector &b);
41   - friend GVector operator* (const GVector &a, const double &b);
42   - friend GVector operator* (const double &a, const GVector &b);
43   - friend GVector operator/ (const GVector &a, const double &b);
44   - friend int operator== (const GVector &a, const GVector &b);
45   - friend int operator!= (const GVector &a, const GVector &b);
46   -
47   - // I/O friends
48   - friend std::ostream& operator<< (std::ostream& os, const GVector& v);
49   -
50   - // Friend functions
51   - friend GVector cross(const GVector &a, const GVector &b);
52   - friend double norm(const GVector &v);
53   - friend double min(const GVector &v);
54   - friend double max(const GVector &v);
55   - friend double sum(const GVector &v);
56   - friend GVector perm(const GVector &v, const int *p);
57   - friend GVector iperm(const GVector &v, const int *p);
58   - friend GVector acos(const GVector &v);
59   - friend GVector acosh(const GVector &v);
60   - friend GVector asin(const GVector &v);
61   - friend GVector asinh(const GVector &v);
62   - friend GVector atan(const GVector &v);
63   - friend GVector atanh(const GVector &v);
64   - friend GVector cos(const GVector &v);
65   - friend GVector cosh(const GVector &v);
66   - friend GVector exp(const GVector &v);
67   - friend GVector fabs(const GVector &v);
68   - friend GVector log(const GVector &v);
69   - friend GVector log10(const GVector &v);
70   - friend GVector sin(const GVector &v);
71   - friend GVector sinh(const GVector &v);
72   - friend GVector sqrt(const GVector &v);
73   - friend GVector tan(const GVector &v);
74   - friend GVector tanh(const GVector &v);
75   - friend GVector pow(const GVector &v, const double &q);
  34 +
  35 + // Friend classes
  36 + friend class GMatrix;
  37 + friend class GSymMatrix;
  38 + friend class GSparseMatrix;
  39 +
  40 + // Operator friends
  41 + friend GVector operator+ (const GVector &a, const GVector &b);
  42 + friend GVector operator+ (const GVector &a, const double &b);
  43 + friend GVector operator+ (const double &a, const GVector &b);
  44 + friend GVector operator- (const GVector &a, const GVector &b);
  45 + friend GVector operator- (const GVector &a, const double &b);
  46 + friend GVector operator- (const double &a, const GVector &b);
  47 + friend double operator* (const GVector &a, const GVector &b);
  48 + friend GVector operator* (const GVector &a, const double &b);
  49 + friend GVector operator* (const double &a, const GVector &b);
  50 + friend GVector operator/ (const GVector &a, const double &b);
  51 + friend int operator== (const GVector &a, const GVector &b);
  52 + friend int operator!= (const GVector &a, const GVector &b);
  53 +
  54 + // I/O friends
  55 + friend std::ostream& operator<< (std::ostream& os, const GVector& v);
  56 +
  57 + // Friend functions
  58 + friend GVector cross(const GVector &a, const GVector &b);
  59 + friend double norm(const GVector &v);
  60 + friend double min(const GVector &v);
  61 + friend double max(const GVector &v);
  62 + friend double sum(const GVector &v);
  63 + friend GVector perm(const GVector &v, const int *p);
  64 + friend GVector iperm(const GVector &v, const int *p);
  65 + friend GVector acos(const GVector &v);
  66 + friend GVector acosh(const GVector &v);
  67 + friend GVector asin(const GVector &v);
  68 + friend GVector asinh(const GVector &v);
  69 + friend GVector atan(const GVector &v);
  70 + friend GVector atanh(const GVector &v);
  71 + friend GVector cos(const GVector &v);
  72 + friend GVector cosh(const GVector &v);
  73 + friend GVector exp(const GVector &v);
  74 + friend GVector fabs(const GVector &v);
  75 + friend GVector log(const GVector &v);
  76 + friend GVector log10(const GVector &v);
  77 + friend GVector sin(const GVector &v);
  78 + friend GVector sinh(const GVector &v);
  79 + friend GVector sqrt(const GVector &v);
  80 + friend GVector tan(const GVector &v);
  81 + friend GVector tanh(const GVector &v);
  82 + friend GVector pow(const GVector &v, const double &q);
76 83  
77 84 public:
78   - // Constructors and destructors
79   - explicit GVector(int num);
80   - GVector(const GVector& v);
81   - ~GVector();
82   -
83   - // Vector element access operators
84   - double& operator() (int inx);
85   - const double& operator() (int inx) const;
86   -
87   - // Vector operators
88   - GVector& operator= (const GVector& v);
89   - GVector& operator+= (const GVector& v);
90   - GVector& operator-= (const GVector& v);
91   - GVector& operator= (const double& v);
92   - GVector& operator+= (const double& v);
93   - GVector& operator-= (const double& v);
94   - GVector& operator*= (const double& v);
95   - GVector& operator/= (const double& v);
96   - GVector operator- () const;
97   -
98   - // Vector functions
99   - int size() const; // Return dimension of vector
100   - int non_zeros() const; // Returns # of non zeros in vector
  85 + // Constructors and destructors
  86 + GVector(void);
  87 + explicit GVector(int num);
  88 + GVector(const GVector& v);
  89 + ~GVector(void);
  90 +
  91 + // Vector element access operators
  92 + double& operator() (int inx);
  93 + const double& operator() (int inx) const;
  94 +
  95 + // Vector operators
  96 + GVector& operator= (const GVector& v);
  97 + GVector& operator+= (const GVector& v);
  98 + GVector& operator-= (const GVector& v);
  99 + GVector& operator= (const double& v);
  100 + GVector& operator+= (const double& v);
  101 + GVector& operator-= (const double& v);
  102 + GVector& operator*= (const double& v);
  103 + GVector& operator/= (const double& v);
  104 + GVector operator- () const;
  105 +
  106 + // Vector methods
  107 + int size() const; // Return dimension of vector
  108 + int non_zeros() const; // Returns # of non zeros in vector
101 109  
102 110 private:
103   - int m_num;
104   - double* m_data;
  111 + // Private methods
  112 + void init_members(void);
  113 + void alloc_members(void);
  114 + void copy_members(const GVector& v);
  115 + void free_members(void);
  116 +
  117 + // Private data area
  118 + int m_num; //!< Number of elements in vector
  119 + double* m_data; //!< Vector array
105 120 };
106 121  
107 122  
... ... @@ -113,108 +128,110 @@ private:
113 128 inline
114 129 double& GVector::operator() (int inx)
115 130 {
116   - #if defined(G_RANGE_CHECK)
117   - if (inx < 0 || inx >= m_num)
118   - throw GException::out_of_range("GVector::operator(int)", inx, m_num);
119   - #endif
120   - return m_data[inx];
  131 + #if defined(G_RANGE_CHECK)
  132 + if (inx < 0 || inx >= m_num)
  133 + throw GException::out_of_range("GVector::operator(int)", inx, m_num);
  134 + #endif
  135 + return m_data[inx];
121 136 }
122 137  
123 138 // Vector element access operator (const version)
124 139 inline
125 140 const double& GVector::operator() (int inx) const
126 141 {
127   - #if defined(G_RANGE_CHECK)
128   - if (inx < 0 || inx >= m_num)
129   - throw GException::out_of_range("GVector::operator(int)", inx, m_num);
130   - #endif
131   - return m_data[inx];
  142 + #if defined(G_RANGE_CHECK)
  143 + if (inx < 0 || inx >= m_num)
  144 + throw GException::out_of_range("GVector::operator(int)", inx, m_num);
  145 + #endif
  146 + return m_data[inx];
132 147 }
133 148  
134 149 // Vector unary addition operator
135 150 inline
136 151 GVector& GVector::operator+= (const GVector& v)
137 152 {
138   - if (m_num != v.m_num)
139   - throw GException::vector_mismatch("GVector::operator+=(GVector)", m_num, v.m_num);
140   - for (int i = 0; i < m_num; ++i)
141   - m_data[i] += v.m_data[i];
142   - return *this;
  153 + if (m_num != v.m_num)
  154 + throw GException::vector_mismatch("GVector::operator+=(GVector)",
  155 + m_num, v.m_num);
  156 + for (int i = 0; i < m_num; ++i)
  157 + m_data[i] += v.m_data[i];
  158 + return *this;
143 159 }
144 160  
145 161 // Vector unary subtraction operator
146 162 inline
147 163 GVector& GVector::operator-= (const GVector& v)
148 164 {
149   - if (m_num != v.m_num)
150   - throw GException::vector_mismatch("GVector::operator-=(GVector)", m_num, v.m_num);
151   - for (int i = 0; i < m_num; ++i)
152   - m_data[i] -= v.m_data[i];
153   - return *this;
  165 + if (m_num != v.m_num)
  166 + throw GException::vector_mismatch("GVector::operator-=(GVector)",
  167 + m_num, v.m_num);
  168 + for (int i = 0; i < m_num; ++i)
  169 + m_data[i] -= v.m_data[i];
  170 + return *this;
154 171 }
155 172  
156 173 // Scalar assignment operator
157 174 inline
158 175 GVector& GVector::operator= (const double& v)
159 176 {
160   - for (int i = 0; i < m_num; ++i)
161   - m_data[i] = v;
162   - return *this;
  177 + for (int i = 0; i < m_num; ++i)
  178 + m_data[i] = v;
  179 + return *this;
163 180 }
164 181  
165 182 // Scalar unary addition operator
166 183 inline
167 184 GVector& GVector::operator+= (const double& v)
168 185 {
169   - for (int i = 0; i < m_num; ++i)
170   - m_data[i] += v;
171   - return *this;
  186 + for (int i = 0; i < m_num; ++i)
  187 + m_data[i] += v;
  188 + return *this;
172 189 }
173 190  
174 191 // Scalar unary subtraction operator
175 192 inline
176 193 GVector& GVector::operator-= (const double& v)
177 194 {
178   - for (int i = 0; i < m_num; ++i)
179   - m_data[i] -= v;
180   - return *this;
  195 + for (int i = 0; i < m_num; ++i)
  196 + m_data[i] -= v;
  197 + return *this;
181 198 }
182 199  
183 200 // Scalar unary multiplication operator
184 201 inline
185 202 GVector& GVector::operator*= (const double& v)
186 203 {
187   - for (int i = 0; i < m_num; ++i)
188   - m_data[i] *= v;
189   - return *this;
  204 + for (int i = 0; i < m_num; ++i)
  205 + m_data[i] *= v;
  206 + return *this;
190 207 }
191 208  
192 209 // Scalar unary division operator
193 210 inline
194 211 GVector& GVector::operator/= (const double& v)
195 212 {
196   - for (int i = 0; i < m_num; ++i)
197   - m_data[i] /= v;
198   - return *this;
  213 + for (int i = 0; i < m_num; ++i)
  214 + m_data[i] /= v;
  215 + return *this;
199 216 }
200 217  
201 218 // Unary minus operator
202 219 inline
203 220 GVector GVector::operator- ( ) const
204 221 {
205   - GVector result = *this;
206   - for (int i = 0; i < m_num; ++i)
207   - result.m_data[i] = -result.m_data[i];
208   - return result;
  222 + GVector result = *this;
  223 + for (int i = 0; i < m_num; ++i)
  224 + result.m_data[i] = -result.m_data[i];
  225 + return result;
209 226 }
210 227  
211 228 // Return size of vector
212 229 inline
213   -int GVector::size() const
  230 +int GVector::size() const
214 231 {
215   - return m_num;
  232 + return m_num;
216 233 }
217   -
  234 +
218 235  
219 236 /***************************************************************************
220 237 * Inline friends *
... ... @@ -223,114 +240,115 @@ int GVector::size() const
223 240 inline
224 241 GVector operator+ (const GVector& a, const GVector& b)
225 242 {
226   - GVector result = a;
227   - result += b;
228   - return result;
  243 + GVector result = a;
  244 + result += b;
  245 + return result;
229 246 }
230 247  
231 248 // Vector + Scalar
232 249 inline
233 250 GVector operator+ (const GVector& a, const double& b)
234 251 {
235   - GVector result = a;
236   - result += b;
237   - return result;
  252 + GVector result = a;
  253 + result += b;
  254 + return result;
238 255 }
239 256  
240 257 // Scalar + Vector
241 258 inline
242 259 GVector operator+ (const double& a, const GVector& b)
243 260 {
244   - GVector result = b;
245   - result += a;
246   - return result;
  261 + GVector result = b;
  262 + result += a;
  263 + return result;
247 264 }
248 265  
249 266 // Vector - Vector
250 267 inline
251 268 GVector operator- (const GVector& a, const GVector& b)
252 269 {
253   - GVector result = a;
254   - result -= b;
255   - return result;
  270 + GVector result = a;
  271 + result -= b;
  272 + return result;
256 273 }
257 274  
258 275 // Vector - Scalar
259 276 inline
260 277 GVector operator- (const GVector& a, const double& b)
261 278 {
262   - GVector result = a;
263   - result -= b;
264   - return result;
  279 + GVector result = a;
  280 + result -= b;
  281 + return result;
265 282 }
266 283  
267 284 // Scalar - Vector
268 285 inline
269 286 GVector operator- (const double& a, const GVector& b)
270 287 {
271   - GVector result = -b;
272   - result += a;
273   - return result;
  288 + GVector result = -b;
  289 + result += a;
  290 + return result;
274 291 }
275 292  
276 293 // Scalar product
277 294 inline
278 295 double operator* (const GVector& a, const GVector& b)
279 296 {
280   - if (a.m_num != b.m_num)
281   - throw GException::vector_mismatch("operator*(GVector, GVector)", a.m_num, b.m_num);
282   - double result = 0.0;
283   - for (int i = 0; i < a.m_num; ++i)
284   - result += (a.m_data[i] * b.m_data[i]);
285   - return result;
  297 + if (a.m_num != b.m_num)
  298 + throw GException::vector_mismatch("operator*(GVector, GVector)",
  299 + a.m_num, b.m_num);
  300 + double result = 0.0;
  301 + for (int i = 0; i < a.m_num; ++i)
  302 + result += (a.m_data[i] * b.m_data[i]);
  303 + return result;
286 304 }
287 305  
288 306 // Vector * double
289 307 inline
290 308 GVector operator* (const GVector& a, const double& b)
291 309 {
292   - GVector result = a;
293   - result *= b;
294   - return result;
  310 + GVector result = a;
  311 + result *= b;
  312 + return result;
295 313 }
296 314  
297 315 // double * Vector
298 316 inline
299 317 GVector operator* (const double& a, const GVector& b)
300 318 {
301   - GVector result = b;
302   - result *= a;
303   - return result;
  319 + GVector result = b;
  320 + result *= a;
  321 + return result;
304 322 }
305 323  
306 324 // Vector / double
307 325 inline
308 326 GVector operator/ (const GVector& a, const double& b)
309 327 {
310   - GVector result = a;
311   - result /= b;
312   - return result;
  328 + GVector result = a;
  329 + result /= b;
  330 + return result;
313 331 }
314 332  
315 333 // Vector norm
316 334 inline
317 335 double norm(const GVector &v)
318 336 {
319   - double result = 0.0;
320   - for (int i = 0; i < v.m_num; ++i)
321   - result += (v.m_data[i] * v.m_data[i]);
322   - result = (result > 0.0) ? sqrt(result) : 0.0;
323   - return result;
  337 + double result = 0.0;
  338 + for (int i = 0; i < v.m_num; ++i)
  339 + result += (v.m_data[i] * v.m_data[i]);
  340 + result = (result > 0.0) ? sqrt(result) : 0.0;
  341 + return result;
324 342 }
325 343  
326 344 // Vector minimum
327 345 inline
328 346 double min(const GVector &v)
329 347 {
330   - double result = v.m_data[0];
331   - for (int i = 1; i < v.m_num; ++i) {
332   - if (v.m_data[i] < result)
333   - result = v.m_data[i];
  348 + double result = v.m_data[0];
  349 + for (int i = 1; i < v.m_num; ++i) {
  350 + if (v.m_data[i] < result)
  351 + result = v.m_data[i];
334 352 }
335 353 return result;
336 354 }
... ... @@ -339,220 +357,220 @@ double min(const GVector &amp;v)
339 357 inline
340 358 double max(const GVector &v)
341 359 {
342   - double result = v.m_data[0];
343   - for (int i = 1; i < v.m_num; ++i) {
344   - if (v.m_data[i] > result)
345   - result = v.m_data[i];
346   - }
347   - return result;
  360 + double result = v.m_data[0];
  361 + for (int i = 1; i < v.m_num; ++i) {
  362 + if (v.m_data[i] > result)
  363 + result = v.m_data[i];
  364 + }
  365 + return result;
348 366 }
349 367  
350 368 // Vector sum
351 369 inline
352 370 double sum(const GVector &v)
353 371 {
354   - double result = 0.0;
355   - for (int i = 0; i < v.m_num; ++i)
356   - result += v.m_data[i];
357   - return result;
  372 + double result = 0.0;
  373 + for (int i = 0; i < v.m_num; ++i)
  374 + result += v.m_data[i];
  375 + return result;
358 376 }
359 377  
360 378 // Vector permutation
361 379 inline
362 380 GVector perm(const GVector &v, const int *p)
363 381 {
364   - GVector result(v.m_num);
365   - if (p == NULL)
366   - result = v;
367   - else {
368   - for (int i = 0; i < v.m_num; ++i)
369   - result.m_data[i] = v.m_data[p[i]];
370   - }
371   - return result;
  382 + GVector result(v.m_num);
  383 + if (p == NULL)
  384 + result = v;
  385 + else {
  386 + for (int i = 0; i < v.m_num; ++i)
  387 + result.m_data[i] = v.m_data[p[i]];
  388 + }
  389 + return result;
372 390 }
373 391  
374 392 // Inverse vector permutation
375 393 inline
376 394 GVector iperm(const GVector &v, const int *p)
377 395 {
378   - GVector result(v.m_num);
379   - if (p == NULL)
380   - result = v;
381   - else {
382   - for (int i = 0; i < v.m_num; ++i)
383   - result.m_data[p[i]] = v.m_data[i];
384   - }
385   - return result;
  396 + GVector result(v.m_num);
  397 + if (p == NULL)
  398 + result = v;
  399 + else {
  400 + for (int i = 0; i < v.m_num; ++i)
  401 + result.m_data[p[i]] = v.m_data[i];
  402 + }
  403 + return result;
386 404 }
387 405  
388 406 // Vector acos
389 407 inline
390 408 GVector acos(const GVector &v)
391 409 {
392   - GVector result(v.m_num);
393   - for (int i = 0; i < v.m_num; ++i)
394   - result.m_data[i] = acos(v.m_data[i]);
395   - return result;
  410 + GVector result(v.m_num);
  411 + for (int i = 0; i < v.m_num; ++i)
  412 + result.m_data[i] = acos(v.m_data[i]);
  413 + return result;
396 414 }
397 415  
398 416 // Vector acosh
399 417 inline
400 418 GVector acosh(const GVector &v)
401 419 {
402   - GVector result(v.m_num);
403   - for (int i = 0; i < v.m_num; ++i)
404   - result.m_data[i] = acosh(v.m_data[i]);
405   - return result;
  420 + GVector result(v.m_num);
  421 + for (int i = 0; i < v.m_num; ++i)
  422 + result.m_data[i] = acosh(v.m_data[i]);
  423 + return result;
406 424 }
407 425  
408 426 // Vector asin
409 427 inline
410 428 GVector asin(const GVector &v)
411 429 {
412   - GVector result(v.m_num);
413   - for (int i = 0; i < v.m_num; ++i)
414   - result.m_data[i] = asin(v.m_data[i]);
415   - return result;
  430 + GVector result(v.m_num);
  431 + for (int i = 0; i < v.m_num; ++i)
  432 + result.m_data[i] = asin(v.m_data[i]);
  433 + return result;
416 434 }
417 435  
418 436 // Vector asinh
419 437 inline
420 438 GVector asinh(const GVector &v)
421 439 {
422   - GVector result(v.m_num);
423   - for (int i = 0; i < v.m_num; ++i)
424   - result.m_data[i] = asinh(v.m_data[i]);
425   - return result;
  440 + GVector result(v.m_num);
  441 + for (int i = 0; i < v.m_num; ++i)
  442 + result.m_data[i] = asinh(v.m_data[i]);
  443 + return result;
426 444 }
427 445  
428 446 // Vector atan
429 447 inline
430 448 GVector atan(const GVector &v)
431 449 {
432   - GVector result(v.m_num);
433   - for (int i = 0; i < v.m_num; ++i)
434   - result.m_data[i] = atan(v.m_data[i]);
435   - return result;
  450 + GVector result(v.m_num);
  451 + for (int i = 0; i < v.m_num; ++i)
  452 + result.m_data[i] = atan(v.m_data[i]);
  453 + return result;
436 454 }
437 455  
438 456 // Vector atanh
439 457 inline
440 458 GVector atanh(const GVector &v)
441 459 {
442   - GVector result(v.m_num);
443   - for (int i = 0; i < v.m_num; ++i)
444   - result.m_data[i] = atanh(v.m_data[i]);
445   - return result;
  460 + GVector result(v.m_num);
  461 + for (int i = 0; i < v.m_num; ++i)
  462 + result.m_data[i] = atanh(v.m_data[i]);
  463 + return result;
446 464 }
447 465  
448 466 // Vector cos
449 467 inline
450 468 GVector cos(const GVector &v)
451 469 {
452   - GVector result(v.m_num);
453   - for (int i = 0; i < v.m_num; ++i)
454   - result.m_data[i] = cos(v.m_data[i]);
455   - return result;
  470 + GVector result(v.m_num);
  471 + for (int i = 0; i < v.m_num; ++i)
  472 + result.m_data[i] = cos(v.m_data[i]);
  473 + return result;
456 474 }
457 475  
458 476 // Vector cosh
459 477 inline
460 478 GVector cosh(const GVector &v)
461 479 {
462   - GVector result(v.m_num);
463   - for (int i = 0; i < v.m_num; ++i)
464   - result.m_data[i] = cosh(v.m_data[i]);
465   - return result;
  480 + GVector result(v.m_num);
  481 + for (int i = 0; i < v.m_num; ++i)
  482 + result.m_data[i] = cosh(v.m_data[i]);
  483 + return result;
466 484 }
467 485  
468 486 // Vector exp
469 487 inline
470 488 GVector exp(const GVector &v)
471 489 {
472   - GVector result(v.m_num);
473   - for (int i = 0; i < v.m_num; ++i)
474   - result.m_data[i] = exp(v.m_data[i]);
475   - return result;
  490 + GVector result(v.m_num);
  491 + for (int i = 0; i < v.m_num; ++i)
  492 + result.m_data[i] = exp(v.m_data[i]);
  493 + return result;
476 494 }
477 495  
478 496 // Vector fabs
479 497 inline
480 498 GVector fabs(const GVector &v)
481 499 {
482   - GVector result(v.m_num);
483   - for (int i = 0; i < v.m_num; ++i)
484   - result.m_data[i] = fabs(v.m_data[i]);
485   - return result;
  500 + GVector result(v.m_num);
  501 + for (int i = 0; i < v.m_num; ++i)
  502 + result.m_data[i] = fabs(v.m_data[i]);
  503 + return result;
486 504 }
487 505  
488 506 // Vector log
489 507 inline
490 508 GVector log(const GVector &v)
491 509 {
492   - GVector result(v.m_num);
493   - for (int i = 0; i < v.m_num; ++i)
494   - result.m_data[i] = log(v.m_data[i]);
495   - return result;
  510 + GVector result(v.m_num);
  511 + for (int i = 0; i < v.m_num; ++i)
  512 + result.m_data[i] = log(v.m_data[i]);
  513 + return result;
496 514 }
497 515  
498 516 // Vector log10
499 517 inline
500 518 GVector log10(const GVector &v)
501 519 {
502   - GVector result(v.m_num);
503   - for (int i = 0; i < v.m_num; ++i)
504   - result.m_data[i] = log10(v.m_data[i]);
505   - return result;
  520 + GVector result(v.m_num);
  521 + for (int i = 0; i < v.m_num; ++i)
  522 + result.m_data[i] = log10(v.m_data[i]);
  523 + return result;
506 524 }
507 525  
508 526 // Vector sin
509 527 inline
510 528 GVector sin(const GVector &v)
511 529 {
512   - GVector result(v.m_num);
513   - for (int i = 0; i < v.m_num; ++i)
514   - result.m_data[i] = sin(v.m_data[i]);
515   - return result;
  530 + GVector result(v.m_num);
  531 + for (int i = 0; i < v.m_num; ++i)
  532 + result.m_data[i] = sin(v.m_data[i]);
  533 + return result;
516 534 }
517 535  
518 536 // Vector sinh
519 537 inline
520 538 GVector sinh(const GVector &v)
521 539 {
522   - GVector result(v.m_num);
523   - for (int i = 0; i < v.m_num; ++i)
524   - result.m_data[i] = sinh(v.m_data[i]);
525   - return result;
  540 + GVector result(v.m_num);
  541 + for (int i = 0; i < v.m_num; ++i)
  542 + result.m_data[i] = sinh(v.m_data[i]);
  543 + return result;
526 544 }
527 545  
528 546 // Vector sqrt
529 547 inline
530 548 GVector sqrt(const GVector &v)
531 549 {
532   - GVector result(v.m_num);
533   - for (int i = 0; i < v.m_num; ++i)
534   - result.m_data[i] = sqrt(v.m_data[i]);
535   - return result;
  550 + GVector result(v.m_num);
  551 + for (int i = 0; i < v.m_num; ++i)
  552 + result.m_data[i] = sqrt(v.m_data[i]);
  553 + return result;
536 554 }
537 555  
538 556 // Vector tan
539 557 inline
540 558 GVector tan(const GVector &v)
541 559 {
542   - GVector result(v.m_num);
543   - for (int i = 0; i < v.m_num; ++i)
544   - result.m_data[i] = tan(v.m_data[i]);
545   - return result;
  560 + GVector result(v.m_num);
  561 + for (int i = 0; i < v.m_num; ++i)
  562 + result.m_data[i] = tan(v.m_data[i]);
  563 + return result;
546 564 }
547 565  
548 566 // Vector tanh
549 567 inline
550 568 GVector tanh(const GVector &v)
551 569 {
552   - GVector result(v.m_num);
553   - for (int i = 0; i < v.m_num; ++i)
554   - result.m_data[i] = tanh(v.m_data[i]);
555   - return result;
  570 + GVector result(v.m_num);
  571 + for (int i = 0; i < v.m_num; ++i)
  572 + result.m_data[i] = tanh(v.m_data[i]);
  573 + return result;
556 574 }
557 575  
558 576 // Vector pow
... ... @@ -569,36 +587,36 @@ GVector pow(const GVector &amp;v, const double &amp;a)
569 587 inline
570 588 int operator== (const GVector &a, const GVector &b)
571 589 {
572   - int result = 1;
573   - if (a.m_num == b.m_num) {
574   - for (int i = 0; i < a.m_num; ++i) {
575   - if (a.m_data[i] != b.m_data[i]) {
576   - result = 0;
577   - break;
578   - }
579   - }
580   - }
581   - else
582   - result = 0;
583   - return result;
  590 + int result = 1;
  591 + if (a.m_num == b.m_num) {
  592 + for (int i = 0; i < a.m_num; ++i) {
  593 + if (a.m_data[i] != b.m_data[i]) {
  594 + result = 0;
  595 + break;
  596 + }
  597 + }
  598 + }
  599 + else
  600 + result = 0;
  601 + return result;
584 602 }
585 603  
586 604 // Not equality operator
587 605 inline
588 606 int operator!= (const GVector &a, const GVector &b)
589 607 {
590   - int result = 0;
591   - if (a.m_num == b.m_num) {
592   - for (int i = 0; i < a.m_num; ++i) {
593   - if (a.m_data[i] != b.m_data[i]) {
594   - result = 1;
595   - break;
596   - }
597   - }
598   - }
599   - else
600   - result = 1;
601   - return result;
  608 + int result = 0;
  609 + if (a.m_num == b.m_num) {
  610 + for (int i = 0; i < a.m_num; ++i) {
  611 + if (a.m_data[i] != b.m_data[i]) {
  612 + result = 1;
  613 + break;
  614 + }
  615 + }
  616 + }
  617 + else
  618 + result = 1;
  619 + return result;
602 620 }
603 621  
604 622 #endif /* GVECTOR_HPP */
... ...
src/linalg/GVector.cpp
  View file @4078961
1 1 /***************************************************************************
2 2 * GVector.cpp - vector class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright : (C) 2006 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,71 +10,90 @@
10 10 * (at your option) any later version. *
11 11 * *
12 12 ***************************************************************************/
  13 +/**
  14 + * @file GVector.cpp
  15 + * @brief GVector class implementation.
  16 + * @author J. Knodlseder
  17 + */
13 18  
14 19 /* __ Includes ___________________________________________________________ */
15 20 #include "GVector.hpp"
16 21  
17 22  
  23 +/* __ Method name definitions ____________________________________________ */
  24 +#define G_CROSS "cross(GVector,GVector)"
  25 +
  26 +
18 27 /*==========================================================================
19 28 = =
20 29 = GVector constructors/destructors =
21 30 = =
22 31 ==========================================================================*/
23 32  
24   -/***************************************************************************
25   - * Constructor *
  33 +/***********************************************************************//**
  34 + * @brief Void vector constructor (contains no elements)
  35 + ***************************************************************************/
  36 +GVector::GVector(void)
  37 +{
  38 + // Initialise class members
  39 + init_members();
  40 +
  41 + // Return
  42 + return;
  43 +}
  44 +
  45 +
  46 +/***********************************************************************//**
  47 + * @brief Vector constructor
  48 + *
  49 + * @param[in] num Number of elements in vector
  50 + *
  51 + * Initialises a vector with num elements (values are set to 0)
26 52 ***************************************************************************/
27 53 GVector::GVector(int num)
28 54 {
29   - // Throw exception if requested vector length is zero
30   - if (num == 0)
31   - throw GException::empty("GVector constructor");
32   -
33   - // Allocate vector array. Throw exception if allocation failed
34   - m_data = new double[num];
35   - if (m_data == NULL)
36   - throw GException::mem_alloc("GVector constructor", num);
37   -
38   - // Store vector size and initialize elements to 0.0
39   - m_num = num;
40   - for (int i = 0; i < m_num; ++i)
41   - m_data[i] = 0.0;
42   -
43   - // Return
44   - return;
  55 + // Initialise class members
  56 + init_members();
  57 +
  58 + // Store vector size
  59 + m_num = num;
  60 +
  61 + // Allocate vector (filled with 0)
  62 + alloc_members();
  63 +
  64 + // Return
  65 + return;
45 66 }
46 67  
47 68  
48   -/***************************************************************************
49   - * Copy constructor *
  69 +/***********************************************************************//**
  70 + * @brief Copy constructor
  71 + *
  72 + * @param[in] v Vector from which class should be instantiated.
50 73 ***************************************************************************/
51 74 GVector::GVector(const GVector& v)
52 75 {
53   - // Allocate vector array. Throw exception if allocation failed
54   - m_data = new double[v.m_num];
55   - if (m_data == NULL)
56   - throw GException::mem_alloc("GVector copy constructor", v.m_num);
57   -
58   - // Store vector size and copy elements
59   - m_num = v.m_num;
60   - for (int i = 0; i < m_num; ++i)
61   - m_data[i] = v.m_data[i];
62   -
63   - // Return
64   - return;
  76 + // Initialise class members
  77 + init_members();
  78 +
  79 + // Copy members
  80 + copy_members(v);
  81 +
  82 + // Return
  83 + return;
65 84 }
66 85  
67 86  
68   -/***************************************************************************
69   - * Destructor *
  87 +/***********************************************************************//**
  88 + * @brief Destructor
70 89 ***************************************************************************/
71   -GVector::~GVector()
  90 +GVector::~GVector(void)
72 91 {
73   - // Deallocate only if memory has indeed been allocated
74   - if (m_data != NULL) delete[] m_data;
  92 + // Free members
  93 + free_members();
75 94  
76   - // Return
77   - return;
  95 + // Return
  96 + return;
78 97 }
79 98  
80 99  
... ... @@ -84,96 +103,192 @@ GVector::~GVector()
84 103 = =
85 104 ==========================================================================*/
86 105  
87   -/***************************************************************************
88   - * GVector assignment operator *
  106 +/***********************************************************************//**
  107 + * @brief Assignment operator
  108 + *
  109 + * @param[in] v GVector instance to be assigned
89 110 ***************************************************************************/
90 111 GVector& GVector::operator= (const GVector& v)
91 112 {
92   - // Assign only if not identical
93   - if (this != &v) {
94   -
95   - // First delete old data if it exists
96   - if (m_data != NULL) delete [] m_data;
97   -
98   - // Allocate vector array. Throw exception if allocation failed
99   - m_data = new double[v.m_num];
100   - if (m_data == NULL)
101   - throw GException::mem_alloc("GVector assignment operator", v.m_num);
102   -
103   - // Store vector size and copy elements
104   - m_num = v.m_num;
105   - for (int i = 0; i < m_num; ++i)
106   - m_data[i] = v.m_data[i];
107   - }
108   -
109   - // Return vector
110   - return *this;
  113 + // Execute only if object is not identical
  114 + if (this != &v) {
  115 +
  116 + // Free members
  117 + free_members();
  118 +
  119 + // Initialise private members
  120 + init_members();
  121 +
  122 + // Copy members
  123 + copy_members(v);
  124 +
  125 + } // endif: object was not identical
  126 +
  127 + // Return this object
  128 + return *this;
111 129 }
112 130  
113 131  
114 132 /*==========================================================================
115 133 = =
116   - = GVector member functions =
  134 + = GVector public methods =
117 135 = =
118 136 ==========================================================================*/
119 137  
120   -/***************************************************************************
121   - * Returns number of non-zeros in vector *
  138 +/***********************************************************************//**
  139 + * @brief Returns number of non-zero elements in vector
122 140 ***************************************************************************/
123 141 int GVector::non_zeros() const
124 142 {
125   - // Initialise number of non-zeros
126   - int non_zeros = 0;
127   -
128   - // Gather all non-zero elements
129   - for (int i = 0; i < m_num; ++i) {
130   - if (m_data[i] != 0.0)
131   - non_zeros++;
132   - }
133   -
134   - // Return number of non-zero elements
135   - return non_zeros;
  143 + // Initialise number of non-zeros
  144 + int non_zeros = 0;
  145 +
  146 + // Gather all non-zero elements
  147 + for (int i = 0; i < m_num; ++i) {
  148 + if (m_data[i] != 0.0)
  149 + non_zeros++;
  150 + }
  151 +
  152 + // Return number of non-zero elements
  153 + return non_zeros;
136 154 }
137 155  
138 156  
139 157 /*==========================================================================
140 158 = =
141   - = GVector private functions =
  159 + = GVector private methods =
142 160 = =
143 161 ==========================================================================*/
144 162  
  163 +/***********************************************************************//**
  164 + * @brief Initialise class members
  165 + ***************************************************************************/
  166 +void GVector::init_members(void)
  167 +{
  168 + // Initialise members
  169 + m_num = 0;
  170 + m_data = NULL;
  171 +
  172 + // Return
  173 + return;
  174 +}
  175 +
  176 +
  177 +/***********************************************************************//**
  178 + * @brief Allocate vector
  179 + ***************************************************************************/
  180 +void GVector::alloc_members(void)
  181 +{
  182 + // Continue only if vector has non-zero length
  183 + if (m_num > 0) {
  184 +
  185 + // Allocate vector and initialize elements to 0
  186 + m_data = new double[m_num];
  187 + for (int i = 0; i < m_num; ++i)
  188 + m_data[i] = 0.0;
  189 +
  190 + } // endif: vector had non-zero length
  191 +
  192 + // Return
  193 + return;
  194 +}
  195 +
  196 +
  197 +/***********************************************************************//**
  198 + * @brief Copy class members
  199 + *
  200 + * @param[in] v Vector from which members should be copied.
  201 + ***************************************************************************/
  202 +void GVector::copy_members(const GVector& v)
  203 +{
  204 + // Copy attributes
  205 + m_num = v.m_num;
  206 +
  207 + // Copy elements
  208 + if (m_num > 0) {
  209 + alloc_members();
  210 + for (int i = 0; i < m_num; ++i)
  211 + m_data[i] = v.m_data[i];
  212 + }
  213 +
  214 + // Return
  215 + return;
  216 +}
  217 +
  218 +
  219 +/***********************************************************************//**
  220 + * @brief Delete class members
  221 + ***************************************************************************/
  222 +void GVector::free_members(void)
  223 +{
  224 + // Free memory
  225 + if (m_data != NULL) delete m_data;
  226 +
  227 + // Signal free pointers
  228 + m_data = NULL;
  229 +
  230 + // Return
  231 + return;
  232 +}
  233 +
  234 +
145 235 /*==========================================================================
146 236 = =
147 237 = GVector friends =
148 238 = =
149 239 ==========================================================================*/
150 240  
151   -/***************************************************************************
152   - * GVector output operator *
  241 +/***********************************************************************//**
  242 + * @brief Output operator
  243 + *
  244 + * @param[in] os Output stream.
  245 + * @param[in] v Vector to put in output stream.
153 246 ***************************************************************************/
154 247 std::ostream& operator<< (std::ostream& os, const GVector& v)
155 248 {
156   - for (int i = 0; i < v.m_num; ++i) {
157   - os << v(i);
158   - if (i != v.m_num-1)
159   - os << ", ";
160   - }
161   - return os;
  249 + // Prepend (
  250 + os << "(";
  251 +
  252 + // Put all elements in stream
  253 + for (int i = 0; i < v.m_num; ++i) {
  254 + os << v(i);
  255 + if (i != v.m_num-1)
  256 + os << ", ";
  257 + }
  258 +
  259 + // Append )
  260 + os << ")";
  261 +
  262 + // Return output stream
  263 + return os;
162 264 }
163 265  
164 266  
165   -/***************************************************************************
166   - * GVector cross product *
  267 +/***********************************************************************//**
  268 + * @brief GVector cross product
  269 + *
  270 + * @param[in] a First vector for cross product.
  271 + * @param[in] b Second vector for cross product.
  272 + *
  273 + * Computes the cross product between two 3-element vectors (note that the
  274 + * cross product is only defined for 3-element vectors).
167 275 ***************************************************************************/
168 276 GVector cross(const GVector &a, const GVector &b)
169 277 {
170   - if (a.m_num != b.m_num)
171   - throw GException::vector_mismatch("cross(GVector&, GVector&)", a.m_num, b.m_num);
172   - if (a.m_num != 3)
173   - throw GException::vector_bad_cross_dim("cross(GVector&, GVector&)", a.m_num);
174   - GVector result(3);
175   - result(0) = a.m_data[1]*b.m_data[2] - a.m_data[2]*b.m_data[1];
176   - result(1) = a.m_data[2]*b.m_data[0] - a.m_data[0]*b.m_data[2];
177   - result(2) = a.m_data[0]*b.m_data[1] - a.m_data[1]*b.m_data[0];
178   - return result;
  278 + // Verify that vectors have same dimensions
  279 + if (a.m_num != b.m_num)
  280 + throw GException::vector_mismatch(G_CROSS, a.m_num, b.m_num);
  281 +
  282 + // Verify that vectors have 3 elements
  283 + if (a.m_num != 3)
  284 + throw GException::vector_bad_cross_dim(G_CROSS, a.m_num);
  285 +
  286 + // Compute cross product
  287 + GVector result(3);
  288 + result(0) = a.m_data[1]*b.m_data[2] - a.m_data[2]*b.m_data[1];
  289 + result(1) = a.m_data[2]*b.m_data[0] - a.m_data[0]*b.m_data[2];
  290 + result(2) = a.m_data[0]*b.m_data[1] - a.m_data[1]*b.m_data[0];
  291 +
  292 + // Return result
  293 + return result;
179 294 }
... ...