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=hm2.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=hm1.m.begin(); \
37  HepGenMatrix::mcIter b=hm2.m.begin(); \
38  HepGenMatrix::mIter t=mret.m.begin(); \
39  HepGenMatrix::mcIter e=hm1.m.begin()+hm1.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  : HepGenMatrix(hm1), m(hm1.nrow), nrow(hm1.nrow)
97 {
98  m = hm1.m;
99 }
100 
101 //
102 // Copy constructor from the class of other precision
103 //
104 
105 
107  : m(hm1.nrow), nrow(hm1.nrow)
108 {
109  if (hm1.num_col() != 1)
110  error("Vector::Vector(Matrix) : Matrix is not Nx1");
111 
112  m = hm1.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 &hm2)
196 #ifdef HEP_GNU_OPTIMIZED_RETURN
197  return mret(hm1.num_row() + hm2.num_row(), 0);
198 {
199 #else
200 {
201  HepVector mret(hm1.num_row() + hm2.num_row(),
202  0);
203 #endif
204  mret.sub(1,hm1);
205  mret.sub(hm1.num_row()+1,hm2);
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 hm2(nrow);
216 {
217 #else
218 {
219  HepVector hm2(nrow);
220 #endif
221  HepGenMatrix::mcIter a=m.begin();
222  HepGenMatrix::mIter b=hm2.m.begin();
223  HepGenMatrix::mcIter e=m.begin()+num_size();
224  for(;a<e; a++, b++) (*b) = -(*a);
225  return hm2;
226 }
227 
228 
229 
230 HepVector operator+(const HepMatrix &hm1,const HepVector &hm2)
231 #ifdef HEP_GNU_OPTIMIZED_RETURN
232  return mret(hm2);
233 {
234 #else
235 {
236  HepVector mret(hm2);
237 #endif
238  CHK_DIM_2(hm1.num_row(),hm2.num_row(),hm1.num_col(),1,+);
239  mret += hm1;
240  return mret;
241 }
242 
243 HepVector operator+(const HepVector &hm1,const HepMatrix &hm2)
244 #ifdef HEP_GNU_OPTIMIZED_RETURN
245  return mret(hm1);
246 {
247 #else
248 {
249  HepVector mret(hm1);
250 #endif
251  CHK_DIM_2(hm1.num_row(),hm2.num_row(),1,hm2.num_col(),+);
252  mret += hm2;
253  return mret;
254 }
255 
256 HepVector operator+(const HepVector &hm1,const HepVector &hm2)
257 #ifdef HEP_GNU_OPTIMIZED_RETURN
258  return mret(hm1.num_row());
259 {
260 #else
261 {
262  HepVector mret(hm1.num_row());
263 #endif
264  CHK_DIM_1(hm1.num_row(),hm2.num_row(),+);
265  SIMPLE_TOP(+)
266  return mret;
267 }
268 
269 //
270 // operator -
271 //
272 
273 HepVector operator-(const HepMatrix &hm1,const HepVector &hm2)
274 #ifdef HEP_GNU_OPTIMIZED_RETURN
275  return mret;
276 {
277 #else
278 {
279  HepVector mret;
280 #endif
281  CHK_DIM_2(hm1.num_row(),hm2.num_row(),hm1.num_col(),1,-);
282  mret = hm1;
283  mret -= hm2;
284  return mret;
285 }
286 
287 HepVector operator-(const HepVector &hm1,const HepMatrix &hm2)
288 #ifdef HEP_GNU_OPTIMIZED_RETURN
289  return mret(hm1);
290 {
291 #else
292 {
293  HepVector mret(hm1);
294 #endif
295  CHK_DIM_2(hm1.num_row(),hm2.num_row(),1,hm2.num_col(),-);
296  mret -= hm2;
297  return mret;
298 }
299 
300 HepVector operator-(const HepVector &hm1,const HepVector &hm2)
301 #ifdef HEP_GNU_OPTIMIZED_RETURN
302  return mret(hm1.num_row());
303 {
304 #else
305 {
306  HepVector mret(hm1.num_row());
307 #endif
308  CHK_DIM_1(hm1.num_row(),hm2.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 
318 HepVector operator/(
319 const HepVector &hm1,double t)
320 #ifdef HEP_GNU_OPTIMIZED_RETURN
321  return mret(hm1);
322 {
323 #else
324 {
325  HepVector mret(hm1);
326 #endif
327  mret /= t;
328  return mret;
329 }
330 
331 HepVector operator*(const HepVector &hm1,double t)
332 #ifdef HEP_GNU_OPTIMIZED_RETURN
333  return mret(hm1);
334 {
335 #else
336 {
337  HepVector mret(hm1);
338 #endif
339  mret *= t;
340  return mret;
341 }
342 
343 HepVector operator*(double t,const HepVector &hm1)
344 #ifdef HEP_GNU_OPTIMIZED_RETURN
345  return mret(hm1);
346 {
347 #else
348 {
349  HepVector mret(hm1);
350 #endif
351  mret *= t;
352  return mret;
353 }
354 
355 HepVector operator*(const HepMatrix &hm1,const HepVector &hm2)
356 #ifdef HEP_GNU_OPTIMIZED_RETURN
357  return mret(hm1.num_row());
358 {
359 #else
360 {
361  HepVector mret(hm1.num_row());
362 #endif
363  CHK_DIM_1(hm1.num_col(),hm2.num_row(),*);
364  HepGenMatrix::mcIter hm1p,hm2p,vp;
366  double temp;
367  m3p=mret.m.begin();
368  for(hm1p=hm1.m.begin();hm1p<hm1.m.begin()+hm1.num_row()*hm1.num_col();hm1p=hm2p)
369  {
370  temp=0;
371  vp=hm2.m.begin();
372  hm2p=hm1p;
373  while(hm2p<hm1p+hm1.num_col())
374  temp+=(*(hm2p++))*(*(vp++));
375  *(m3p++)=temp;
376  }
377  return mret;
378 }
379 
380 HepMatrix operator*(const HepVector &hm1,const HepMatrix &hm2)
381 #ifdef HEP_GNU_OPTIMIZED_RETURN
382  return mret(hm1.num_row(),hm2.num_col());
383 {
384 #else
385 {
386  HepMatrix mret(hm1.num_row(),hm2.num_col());
387 #endif
388  CHK_DIM_1(1,hm2.num_row(),*);
390  HepMatrix::mcIter hm2p;
391  HepMatrix::mIter mrp=mret.m.begin();
392  for(hm1p=hm1.m.begin();hm1p<hm1.m.begin()+hm1.num_row();hm1p++)
393  for(hm2p=hm2.m.begin();hm2p<hm2.m.begin()+hm2.num_col();hm2p++)
394  *(mrp++)=*hm1p*(*hm2p);
395  return mret;
396 }
397 
398 /* -----------------------------------------------------------------------
399  This section contains the assignment and inplace operators =,+=,-=,*=,/=.
400  ----------------------------------------------------------------------- */
401 
403 {
404  CHK_DIM_2(num_row(),hm2.num_row(),num_col(),1,+=);
405  SIMPLE_BOP(+=)
406  return (*this);
407 }
408 
410 {
411  CHK_DIM_2(num_row(),hm2.num_row(),1,hm2.num_col(),+=);
412  SIMPLE_BOP(+=)
413  return (*this);
414 }
415 
417 {
418  CHK_DIM_1(num_row(),hm2.num_row(),+=);
419  SIMPLE_BOP(+=)
420  return (*this);
421 }
422 
424 {
425  CHK_DIM_2(num_row(),hm2.num_row(),num_col(),1,-=);
426  SIMPLE_BOP(-=)
427  return (*this);
428 }
429 
431 {
432  CHK_DIM_2(num_row(),hm2.num_row(),1,hm2.num_col(),-=);
433  SIMPLE_BOP(-=)
434  return (*this);
435 }
436 
438 {
439  CHK_DIM_1(num_row(),hm2.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(hm1.nrow != size_)
459  {
460  size_ = hm1.nrow;
461  m.resize(size_);
462  }
463  nrow = hm1.nrow;
464  ncol = 1;
465  m = hm1.m;
466  return (*this);
467 }
468 
470 {
471  if(hm1.nrow != nrow)
472  {
473  nrow = hm1.nrow;
474  m.resize(nrow);
475  }
476  m = hm1.m;
477  return (*this);
478 }
479 
481 {
482  if (hm1.num_col() != 1)
483  error("Vector::operator=(Matrix) : Matrix is not Nx1");
484 
485  if(hm1.nrow != nrow)
486  {
487  nrow = hm1.nrow;
488  m.resize(nrow);
489  }
490  m = hm1.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 &os, const HepVector &q)
515 {
516  os << std::endl;
517 /* Fixed format needs 3 extra characters for field, while scientific needs 7 */
518  int width;
519  if(os.flags() & std::ios::fixed)
520  width = os.precision()+3;
521  else
522  width = os.precision()+7;
523  for(int irow = 1; irow<= q.num_row(); irow++)
524  {
525  os.width(width);
526  os << q(irow) << std::endl;
527  }
528  return os;
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 
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 hmm=1;hmm<=nxch;hmm++) {
609  int ij = ir[hmm];
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
a
@ a
Definition: testCategories.cc:125
b
@ b
Definition: testCategories.cc:125
CLHEP::HepMatrix::num_row
virtual int num_row() const
Definition: Matrix.cc:120
CLHEP::HepGenMatrix
Definition: Matrix/CLHEP/Matrix/GenMatrix.h:36
CHK_DIM_2
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: Vector.cc:42
init
void init()
Definition: testRandom.cc:27
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: Vector.cc:24
CLHEP::HepVector::HepVector
HepVector()
CLHEP::HepVector::operator-
HepVector operator-() const
Definition: Vector.cc:213
CLHEP::HepVector::operator-=
HepVector & operator-=(const HepMatrix &v2)
Definition: Vector.cc:430
CLHEP::HepMatrix::operator=
HepMatrix & operator=(const HepMatrix &)
Definition: Matrix.cc:417
CLHEP::detail::n
n
Definition: Ranlux64Engine.cc:85
CLHEP::HepMatrix
Definition: Matrix/CLHEP/Matrix/Matrix.h:209
CLHEP::HepVector::operator=
HepVector & operator=(const HepVector &hm2)
Definition: Vector.cc:469
CLHEP::HepMatrix::operator+=
HepMatrix & operator+=(const HepMatrix &)
Definition: Matrix.cc:391
CLHEP::HepVector::num_row
virtual int num_row() const
Definition: Vector.cc:117
CLHEP::HepVector::dot
friend double dot(const HepVector &v1, const HepVector &v2)
Definition: Vector.cc:543
f
void f(void g())
Definition: excDblThrow.cc:38
SIMPLE_TOP
#define SIMPLE_TOP(OPER)
Definition: Vector.cc:35
CLHEP
Definition: ClhepVersion.h:13
ij
ij
Definition: JamesRandomSeeding.txt:50
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::HepGenMatrix::mcIter
std::vector< double, Alloc< double, 25 > >::const_iterator mcIter
Definition: Matrix/CLHEP/Matrix/GenMatrix.h:78
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: Vector.cc:47
CLHEP::HepVector::operator()
const double & operator()(int row) const
CLHEP::HepVector::solve
friend HepVector solve(const HepMatrix &a, const HepVector &v)
Definition: Vector.cc:576
CLHEP::Hep3Vector
Definition: Geometry/CLHEP/Vector/ThreeVector.h:41
CLHEP::HepVector::operator+
friend HepVector operator+(const HepVector &v1, const HepVector &v2)
Definition: Vector.cc:256
CLHEP::HepVector::num_col
virtual int num_col() const
Definition: Vector.cc:119
j
long j
Definition: JamesRandomSeeding.txt:28
CLHEP::HepVector::T
HepMatrix T() const
Definition: Vector.cc:531
CLHEP::HepVector::HepMatrix
friend class HepMatrix
Definition: Matrix/CLHEP/Matrix/Vector.h:138
CLHEP::HepVector::sub
HepVector sub(int min_row, int max_row) const
Definition: Vector.cc:151
CLHEP::HepMatrix::operator-=
HepMatrix & operator-=(const HepMatrix &)
Definition: Matrix.cc:398
CLHEP::HepVector::operator*
friend HepVector operator*(const HepSymMatrix &hm1, const HepVector &hm2)
Definition: SymMatrix.cc:510
CLHEP::HepVector::~HepVector
virtual ~HepVector()
Definition: Vector.cc:92
CLHEP::HepVector::num_size
virtual int num_size() const
Definition: Vector.cc:118
CLHEP::HepVector::operator+=
HepVector & operator+=(const HepMatrix &v2)
Definition: Vector.cc:409
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::HepGenMatrix::mIter
std::vector< double, Alloc< double, 25 > >::iterator mIter
Definition: Matrix/CLHEP/Matrix/GenMatrix.h:77
CLHEP::HepVector::operator*=
HepVector & operator*=(double t)
Definition: Vector.cc:450
CLHEP::HepVector
Definition: Matrix/CLHEP/Matrix/Vector.h:39
CLHEP::HepGenMatrix::error
static void error(const char *s)
Definition: GenMatrix.cc:73
CLHEP::operator<<
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:86
CLHEP::HepVector::apply
HepVector apply(double(*f)(double, int)) const
Definition: Vector.cc:556
CLHEP::HepRandom
Definition: Matrix/CLHEP/Random/Random.h:50
CLHEP::dsum
HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2)
Definition: DiagMatrix.cc:164
CLHEP::HepMatrix::num_col
virtual int num_col() const
Definition: Matrix.cc:122
CLHEP::HepVector::operator/=
HepVector & operator/=(double t)
Definition: Vector.cc:444
SIMPLE_BOP
#define SIMPLE_BOP(OPER)
Definition: Vector.cc:29