CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

Vector.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // ---------------------------------------------------------------------------
3 //
4 #ifdef GNUPRAGMA
5 #pragma implementation
6 #endif
7 
8 #include <string.h>
9 
10 #include "CLHEP/Matrix/defs.h"
11 #include "CLHEP/Random/Random.h"
12 #include "CLHEP/Vector/ThreeVector.h"
13 #include "CLHEP/Matrix/Vector.h"
14 #include "CLHEP/Matrix/Matrix.h"
15 
16 #ifdef HEP_DEBUG_INLINE
17 #include "CLHEP/Matrix/Vector.icc"
18 #endif
19 
20 namespace CLHEP {
21 
22 // Simple operation for all elements
23 
24 #define SIMPLE_UOP(OPER) \
25  HepGenMatrix::mIter a=m.begin(); \
26  HepGenMatrix::mIter e=m.begin()+num_size(); \
27  for(;a<e; a++) (*a) OPER t;
28 
29 #define SIMPLE_BOP(OPER) \
30  mIter a=m.begin(); \
31  mcIter b=m2.m.begin(); \
32  mcIter e=m.begin()+num_size(); \
33  for(;a<e; a++, b++) (*a) OPER (*b);
34 
35 #define SIMPLE_TOP(OPER) \
36  HepGenMatrix::mcIter a=m1.m.begin(); \
37  HepGenMatrix::mcIter b=m2.m.begin(); \
38  HepGenMatrix::mIter t=mret.m.begin(); \
39  HepGenMatrix::mcIter e=m1.m.begin()+m1.num_size(); \
40  for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b);
41 
42 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
43  if (r1!=r2 || c1!=c2) { \
44  HepGenMatrix::error("Range error in Vector function " #fun "(1)."); \
45  }
46 
47 #define CHK_DIM_1(c1,r2,fun) \
48  if (c1!=r2) { \
49  HepGenMatrix::error("Range error in Vector function " #fun "(2)."); \
50  }
51 
52 // Constructors. (Default constructors are inlined and in .icc file)
53 
55  : m(p), nrow(p)
56 {
57 }
58 
60  : m(p), nrow(p)
61 {
62  switch (init)
63  {
64  case 0:
65  m.assign(p,0);
66  break;
67 
68  case 1:
69  {
70  mIter e = m.begin() + nrow;
71  for (mIter i=m.begin(); i<e; i++) *i = 1.0;
72  break;
73  }
74 
75  default:
76  error("Vector: initialization must be either 0 or 1.");
77  }
78 }
79 
81  : m(p), nrow(p)
82 {
83  HepGenMatrix::mIter a = m.begin();
84  HepGenMatrix::mIter b = m.begin() + nrow;
85  for(;a<b;a++) *a = r();
86 }
87 
88 
89 //
90 // Destructor
91 //
93 }
94 
96  : m(m1.nrow), nrow(m1.nrow)
97 {
98  m = m1.m;
99 }
100 
101 //
102 // Copy constructor from the class of other precision
103 //
104 
105 
107  : m(m1.nrow), nrow(m1.nrow)
108 {
109  if (m1.num_col() != 1)
110  error("Vector::Vector(Matrix) : Matrix is not Nx1");
111 
112  m = m1.m;
113 }
114 
115 // trivial methods
116 
117 inline int HepVector::num_row() const {return nrow;}
118 inline int HepVector::num_size() const {return nrow;}
119 inline int HepVector::num_col() const { return 1; }
120 
121 // operator()
122 
123 #ifdef MATRIX_BOUND_CHECK
124 inline double & HepVector::operator()(int row, int col)
125 {
126  if( col!=1 || row<1 || row>nrow)
127  error("Range error in HepVector::operator(i,j)");
128 #else
129 inline double & HepVector::operator()(int row, int)
130 {
131 #endif
132 
133  return *(m.begin()+(row-1));
134 }
135 
136 #ifdef MATRIX_BOUND_CHECK
137 inline const double & HepVector::operator()(int row, int col) const
138 {
139  if( col!=1 || row<1 || row>nrow)
140  error("Range error in HepVector::operator(i,j)");
141 #else
142 inline const double & HepVector::operator()(int row, int) const
143 {
144 #endif
145 
146  return *(m.begin()+(row-1));
147 }
148 
149 // Sub matrix
150 
151 HepVector HepVector::sub(int min_row, int max_row) const
152 #ifdef HEP_GNU_OPTIMIZED_RETURN
153 return vret(max_row-min_row+1);
154 {
155 #else
156 {
157  HepVector vret(max_row-min_row+1);
158 #endif
159  if(max_row > num_row())
160  error("HepVector::sub: Index out of range");
161  HepGenMatrix::mIter a = vret.m.begin();
162  HepGenMatrix::mcIter b = m.begin() + min_row - 1;
163  HepGenMatrix::mIter e = vret.m.begin() + vret.num_row();
164  for(;a<e;) *(a++) = *(b++);
165  return vret;
166 }
167 
168 HepVector HepVector::sub(int min_row, int max_row)
169 {
170  HepVector vret(max_row-min_row+1);
171  if(max_row > num_row())
172  error("HepVector::sub: Index out of range");
173  HepGenMatrix::mIter a = vret.m.begin();
174  HepGenMatrix::mIter b = m.begin() + min_row - 1;
175  HepGenMatrix::mIter e = vret.m.begin() + vret.num_row();
176  for(;a<e;) *(a++) = *(b++);
177  return vret;
178 }
179 
180 void HepVector::sub(int row,const HepVector &v1)
181 {
182  if(row <1 || row+v1.num_row()-1 > num_row())
183  error("HepVector::sub: Index out of range");
184  HepGenMatrix::mcIter a = v1.m.begin();
185  HepGenMatrix::mIter b = m.begin() + row - 1;
186  HepGenMatrix::mcIter e = v1.m.begin() + v1.num_row();
187  for(;a<e;) *(b++) = *(a++);
188 }
189 
190 //
191 // Direct sum of two matricies
192 //
193 
195  const HepVector &m2)
196 #ifdef HEP_GNU_OPTIMIZED_RETURN
197  return mret(m1.num_row() + m2.num_row(), 0);
198 {
199 #else
200 {
201  HepVector mret(m1.num_row() + m2.num_row(),
202  0);
203 #endif
204  mret.sub(1,m1);
205  mret.sub(m1.num_row()+1,m2);
206  return mret;
207 }
208 
209 /* -----------------------------------------------------------------------
210  This section contains support routines for matrix.h. This section contains
211  The two argument functions +,-. They call the copy constructor and +=,-=.
212  ----------------------------------------------------------------------- */
214 #ifdef HEP_GNU_OPTIMIZED_RETURN
215  return m2(nrow);
216 {
217 #else
218 {
219  HepVector m2(nrow);
220 #endif
221  HepGenMatrix::mcIter a=m.begin();
222  HepGenMatrix::mIter b=m2.m.begin();
223  HepGenMatrix::mcIter e=m.begin()+num_size();
224  for(;a<e; a++, b++) (*b) = -(*a);
225  return m2;
226 }
227 
228 
229 
231 #ifdef HEP_GNU_OPTIMIZED_RETURN
232  return mret(m2);
233 {
234 #else
235 {
236  HepVector mret(m2);
237 #endif
238  CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),1,+);
239  mret += m1;
240  return mret;
241 }
242 
244 #ifdef HEP_GNU_OPTIMIZED_RETURN
245  return mret(m1);
246 {
247 #else
248 {
249  HepVector mret(m1);
250 #endif
251  CHK_DIM_2(m1.num_row(),m2.num_row(),1,m2.num_col(),+);
252  mret += m2;
253  return mret;
254 }
255 
257 #ifdef HEP_GNU_OPTIMIZED_RETURN
258  return mret(m1.num_row());
259 {
260 #else
261 {
262  HepVector mret(m1.num_row());
263 #endif
264  CHK_DIM_1(m1.num_row(),m2.num_row(),+);
265  SIMPLE_TOP(+)
266  return mret;
267 }
268 
269 //
270 // operator -
271 //
272 
274 #ifdef HEP_GNU_OPTIMIZED_RETURN
275  return mret;
276 {
277 #else
278 {
279  HepVector mret;
280 #endif
281  CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),1,-);
282  mret = m1;
283  mret -= m2;
284  return mret;
285 }
286 
288 #ifdef HEP_GNU_OPTIMIZED_RETURN
289  return mret(m1);
290 {
291 #else
292 {
293  HepVector mret(m1);
294 #endif
295  CHK_DIM_2(m1.num_row(),m2.num_row(),1,m2.num_col(),-);
296  mret -= m2;
297  return mret;
298 }
299 
301 #ifdef HEP_GNU_OPTIMIZED_RETURN
302  return mret(m1.num_row());
303 {
304 #else
305 {
306  HepVector mret(m1.num_row());
307 #endif
308  CHK_DIM_1(m1.num_row(),m2.num_row(),-);
309  SIMPLE_TOP(-)
310  return mret;
311 }
312 
313 /* -----------------------------------------------------------------------
314  This section contains support routines for matrix.h. This file contains
315  The two argument functions *,/. They call copy constructor and then /=,*=.
316  ----------------------------------------------------------------------- */
317 
319 const HepVector &m1,double t)
320 #ifdef HEP_GNU_OPTIMIZED_RETURN
321  return mret(m1);
322 {
323 #else
324 {
325  HepVector mret(m1);
326 #endif
327  mret /= t;
328  return mret;
329 }
330 
331 HepVector operator*(const HepVector &m1,double t)
332 #ifdef HEP_GNU_OPTIMIZED_RETURN
333  return mret(m1);
334 {
335 #else
336 {
337  HepVector mret(m1);
338 #endif
339  mret *= t;
340  return mret;
341 }
342 
343 HepVector operator*(double t,const HepVector &m1)
344 #ifdef HEP_GNU_OPTIMIZED_RETURN
345  return mret(m1);
346 {
347 #else
348 {
349  HepVector mret(m1);
350 #endif
351  mret *= t;
352  return mret;
353 }
354 
356 #ifdef HEP_GNU_OPTIMIZED_RETURN
357  return mret(m1.num_row());
358 {
359 #else
360 {
361  HepVector mret(m1.num_row());
362 #endif
363  CHK_DIM_1(m1.num_col(),m2.num_row(),*);
364  HepGenMatrix::mcIter m1p,m2p,vp;
366  double temp;
367  m3p=mret.m.begin();
368  for(m1p=m1.m.begin();m1p<m1.m.begin()+m1.num_row()*m1.num_col();m1p=m2p)
369  {
370  temp=0;
371  vp=m2.m.begin();
372  m2p=m1p;
373  while(m2p<m1p+m1.num_col())
374  temp+=(*(m2p++))*(*(vp++));
375  *(m3p++)=temp;
376  }
377  return mret;
378 }
379 
381 #ifdef HEP_GNU_OPTIMIZED_RETURN
382  return mret(m1.num_row(),m2.num_col());
383 {
384 #else
385 {
386  HepMatrix mret(m1.num_row(),m2.num_col());
387 #endif
388  CHK_DIM_1(1,m2.num_row(),*);
390  HepMatrix::mcIter m2p;
391  HepMatrix::mIter mrp=mret.m.begin();
392  for(m1p=m1.m.begin();m1p<m1.m.begin()+m1.num_row();m1p++)
393  for(m2p=m2.m.begin();m2p<m2.m.begin()+m2.num_col();m2p++)
394  *(mrp++)=*m1p*(*m2p);
395  return mret;
396 }
397 
398 /* -----------------------------------------------------------------------
399  This section contains the assignment and inplace operators =,+=,-=,*=,/=.
400  ----------------------------------------------------------------------- */
401 
403 {
404  CHK_DIM_2(num_row(),m2.num_row(),num_col(),1,+=);
405  SIMPLE_BOP(+=)
406  return (*this);
407 }
408 
410 {
411  CHK_DIM_2(num_row(),m2.num_row(),1,m2.num_col(),+=);
412  SIMPLE_BOP(+=)
413  return (*this);
414 }
415 
417 {
418  CHK_DIM_1(num_row(),m2.num_row(),+=);
419  SIMPLE_BOP(+=)
420  return (*this);
421 }
422 
424 {
425  CHK_DIM_2(num_row(),m2.num_row(),num_col(),1,-=);
426  SIMPLE_BOP(-=)
427  return (*this);
428 }
429 
431 {
432  CHK_DIM_2(num_row(),m2.num_row(),1,m2.num_col(),-=);
433  SIMPLE_BOP(-=)
434  return (*this);
435 }
436 
438 {
439  CHK_DIM_1(num_row(),m2.num_row(),-=);
440  SIMPLE_BOP(-=)
441  return (*this);
442 }
443 
445 {
446  SIMPLE_UOP(/=)
447  return (*this);
448 }
449 
451 {
452  SIMPLE_UOP(*=)
453  return (*this);
454 }
455 
457 {
458  if(m1.nrow != size_)
459  {
460  size_ = m1.nrow;
461  m.resize(size_);
462  }
463  nrow = m1.nrow;
464  ncol = 1;
465  m = m1.m;
466  return (*this);
467 }
468 
470 {
471  if(m1.nrow != nrow)
472  {
473  nrow = m1.nrow;
474  m.resize(nrow);
475  }
476  m = m1.m;
477  return (*this);
478 }
479 
481 {
482  if (m1.num_col() != 1)
483  error("Vector::operator=(Matrix) : Matrix is not Nx1");
484 
485  if(m1.nrow != nrow)
486  {
487  nrow = m1.nrow;
488  m.resize(nrow);
489  }
490  m = m1.m;
491  return (*this);
492 }
493 
495 {
496  if(nrow != 3)
497  {
498  nrow = 3;
499  m.resize(nrow);
500  }
501  m[0] = v.x();
502  m[1] = v.y();
503  m[2] = v.z();
504  return (*this);
505 }
506 
507 //
508 // Copy constructor from the class of other precision
509 //
510 
511 
512 // Print the Matrix.
513 
514 std::ostream& operator<<(std::ostream &s, const HepVector &q)
515 {
516  s << std::endl;
517 /* Fixed format needs 3 extra characters for field, while scientific needs 7 */
518  int width;
519  if(s.flags() & std::ios::fixed)
520  width = s.precision()+3;
521  else
522  width = s.precision()+7;
523  for(int irow = 1; irow<= q.num_row(); irow++)
524  {
525  s.width(width);
526  s << q(irow) << std::endl;
527  }
528  return s;
529 }
530 
532 #ifdef HEP_GNU_OPTIMIZED_RETURN
533 return mret(1,num_row());
534 {
535 #else
536 {
537  HepMatrix mret(1,num_row());
538 #endif
539  mret.m = m;
540  return mret;
541 }
542 
543 double dot(const HepVector &v1,const HepVector &v2)
544 {
545  if(v1.num_row()!=v2.num_row())
546  HepGenMatrix::error("v1 and v2 need to be the same size in dot(HepVector, HepVector)");
547  double d= 0;
548  HepGenMatrix::mcIter a = v1.m.begin();
549  HepGenMatrix::mcIter b = v2.m.begin();
550  HepGenMatrix::mcIter e = a + v1.num_size();
551  for(;a<e;) d += (*(a++)) * (*(b++));
552  return d;
553 }
554 
556 apply(double (*f)(double, int)) const
557 #ifdef HEP_GNU_OPTIMIZED_RETURN
558 return mret(num_row());
559 {
560 #else
561 {
562  HepVector mret(num_row());
563 #endif
564  HepGenMatrix::mcIter a = m.begin();
565  HepGenMatrix::mIter b = mret.m.begin();
566  for(int ir=1;ir<=num_row();ir++) {
567  *(b++) = (*f)(*(a++), ir);
568  }
569  return mret;
570 }
571 
572 void HepVector::invert(int &) {
573  error("HepVector::invert: You cannot invert a Vector");
574 }
575 
576 HepVector solve(const HepMatrix &a, const HepVector &v)
577 #ifdef HEP_GNU_OPTIMIZED_RETURN
578  return vret(v);
579 {
580 #else
581 {
582  HepVector vret(v);
583 #endif
584  static int max_array = 20;
585  static int *ir = new int [max_array+1];
586 
587  if(a.ncol != a.nrow)
588  HepGenMatrix::error("Matrix::solve Matrix is not NxN");
589  if(a.ncol != v.nrow)
590  HepGenMatrix::error("Matrix::solve Vector has wrong number of rows");
591 
592  int n = a.ncol;
593  if (n > max_array) {
594  delete [] ir;
595  max_array = n;
596  ir = new int [max_array+1];
597  }
598  double det;
599  HepMatrix mt(a);
600  int i = mt.dfact_matrix(det, ir);
601  if (i!=0) {
602  for (i=1;i<=n;i++) vret(i) = 0;
603  return vret;
604  }
605  double s21, s22;
606  int nxch = ir[n];
607  if (nxch!=0) {
608  for (int mm=1;mm<=nxch;mm++) {
609  int ij = ir[mm];
610  i = ij >> 12;
611  int j = ij%4096;
612  double te = vret(i);
613  vret(i) = vret(j);
614  vret(j) = te;
615  }
616  }
617  vret(1) = mt(1,1) * vret(1);
618  if (n!=1) {
619  for (i=2;i<=n;i++) {
620  s21 = -vret(i);
621  for (int j=1;j<i;j++) {
622  s21 += mt(i,j) * vret(j);
623  }
624  vret(i) = -mt(i,i)*s21;
625  }
626  for (i=1;i<n;i++) {
627  int nmi = n-i;
628  s22 = -vret(nmi);
629  for (int j=1;j<=i;j++) {
630  s22 += mt(nmi,n-j+1) * vret(n-j+1);
631  }
632  vret(nmi) = -s22;
633  }
634  }
635  return vret;
636 }
637 
638 } // namespace CLHEP