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
20namespace 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//
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
117inline int HepVector::num_row() const {return nrow;}
118inline int HepVector::num_size() const {return nrow;}
119inline int HepVector::num_col() const { return 1; }
120
121// operator()
122
123#ifdef MATRIX_BOUND_CHECK
124inline 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
129inline double & HepVector::operator()(int row, int)
130{
131#endif
132
133 return *(m.begin()+(row-1));
134}
135
136#ifdef MATRIX_BOUND_CHECK
137inline 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
142inline const double & HepVector::operator()(int row, int) const
143{
144#endif
145
146 return *(m.begin()+(row-1));
147}
148
149// Sub matrix
150
151HepVector HepVector::sub(int min_row, int max_row) const
152#ifdef HEP_GNU_OPTIMIZED_RETURN
153return 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
168HepVector 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
180void 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
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
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
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
273HepVector 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
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
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
318HepVector operator/(
319const 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
331HepVector 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
343HepVector 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
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
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(),*);
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
514std::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
533return 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
543double 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();
551 for(;a<e;) d += (*(a++)) * (*(b++));
552 return d;
553}
554
556apply(double (*f)(double, int)) const
557#ifdef HEP_GNU_OPTIMIZED_RETURN
558return 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
572void 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
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition DiagMatrix.cc:47
#define SIMPLE_BOP(OPER)
Definition DiagMatrix.cc:34
#define SIMPLE_UOP(OPER)
Definition DiagMatrix.cc:29
#define SIMPLE_TOP(OPER)
Definition DiagMatrix.cc:40
#define CHK_DIM_1(c1, r2, fun)
Definition DiagMatrix.cc:52
double z() const
double x() const
double y() const
std::vector< double, Alloc< double, 25 > >::iterator mIter
std::vector< double, Alloc< double, 25 > >::const_iterator mcIter
static void error(const char *s)
Definition GenMatrix.cc:73
HepMatrix & operator=(const HepMatrix &)
Definition Matrix.cc:417
virtual int num_col() const
Definition Matrix.cc:122
virtual int num_row() const
Definition Matrix.cc:120
HepMatrix & operator+=(const HepMatrix &)
Definition Matrix.cc:391
HepMatrix & operator-=(const HepMatrix &)
Definition Matrix.cc:398
HepVector & operator/=(double t)
Definition Vector.cc:444
HepVector sub(int min_row, int max_row) const
Definition Vector.cc:151
HepMatrix T() const
Definition Vector.cc:531
HepVector & operator-=(const HepMatrix &v2)
Definition Vector.cc:430
virtual ~HepVector()
Definition Vector.cc:92
virtual int num_size() const
Definition Vector.cc:118
virtual int num_col() const
Definition Vector.cc:119
HepVector apply(double(*f)(double, int)) const
Definition Vector.cc:556
HepVector operator-() const
Definition Vector.cc:213
virtual int num_row() const
Definition Vector.cc:117
friend double dot(const HepVector &v1, const HepVector &v2)
Definition Vector.cc:543
friend HepVector solve(const HepMatrix &a, const HepVector &v)
Definition Vector.cc:576
const double & operator()(int row) const
HepVector & operator=(const HepVector &hm2)
Definition Vector.cc:469
friend HepVector operator*(const HepSymMatrix &hm1, const HepVector &hm2)
Definition SymMatrix.cc:510
friend HepVector operator+(const HepVector &v1, const HepVector &v2)
Definition Vector.cc:256
HepVector & operator*=(double t)
Definition Vector.cc:450
HepVector & operator+=(const HepMatrix &v2)
Definition Vector.cc:409
void f(void g())
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition AxisAngle.cc:86
HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2)
void init()
Definition testRandom.cc:27