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

DiagMatrix.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // ---------------------------------------------------------------------------
3 //
4 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 //
6 
7 #ifdef GNUPRAGMA
8 #pragma implementation
9 #endif
10 
11 #include <string.h>
12 #include <cmath>
13 
14 #include "CLHEP/Matrix/defs.h"
15 #include "CLHEP/Random/Random.h"
16 #include "CLHEP/Matrix/DiagMatrix.h"
17 #include "CLHEP/Matrix/Matrix.h"
18 #include "CLHEP/Matrix/SymMatrix.h"
19 #include "CLHEP/Matrix/Vector.h"
20 
21 #ifdef HEP_DEBUG_INLINE
22 #include "CLHEP/Matrix/DiagMatrix.icc"
23 #endif
24 
25 namespace CLHEP {
26 
27 // Simple operation for all elements
28 
29 #define SIMPLE_UOP(OPER) \
30  HepMatrix::mIter a=m.begin(); \
31  HepMatrix::mIter e=m.begin()+num_size(); \
32  for(;a<e; a++) (*a) OPER t;
33 
34 #define SIMPLE_BOP(OPER) \
35  HepMatrix::mIter a=m.begin(); \
36  HepMatrix::mcIter b=m2.m.begin(); \
37  HepMatrix::mIter e=m.begin()+num_size(); \
38  for(;a<e; a++, b++) (*a) OPER (*b);
39 
40 #define SIMPLE_TOP(OPER) \
41  HepMatrix::mcIter a=m1.m.begin(); \
42  HepMatrix::mcIter b=m2.m.begin(); \
43  HepMatrix::mIter t=mret.m.begin(); \
44  HepMatrix::mcIter e=m1.m.begin()+m1.nrow; \
45  for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b);
46 
47 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
48  if (r1!=r2 || c1!=c2) { \
49  HepGenMatrix::error("Range error in DiagMatrix function " #fun "(1)."); \
50  }
51 
52 #define CHK_DIM_1(c1,r2,fun) \
53  if (c1!=r2) { \
54  HepGenMatrix::error("Range error in DiagMatrix function " #fun "(2)."); \
55  }
56 
57 // static constant
58 
59 #if defined(__sun) || !defined(__GNUG__)
60 //
61 // Sun CC 4.0.1 has this bug.
62 //
63 double HepDiagMatrix::zero = 0;
64 #else
65 const double HepDiagMatrix::zero = 0;
66 #endif
67 
68 // Constructors. (Default constructors are inlined and in .icc file)
69 
71  : m(p), nrow(p)
72 {
73 }
74 
76  : m(p), nrow(p)
77 {
78  switch(init)
79  {
80  case 0:
81  m.assign(nrow,0);
82  break;
83 
84  case 1:
85  {
86  HepMatrix::mIter a=m.begin();
87  HepMatrix::mIter b=m.begin() + p;
88  for( ; a<b; a++) *a = 1.0;
89  break;
90  }
91  default:
92  error("DiagMatrix: initialization must be either 0 or 1.");
93  }
94 }
95 
97  : m(p), nrow(p)
98 {
99  HepMatrix::mIter a = m.begin();
100  HepMatrix::mIter b = m.begin() + num_size();
101  for(;a<b;a++) *a = r();
102 }
103 //
104 // Destructor
105 //
107 }
108 
110  : HepGenMatrix(m1), m(m1.nrow), nrow(m1.nrow)
111 {
112  m = m1.m;
113 }
114 
115 //
116 //
117 // Sub matrix
118 //
119 //
120 
121 HepDiagMatrix HepDiagMatrix::sub(int min_row, int max_row) const
122 #ifdef HEP_GNU_OPTIMIZED_RETURN
123 return mret(max_row-min_row+1);
124 {
125 #else
126 {
127  HepDiagMatrix mret(max_row-min_row+1);
128 #endif
129  if(max_row > num_row())
130  error("HepDiagMatrix::sub: Index out of range");
131  HepMatrix::mIter a = mret.m.begin();
132  HepMatrix::mcIter b = m.begin() + min_row - 1;
133  HepMatrix::mIter e = mret.m.begin() + mret.num_row();
134  for(;a<e;) *(a++) = *(b++);
135  return mret;
136 }
137 
138 HepDiagMatrix HepDiagMatrix::sub(int min_row, int max_row)
139 {
140  HepDiagMatrix mret(max_row-min_row+1);
141  if(max_row > num_row())
142  error("HepDiagMatrix::sub: Index out of range");
143  HepMatrix::mIter a = mret.m.begin();
144  HepMatrix::mIter b = m.begin() + min_row - 1;
145  HepMatrix::mIter e = mret.m.begin() + mret.num_row();
146  for(;a<e;) *(a++) = *(b++);
147  return mret;
148 }
149 
150 void HepDiagMatrix::sub(int row,const HepDiagMatrix &m1)
151 {
152  if(row <1 || row+m1.num_row()-1 > num_row() )
153  error("HepDiagMatrix::sub: Index out of range");
154  HepMatrix::mcIter a = m1.m.begin();
155  HepMatrix::mIter b = m.begin() + row - 1;
156  HepMatrix::mcIter e = m1.m.begin() + m1.num_row();
157  for(;a<e;) *(b++) = *(a++);
158 }
159 
160 //
161 // Direct sum of two matricies
162 //
163 
165  const HepDiagMatrix &m2)
166 #ifdef HEP_GNU_OPTIMIZED_RETURN
167  return mret(m1.num_row() + m2.num_row(), 0);
168 {
169 #else
170 {
171  HepDiagMatrix mret(m1.num_row() + m2.num_row(),
172  0);
173 #endif
174  mret.sub(1,m1);
175  mret.sub(m1.num_row()+1,m2);
176  return mret;
177 }
178 
180 #ifdef HEP_GNU_OPTIMIZED_RETURN
181  return m2(nrow);
182 {
183 #else
184 {
185  HepDiagMatrix m2(nrow);
186 #endif
187  HepMatrix::mcIter a=m.begin();
188  HepMatrix::mIter b=m2.m.begin();
189  HepMatrix::mcIter e=m.begin()+num_size();
190  for(;a<e; a++, b++) (*b) = -(*a);
191  return m2;
192 }
193 
194 
195 
197 #ifdef HEP_GNU_OPTIMIZED_RETURN
198  return mret(m1);
199 {
200 #else
201 {
202  HepMatrix mret(m1);
203 #endif
204  CHK_DIM_2(m1.num_row(),m2.num_row(),
205  m1.num_col(),m2.num_col(),+);
206  mret += m2;
207  return mret;
208 }
209 
211 #ifdef HEP_GNU_OPTIMIZED_RETURN
212  return mret(m2);
213 {
214 #else
215 {
216  HepMatrix mret(m2);
217 #endif
218  CHK_DIM_2(m1.num_row(),m2.num_row(),
219  m1.num_col(),m2.num_col(),+);
220  mret += m1;
221  return mret;
222 }
223 
225 #ifdef HEP_GNU_OPTIMIZED_RETURN
226  return mret(m1.nrow);
227 {
228 #else
229 {
230  HepDiagMatrix mret(m1.nrow);
231 #endif
232  CHK_DIM_1(m1.nrow,m2.nrow,+);
233  SIMPLE_TOP(+)
234  return mret;
235 }
236 
237 HepSymMatrix operator+(const HepDiagMatrix &m1,const HepSymMatrix &m2)
238 #ifdef HEP_GNU_OPTIMIZED_RETURN
239  return mret(m2);
240 {
241 #else
242 {
243  HepSymMatrix mret(m2);
244 #endif
245  CHK_DIM_1(m1.num_row(),m2.num_row(),+);
246  mret += m1;
247  return mret;
248 }
249 
251 #ifdef HEP_GNU_OPTIMIZED_RETURN
252  return mret(m2);
253 {
254 #else
255 {
256  HepSymMatrix mret(m2);
257 #endif
258  CHK_DIM_1(m1.num_row(),m2.num_row(),+);
259  mret += m1;
260  return mret;
261 }
262 
263 //
264 // operator -
265 //
266 
268 #ifdef HEP_GNU_OPTIMIZED_RETURN
269  return mret(m1);
270 {
271 #else
272 {
273  HepMatrix mret(m1);
274 #endif
275  CHK_DIM_2(m1.num_row(),m2.num_row(),
276  m1.num_col(),m2.num_col(),-);
277  mret -= m2;
278  return mret;
279 }
281 #ifdef HEP_GNU_OPTIMIZED_RETURN
282  return mret(m1);
283 {
284 #else
285 {
286  HepMatrix mret(m1);
287 #endif
288  CHK_DIM_2(m1.num_row(),m2.num_row(),
289  m1.num_col(),m2.num_col(),-);
290  mret -= m2;
291  return mret;
292 }
293 
295 #ifdef HEP_GNU_OPTIMIZED_RETURN
296  return mret(m1.nrow);
297 {
298 #else
299 {
300  HepDiagMatrix mret(m1.nrow);
301 #endif
302  CHK_DIM_1(m1.num_row(),m2.num_row(),-);
303  SIMPLE_TOP(-)
304  return mret;
305 }
306 HepSymMatrix operator-(const HepDiagMatrix &m1,const HepSymMatrix &m2)
307 #ifdef HEP_GNU_OPTIMIZED_RETURN
308  return mret(m1);
309 {
310 #else
311 {
312  HepSymMatrix mret(m1);
313 #endif
314  CHK_DIM_1(m1.num_row(),m2.num_row(),-);
315  mret -= m2;
316  return mret;
317 }
318 
320 #ifdef HEP_GNU_OPTIMIZED_RETURN
321  return mret(m1);
322 {
323 #else
324 {
325  HepSymMatrix mret(m1);
326 #endif
327  CHK_DIM_1(m1.num_row(),m2.num_row(),-);
328  mret -= m2;
329  return mret;
330 }
331 
332 /* -----------------------------------------------------------------------
333  This section contains support routines for matrix.h. This file contains
334  The two argument functions *,/. They call copy constructor and then /=,*=.
335  Also contains v_times_vT(const HepVector &v).
336  ----------------------------------------------------------------------- */
337 
339 const HepDiagMatrix &m1,double t)
340 #ifdef HEP_GNU_OPTIMIZED_RETURN
341  return mret(m1);
342 {
343 #else
344 {
345  HepDiagMatrix mret(m1);
346 #endif
347  mret /= t;
348  return mret;
349 }
350 
352 #ifdef HEP_GNU_OPTIMIZED_RETURN
353  return mret(m1);
354 {
355 #else
356 {
357  HepDiagMatrix mret(m1);
358 #endif
359  mret *= t;
360  return mret;
361 }
362 
364 #ifdef HEP_GNU_OPTIMIZED_RETURN
365  return mret(m1);
366 {
367 #else
368 {
369  HepDiagMatrix mret(m1);
370 #endif
371  mret *= t;
372  return mret;
373 }
374 
376 #ifdef HEP_GNU_OPTIMIZED_RETURN
377  return mret(m1.num_row(),m2.num_col());
378 {
379 #else
380  {
381  HepMatrix mret(m1.num_row(),m2.num_col());
382 #endif
383  CHK_DIM_1(m1.num_col(),m2.num_row(),*);
384  HepMatrix::mcIter mit1=m1.m.begin();
385  HepMatrix::mIter mir=mret.m.begin();
386  for(int irow=1;irow<=m1.num_row();irow++) {
387  HepMatrix::mcIter mcc = m2.m.begin();
388  for(int icol=1;icol<=m1.num_col();icol++) {
389  *(mir++) = *(mit1++) * (*(mcc++));
390  }
391  }
392  return mret;
393  }
394 
396 #ifdef HEP_GNU_OPTIMIZED_RETURN
397  return mret(m1.num_row(),m2.num_col());
398 {
399 #else
400 {
401  HepMatrix mret(m1.num_row(),m2.num_col());
402 #endif
403  CHK_DIM_1(m1.num_col(),m2.num_row(),*);
404  HepMatrix::mcIter mit1=m2.m.begin();
405  HepMatrix::mIter mir=mret.m.begin();
406  HepMatrix::mcIter mrr = m1.m.begin();
407  for(int irow=1;irow<=m2.num_row();irow++) {
408  for(int icol=1;icol<=m2.num_col();icol++) {
409  *(mir++) = *(mit1++) * (*mrr);
410  }
411  mrr++;
412  }
413  return mret;
414 }
415 
417 #ifdef HEP_GNU_OPTIMIZED_RETURN
418  return mret(m1.num_row());
419 {
420 #else
421 {
422  HepDiagMatrix mret(m1.num_row());
423 #endif
424  CHK_DIM_1(m1.num_col(),m2.num_row(),*);
425  HepMatrix::mIter a = mret.m.begin();
426  HepMatrix::mcIter b = m1.m.begin();
427  HepMatrix::mcIter c = m2.m.begin();
428  HepMatrix::mIter e = mret.m.begin() + m1.num_col();
429  for(;a<e;) *(a++) = *(b++) * (*(c++));
430  return mret;
431 }
432 
434 #ifdef HEP_GNU_OPTIMIZED_RETURN
435  return mret(m1.num_row());
436 {
437 #else
438 {
439  HepVector mret(m1.num_row());
440 #endif
441  CHK_DIM_1(m1.num_col(),m2.num_row(),*);
442  HepGenMatrix::mIter mir=mret.m.begin();
443  HepGenMatrix::mcIter mi1 = m1.m.begin(), mi2 = m2.m.begin();
444  for(int icol=1;icol<=m1.num_col();icol++) {
445  *(mir++) = *(mi1++) * *(mi2++);
446  }
447  return mret;
448 }
449 
450 /* -----------------------------------------------------------------------
451  This section contains the assignment and inplace operators =,+=,-=,*=,/=.
452  ----------------------------------------------------------------------- */
453 
455 {
456  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
457  int n = num_row();
458  mIter mrr = m.begin();
459  HepMatrix::mcIter mr = m2.m.begin();
460  for(int r=1;r<=n;r++) {
461  *mrr += *(mr++);
462  if(r<n) mrr += (n+1);
463  }
464  return (*this);
465 }
466 
468 {
469  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
470  HepMatrix::mIter a=m.begin();
471  HepMatrix::mcIter b=m2.m.begin();
472  for(int i=1;i<=num_row();i++) {
473  *a += *(b++);
474  if(i<num_row()) a += (i+1);
475  }
476  return (*this);
477 }
478 
480 {
481  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
482  SIMPLE_BOP(+=)
483  return (*this);
484 }
485 
487 {
488  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
489  int n = num_row();
490  mIter mrr = m.begin();
491  HepMatrix::mcIter mr = m2.m.begin();
492  for(int r=1;r<=n;r++) {
493  *mrr -= *(mr++);
494  if(r<n) mrr += (n+1);
495  }
496  return (*this);
497 }
498 
500 {
501  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
502  HepMatrix::mIter a=m.begin();
503  HepMatrix::mcIter b=m2.m.begin();
504  for(int i=1;i<=num_row();i++) {
505  *a -= *(b++);
506  if(i<num_row()) a += (i+1);
507  }
508  return (*this);
509 }
510 
512 {
513  CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
514  SIMPLE_BOP(-=)
515  return (*this);
516 }
517 
519 {
520  SIMPLE_UOP(/=)
521  return (*this);
522 }
523 
525 {
526  SIMPLE_UOP(*=)
527  return (*this);
528 }
529 
531 {
532  if(m1.nrow*m1.nrow != size_)
533  {
534  size_ = m1.nrow * m1.nrow;
535  m.resize(size_);
536  }
537  nrow = m1.nrow;
538  ncol = m1.nrow;
539  int n = nrow;
540  m.assign(size_,0);
541  mIter mrr = m.begin();
542  HepMatrix::mcIter mr = m1.m.begin();
543  for(int r=1;r<=n;r++) {
544  *mrr = *(mr++);
545  if(r<n) mrr += (n+1);
546  }
547  return (*this);
548 }
549 
551 {
552  if(m1.nrow != nrow)
553  {
554  nrow = m1.nrow;
555  m.resize(nrow);
556  }
557  m=m1.m;
558  return (*this);
559 }
560 
561 // Print the Matrix.
562 
563 std::ostream& operator<<(std::ostream &s, const HepDiagMatrix &q)
564 {
565  s << "\n";
566 /* Fixed format needs 3 extra characters for field, while scientific needs 7 */
567  int width;
568  if(s.flags() & std::ios::fixed)
569  width = s.precision()+3;
570  else
571  width = s.precision()+7;
572  for(int irow = 1; irow<= q.num_row(); irow++)
573  {
574  for(int icol = 1; icol <= q.num_col(); icol++)
575  {
576  s.width(width);
577  s << q(irow,icol) << " ";
578  }
579  s << std::endl;
580  }
581  return s;
582 }
583 
584 HepDiagMatrix HepDiagMatrix::
585 apply(double (*f)(double, int, int)) const
586 #ifdef HEP_GNU_OPTIMIZED_RETURN
587 return mret(num_row());
588 {
589 #else
590 {
591  HepDiagMatrix mret(num_row());
592 #endif
593  HepMatrix::mcIter a = m.begin();
594  HepMatrix::mIter b = mret.m.begin();
595  for(int ir=1;ir<=num_row();ir++) {
596  *(b++) = (*f)(*(a++), ir, ir);
597  }
598  return mret;
599 }
600 
602 {
603  if(m1.num_row()!=nrow)
604  {
605  nrow = m1.num_row();
606  m.resize(nrow);
607  }
608  HepMatrix::mcIter a = m1.m.begin();
609  HepMatrix::mIter b = m.begin();
610  for(int r=1;r<=nrow;r++) {
611  *(b++) = *a;
612  if(r<nrow) a += (nrow+1);
613  }
614 }
615 
617 {
618  if(m1.num_row()!=nrow)
619  {
620  nrow = m1.num_row();
621  m.resize(nrow);
622  }
623  HepMatrix::mcIter a = m1.m.begin();
624  HepMatrix::mIter b = m.begin();
625  for(int r=1;r<=nrow;r++) {
626  *(b++) = *a;
627  if(r<nrow) a += (r+1);
628  }
629 }
630 
632 #ifdef HEP_GNU_OPTIMIZED_RETURN
633  return mret(m1.num_row());
634 {
635 #else
636 {
637  HepSymMatrix mret(m1.num_row());
638 #endif
639  CHK_DIM_1(num_row(),m1.num_col(),"similarity");
640 // HepMatrix temp = m1*(*this);
641 // If m1*(*this) has correct dimensions, then so will the m1.T multiplication.
642 // So there is no need to check dimensions again.
643  HepMatrix::mIter mrc = mret.m.begin();
644  for(int r=1;r<=mret.num_row();r++) {
645  HepMatrix::mcIter mrr = m1.m.begin()+(r-1)*m1.num_col();
646  HepMatrix::mcIter mc = m1.m.begin();
647  for(int c=1;c<=r;c++) {
648  HepMatrix::mcIter mi = m.begin();
649  register double tmp = 0;
650  HepMatrix::mcIter mr = mrr;
651  for(int i=0;i<m1.num_col();i++)
652  tmp+=*(mr++) * *(mc++) * *(mi++);
653  *(mrc++) = tmp;
654  }
655  }
656  return mret;
657 }
658 
659 double HepDiagMatrix::similarity(const HepVector &m1) const
660 {
661  register double mret;
663  HepMatrix::mcIter mi = m.begin();
664  HepMatrix::mcIter mv = m1.m.begin();
665  mret = *(mv)* *(mv)* *(mi++);
666  mv++;
667  for(int i=2;i<=m1.num_row();i++) {
668  mret+=*(mv)* *(mv)* *(mi++);
669  mv++;
670  }
671  return mret;
672 }
673 
675 #ifdef HEP_GNU_OPTIMIZED_RETURN
676  return mret(m1.num_col());
677 {
678 #else
679 {
680  HepSymMatrix mret(m1.num_col());
681 #endif
682  CHK_DIM_1(num_col(),m1.num_row(),similarityT);
683 // Matrix temp = (*this)*m1;
684 // If m1*(*this) has correct dimensions, then so will the m1.T multiplication.
685 // So there is no need to check dimensions again.
686  for(int r=1;r<=mret.num_row();r++)
687  for(int c=1;c<=r;c++)
688  {
689  HepMatrix::mcIter mi = m.begin();
690  register double tmp = m1(1,r)*m1(1,c)* *(mi++);
691  for(int i=2;i<=m1.num_row();i++)
692  tmp+=m1(i,r)*m1(i,c)* *(mi++);
693  mret.fast(r,c) = tmp;
694  }
695  return mret;
696 }
697 
698 void HepDiagMatrix::invert(int &ierr) {
699  int n = num_row();
700  ierr = 1;
701  HepMatrix::mIter mm = m.begin();
702  int i;
703  for(i=0;i<n;i++) {
704  if(*(mm++)==0) return;
705  }
706  ierr = 0;
707  mm = m.begin();
708  for(i=0;i<n;i++) {
709  *mm = 1.0 / *mm;
710  mm++;
711  }
712 }
713 
715  double d = 1.0;
716  HepMatrix::mcIter end = m.begin() + nrow;
717  for (HepMatrix::mcIter p=m.begin(); p < end; p++)
718  d *= *p;
719  return d;
720 }
721 
722 double HepDiagMatrix::trace() const {
723  double d = 0.0;
724  HepMatrix::mcIter end = m.begin() + nrow;
725  for (HepMatrix::mcIter p=m.begin(); p < end; p++)
726  d += *p;
727  return d;
728 }
729 
730 } // namespace CLHEP