ROOT  6.06/09
Reference Guide
TVectorT.cxx
Go to the documentation of this file.
1 // @(#)root/matrix:$Id$
2 // Authors: Fons Rademakers, Eddy Offermann Nov 2003
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 //////////////////////////////////////////////////////////////////////////
13 // //
14 // TVectorT //
15 // //
16 // Template class of Vectors in the linear algebra package //
17 // //
18 // Unless otherwise specified, vector indices always start with 0, //
19 // spanning up to the specified limit-1. //
20 // //
21 // For (n) vectors where n <= kSizeMax (5 currently) storage space is //
22 // available on the stack, thus avoiding expensive allocation/ //
23 // deallocation of heap space . However, this introduces of course //
24 // kSizeMax overhead for each vector object . If this is an issue //
25 // recompile with a new appropriate value (>=0) for kSizeMax //
26 // //
27 // Another way to assign and store vector data is through Use //
28 // see for instance stressLinear.cxx file . //
29 // //
30 // Note that Constructors/assignments exists for all different matrix //
31 // views //
32 // //
33 // For usage examples see $ROOTSYS/test/stressLinear.cxx //
34 // //
35 //////////////////////////////////////////////////////////////////////////
36 
37 #include "TVectorT.h"
38 #include "TClass.h"
39 #include "TMath.h"
40 #include "TROOT.h"
41 #include "Varargs.h"
42 
43 templateClassImp(TVectorT)
44 
45 
46 ////////////////////////////////////////////////////////////////////////////////
47 /// Delete data pointer m, if it was assigned on the heap
48 
49 template<class Element>
50 void TVectorT<Element>::Delete_m(Int_t size,Element *&m)
51 {
52  if (m) {
53  if (size > kSizeMax)
54  delete [] m;
55  m = 0;
56  }
57 }
58 
59 ////////////////////////////////////////////////////////////////////////////////
60 /// Return data pointer . if requested size <= kSizeMax, assign pointer
61 /// to the stack space
62 
63 template<class Element>
64 Element* TVectorT<Element>::New_m(Int_t size)
65 {
66  if (size == 0) return 0;
67  else {
68  if ( size <= kSizeMax )
69  return fDataStack;
70  else {
71  Element *heap = new Element[size];
72  return heap;
73  }
74  }
75 }
76 
77 ////////////////////////////////////////////////////////////////////////////////
78 /// Add vector v to this vector
79 
80 template<class Element>
81 void TVectorT<Element>::Add(const TVectorT<Element> &v)
82 {
83  if (gMatrixCheck && !AreCompatible(*this,v)) {
84  Error("Add(TVectorT<Element> &)","vector's not compatible");
85  return;
86  }
87 
88  const Element *sp = v.GetMatrixArray();
89  Element *tp = this->GetMatrixArray();
90  const Element * const tp_last = tp+fNrows;
91  while (tp < tp_last)
92  *tp++ += *sp++;
93 }
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Set this vector to v1+v2
97 
98 template<class Element>
99 void TVectorT<Element>::Add(const TVectorT<Element> &v1,const TVectorT<Element> &v2)
100 {
101  if (gMatrixCheck) {
102  if ( !AreCompatible(*this,v1) && !AreCompatible(*this,v2)) {
103  Error("Add(TVectorT<Element> &)","vectors not compatible");
104  return;
105  }
106  }
107 
108  const Element *sv1 = v1.GetMatrixArray();
109  const Element *sv2 = v2.GetMatrixArray();
110  Element *tp = this->GetMatrixArray();
111  const Element * const tp_last = tp+fNrows;
112  while (tp < tp_last)
113  *tp++ = *sv1++ + *sv2++;
114 }
115 
116 ////////////////////////////////////////////////////////////////////////////////
117 /// Copy copySize doubles from *oldp to *newp . However take care of the
118 /// situation where both pointers are assigned to the same stack space
119 
120 template<class Element>
121 Int_t TVectorT<Element>::Memcpy_m(Element *newp,const Element *oldp,Int_t copySize,
122  Int_t newSize,Int_t oldSize)
123 {
124  if (copySize == 0 || oldp == newp) return 0;
125  else {
126  if ( newSize <= kSizeMax && oldSize <= kSizeMax ) {
127  // both pointers are inside fDataStack, be careful with copy direction !
128  if (newp > oldp) {
129  for (Int_t i = copySize-1; i >= 0; i--)
130  newp[i] = oldp[i];
131  } else {
132  for (Int_t i = 0; i < copySize; i++)
133  newp[i] = oldp[i];
134  }
135  }
136  else {
137  memcpy(newp,oldp,copySize*sizeof(Element));
138  }
139  }
140  return 0;
141 }
142 
143 ////////////////////////////////////////////////////////////////////////////////
144 /// Allocate new vector. Arguments are number of rows and row
145 /// lowerbound (0 default).
146 
147 template<class Element>
148 void TVectorT<Element>::Allocate(Int_t nrows,Int_t row_lwb,Int_t init)
149 {
150  fIsOwner = kTRUE;
151  fNrows = 0;
152  fRowLwb = 0;
153  fElements = 0;
154 
155  if (nrows < 0)
156  {
157  Error("Allocate","nrows=%d",nrows);
158  return;
159  }
160 
161  MakeValid();
162  fNrows = nrows;
163  fRowLwb = row_lwb;
164 
165  fElements = New_m(fNrows);
166  if (init)
167  memset(fElements,0,fNrows*sizeof(Element));
168 }
169 
170 ////////////////////////////////////////////////////////////////////////////////
171 /// Constructor n-vector
172 
173 template<class Element>
174 TVectorT<Element>::TVectorT(Int_t n)
175 {
176  Allocate(n,0,1);
177 }
178 
179 ////////////////////////////////////////////////////////////////////////////////
180 /// Constructor [lwb..upb]-vector
181 
182 template<class Element>
183 TVectorT<Element>::TVectorT(Int_t lwb,Int_t upb)
184 {
185  Allocate(upb-lwb+1,lwb,1);
186 }
187 
188 ////////////////////////////////////////////////////////////////////////////////
189 /// Constructor n-vector with data copied from array elements
190 
191 template<class Element>
192 TVectorT<Element>::TVectorT(Int_t n,const Element *elements)
193 {
194  Allocate(n,0);
195  SetElements(elements);
196 }
197 
198 ////////////////////////////////////////////////////////////////////////////////
199 /// Constructor [lwb..upb]-vector with data copied from array elements
200 
201 template<class Element>
202 TVectorT<Element>::TVectorT(Int_t lwb,Int_t upb,const Element *elements)
203 {
204  Allocate(upb-lwb+1,lwb);
205  SetElements(elements);
206 }
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 /// Copy constructor
210 
211 template<class Element>
212 TVectorT<Element>::TVectorT(const TVectorT &another) : TObject(another)
213 {
214  R__ASSERT(another.IsValid());
215  Allocate(another.GetUpb()-another.GetLwb()+1,another.GetLwb());
216  *this = another;
217 }
218 
219 ////////////////////////////////////////////////////////////////////////////////
220 /// Constructor : create vector from matrix row
221 
222 template<class Element>
223 TVectorT<Element>::TVectorT(const TMatrixTRow_const<Element> &mr) : TObject()
224 {
225  const TMatrixTBase<Element> *mt = mr.GetMatrix();
226  R__ASSERT(mt->IsValid());
227  Allocate(mt->GetColUpb()-mt->GetColLwb()+1,mt->GetColLwb());
228  *this = mr;
229 }
230 
231 ////////////////////////////////////////////////////////////////////////////////
232 /// Constructor : create vector from matrix column
233 
234 template<class Element>
235 TVectorT<Element>::TVectorT(const TMatrixTColumn_const<Element> &mc) : TObject()
236 {
237  const TMatrixTBase<Element> *mt = mc.GetMatrix();
238  R__ASSERT(mt->IsValid());
239  Allocate(mt->GetRowUpb()-mt->GetRowLwb()+1,mt->GetRowLwb());
240  *this = mc;
241 }
242 
243 ////////////////////////////////////////////////////////////////////////////////
244 /// Constructor : create vector from matrix diagonal
245 
246 template<class Element>
247 TVectorT<Element>::TVectorT(const TMatrixTDiag_const<Element> &md) : TObject()
248 {
249  const TMatrixTBase<Element> *mt = md.GetMatrix();
250  R__ASSERT(mt->IsValid());
251  Allocate(TMath::Min(mt->GetNrows(),mt->GetNcols()));
252  *this = md;
253 }
254 
255 ////////////////////////////////////////////////////////////////////////////////
256 /// Make a vector and assign initial values. Argument list should contain
257 /// Element values to assign to vector elements. The list must be
258 /// terminated by the string "END". Example:
259 /// TVectorT foo(1,3,0.0,1.0,1.5,"END");
260 
261 template<class Element>
262 TVectorT<Element>::TVectorT(Int_t lwb,Int_t upb,Element va_(iv1), ...)
263 {
264  if (upb < lwb) {
265  Error("TVectorT(Int_t, Int_t, ...)","upb(%d) < lwb(%d)",upb,lwb);
266  return;
267  }
268 
269  Allocate(upb-lwb+1,lwb);
270 
271  va_list args;
272  va_start(args,va_(iv1)); // Init 'args' to the beginning of
273  // the variable length list of args
274 
275  (*this)(lwb) = iv1;
276  for (Int_t i = lwb+1; i <= upb; i++)
277  (*this)(i) = (Element)va_arg(args,Double_t);
278 
279  if (strcmp((char *)va_arg(args,char *),"END"))
280  Error("TVectorT(Int_t, Int_t, ...)","argument list must be terminated by \"END\"");
281 
282  va_end(args);
283 }
284 
285 ////////////////////////////////////////////////////////////////////////////////
286 /// Resize the vector to [lwb:upb] .
287 /// New dynamic elemenst are created, the overlapping part of the old ones are
288 /// copied to the new structures, then the old elements are deleted.
289 
290 template<class Element>
291 TVectorT<Element> &TVectorT<Element>::ResizeTo(Int_t lwb,Int_t upb)
292 {
293  R__ASSERT(IsValid());
294  if (!fIsOwner) {
295  Error("ResizeTo(lwb,upb)","Not owner of data array,cannot resize");
296  return *this;
297  }
298 
299  const Int_t new_nrows = upb-lwb+1;
300 
301  if (fNrows > 0) {
302 
303  if (fNrows == new_nrows && fRowLwb == lwb)
304  return *this;
305  else if (new_nrows == 0) {
306  Clear();
307  return *this;
308  }
309 
310  Element *elements_old = GetMatrixArray();
311  const Int_t nrows_old = fNrows;
312  const Int_t rowLwb_old = fRowLwb;
313 
314  Allocate(new_nrows,lwb);
315  R__ASSERT(IsValid());
316  if (fNrows > kSizeMax || nrows_old > kSizeMax)
317  memset(GetMatrixArray(),0,fNrows*sizeof(Element));
318  else if (fNrows > nrows_old)
319  memset(GetMatrixArray()+nrows_old,0,(fNrows-nrows_old)*sizeof(Element));
320 
321  // Copy overlap
322  const Int_t rowLwb_copy = TMath::Max(fRowLwb,rowLwb_old);
323  const Int_t rowUpb_copy = TMath::Min(fRowLwb+fNrows-1,rowLwb_old+nrows_old-1);
324  const Int_t nrows_copy = rowUpb_copy-rowLwb_copy+1;
325 
326  const Int_t nelems_new = fNrows;
327  Element *elements_new = GetMatrixArray();
328  if (nrows_copy > 0) {
329  const Int_t rowOldOff = rowLwb_copy-rowLwb_old;
330  const Int_t rowNewOff = rowLwb_copy-fRowLwb;
331  Memcpy_m(elements_new+rowNewOff,elements_old+rowOldOff,nrows_copy,nelems_new,nrows_old);
332  }
333 
334  Delete_m(nrows_old,elements_old);
335  } else {
336  Allocate(upb-lwb+1,lwb,1);
337  }
338 
339  return *this;
340 }
341 
342 ////////////////////////////////////////////////////////////////////////////////
343 /// Use the array data to fill the vector lwb..upb]
344 
345 template<class Element>
346 TVectorT<Element> &TVectorT<Element>::Use(Int_t lwb,Int_t upb,Element *data)
347 {
348  if (upb < lwb) {
349  Error("Use","upb(%d) < lwb(%d)",upb,lwb);
350  return *this;
351  }
352 
353  Clear();
354  fNrows = upb-lwb+1;
355  fRowLwb = lwb;
356  fElements = data;
357  fIsOwner = kFALSE;
358 
359  return *this;
360 }
361 
362 ////////////////////////////////////////////////////////////////////////////////
363 /// Get subvector [row_lwb..row_upb]; The indexing range of the
364 /// returned vector depends on the argument option:
365 ///
366 /// option == "S" : return [0..row_upb-row_lwb+1] (default)
367 /// else : return [row_lwb..row_upb]
368 
369 template<class Element>
370 TVectorT<Element> &TVectorT<Element>::GetSub(Int_t row_lwb,Int_t row_upb,TVectorT<Element> &target,Option_t *option) const
371 {
372  if (gMatrixCheck) {
373  R__ASSERT(IsValid());
374  if (row_lwb < fRowLwb || row_lwb > fRowLwb+fNrows-1) {
375  Error("GetSub","row_lwb out of bounds");
376  return target;
377  }
378  if (row_upb < fRowLwb || row_upb > fRowLwb+fNrows-1) {
379  Error("GetSub","row_upb out of bounds");
380  return target;
381  }
382  if (row_upb < row_lwb) {
383  Error("GetSub","row_upb < row_lwb");
384  return target;
385  }
386  }
387 
388  TString opt(option);
389  opt.ToUpper();
390  const Int_t shift = (opt.Contains("S")) ? 1 : 0;
391 
392  Int_t row_lwb_sub;
393  Int_t row_upb_sub;
394  if (shift) {
395  row_lwb_sub = 0;
396  row_upb_sub = row_upb-row_lwb;
397  } else {
398  row_lwb_sub = row_lwb;
399  row_upb_sub = row_upb;
400  }
401 
402  target.ResizeTo(row_lwb_sub,row_upb_sub);
403  const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
404 
405  const Element *ap = this->GetMatrixArray()+(row_lwb-fRowLwb);
406  Element *bp = target.GetMatrixArray();
407 
408  for (Int_t irow = 0; irow < nrows_sub; irow++)
409  *bp++ = *ap++;
410 
411  return target;
412 }
413 
414 ////////////////////////////////////////////////////////////////////////////////
415 /// Insert vector source starting at [row_lwb], thereby overwriting the part
416 /// [row_lwb..row_lwb+nrows_source];
417 
418 template<class Element>
419 TVectorT<Element> &TVectorT<Element>::SetSub(Int_t row_lwb,const TVectorT<Element> &source)
420 {
421  if (gMatrixCheck) {
422  R__ASSERT(IsValid());
423  R__ASSERT(source.IsValid());
424 
425  if (row_lwb < fRowLwb && row_lwb > fRowLwb+fNrows-1) {
426  Error("SetSub","row_lwb outof bounds");
427  return *this;
428  }
429  if (row_lwb+source.GetNrows() > fRowLwb+fNrows) {
430  Error("SetSub","source vector too large");
431  return *this;
432  }
433  }
434 
435  const Int_t nRows_source = source.GetNrows();
436 
437  const Element *bp = source.GetMatrixArray();
438  Element *ap = this->GetMatrixArray()+(row_lwb-fRowLwb);
439 
440  for (Int_t irow = 0; irow < nRows_source; irow++)
441  *ap++ = *bp++;
442 
443  return *this;
444 }
445 
446 ////////////////////////////////////////////////////////////////////////////////
447 /// Set vector elements to zero
448 
449 template<class Element>
450 TVectorT<Element> &TVectorT<Element>::Zero()
451 {
452  R__ASSERT(IsValid());
453  memset(this->GetMatrixArray(),0,fNrows*sizeof(Element));
454  return *this;
455 }
456 
457 ////////////////////////////////////////////////////////////////////////////////
458 /// Take an absolute value of a vector, i.e. apply Abs() to each element.
459 
460 template<class Element>
461 TVectorT<Element> &TVectorT<Element>::Abs()
462 {
463  R__ASSERT(IsValid());
464 
465  Element *ep = this->GetMatrixArray();
466  const Element * const fp = ep+fNrows;
467  while (ep < fp) {
468  *ep = TMath::Abs(*ep);
469  ep++;
470  }
471 
472  return *this;
473 }
474 
475 ////////////////////////////////////////////////////////////////////////////////
476 /// Square each element of the vector.
477 
478 template<class Element>
479 TVectorT<Element> &TVectorT<Element>::Sqr()
480 {
481  R__ASSERT(IsValid());
482 
483  Element *ep = this->GetMatrixArray();
484  const Element * const fp = ep+fNrows;
485  while (ep < fp) {
486  *ep = (*ep) * (*ep);
487  ep++;
488  }
489 
490  return *this;
491 }
492 
493 ////////////////////////////////////////////////////////////////////////////////
494 /// Take square root of all elements.
495 
496 template<class Element>
497 TVectorT<Element> &TVectorT<Element>::Sqrt()
498 {
499  R__ASSERT(IsValid());
500 
501  Element *ep = this->GetMatrixArray();
502  const Element * const fp = ep+fNrows;
503  while (ep < fp) {
504  R__ASSERT(*ep >= 0);
505  if (*ep >= 0)
506  *ep = TMath::Sqrt(*ep);
507  else
508  Error("Sqrt()","v(%ld) = %g < 0",Long_t(ep-this->GetMatrixArray()),(float)*ep);
509  ep++;
510  }
511 
512  return *this;
513 }
514 
515 ////////////////////////////////////////////////////////////////////////////////
516 /// v[i] = 1/v[i]
517 
518 template<class Element>
519 TVectorT<Element> &TVectorT<Element>::Invert()
520 {
521  R__ASSERT(IsValid());
522 
523  Element *ep = this->GetMatrixArray();
524  const Element * const fp = ep+fNrows;
525  while (ep < fp) {
526  R__ASSERT(*ep != 0.0);
527  if (*ep != 0.0)
528  *ep = 1./ *ep;
529  else
530  Error("Invert()","v(%ld) = %g",Long_t(ep-this->GetMatrixArray()),(float)*ep);
531  ep++;
532  }
533 
534  return *this;
535 }
536 
537 ////////////////////////////////////////////////////////////////////////////////
538 /// Keep only element as selected through array select non-zero
539 
540 template<class Element>
541 TVectorT<Element> &TVectorT<Element>::SelectNonZeros(const TVectorT<Element> &select)
542 {
543  if (gMatrixCheck && !AreCompatible(*this,select)) {
544  Error("SelectNonZeros(const TVectorT<Element> &","vector's not compatible");
545  return *this;
546  }
547 
548  const Element *sp = select.GetMatrixArray();
549  Element *ep = this->GetMatrixArray();
550  const Element * const fp = ep+fNrows;
551  while (ep < fp) {
552  if (*sp == 0.0)
553  *ep = 0.0;
554  sp++; ep++;
555  }
556 
557  return *this;
558 }
559 
560 ////////////////////////////////////////////////////////////////////////////////
561 /// Compute the 1-norm of the vector SUM{ |v[i]| }.
562 
563 template<class Element>
564 Element TVectorT<Element>::Norm1() const
565 {
566  R__ASSERT(IsValid());
567 
568  Element norm = 0;
569  const Element *ep = this->GetMatrixArray();
570  const Element * const fp = ep+fNrows;
571  while (ep < fp)
572  norm += TMath::Abs(*ep++);
573 
574  return norm;
575 }
576 
577 ////////////////////////////////////////////////////////////////////////////////
578 /// Compute the square of the 2-norm SUM{ v[i]^2 }.
579 
580 template<class Element>
581 Element TVectorT<Element>::Norm2Sqr() const
582 {
583  R__ASSERT(IsValid());
584 
585  Element norm = 0;
586  const Element *ep = this->GetMatrixArray();
587  const Element * const fp = ep+fNrows;
588  while (ep < fp) {
589  norm += (*ep) * (*ep);
590  ep++;
591  }
592 
593  return norm;
594 }
595 
596 ////////////////////////////////////////////////////////////////////////////////
597 /// Compute the infinity-norm of the vector MAX{ |v[i]| }.
598 
599 template<class Element>
600 Element TVectorT<Element>::NormInf() const
601 {
602  R__ASSERT(IsValid());
603 
604  Element norm = 0;
605  const Element *ep = this->GetMatrixArray();
606  const Element * const fp = ep+fNrows;
607  while (ep < fp)
608  norm = TMath::Max(norm,TMath::Abs(*ep++));
609 
610  return norm;
611 }
612 
613 ////////////////////////////////////////////////////////////////////////////////
614 /// Compute the number of elements != 0.0
615 
616 template<class Element>
617 Int_t TVectorT<Element>::NonZeros() const
618 {
619  R__ASSERT(IsValid());
620 
621  Int_t nr_nonzeros = 0;
622  const Element *ep = this->GetMatrixArray();
623  const Element * const fp = ep+fNrows;
624  while (ep < fp)
625  if (*ep++) nr_nonzeros++;
626 
627  return nr_nonzeros;
628 }
629 
630 ////////////////////////////////////////////////////////////////////////////////
631 /// Compute sum of elements
632 
633 template<class Element>
634 Element TVectorT<Element>::Sum() const
635 {
636  R__ASSERT(IsValid());
637 
638  Element sum = 0.0;
639  const Element *ep = this->GetMatrixArray();
640  const Element * const fp = ep+fNrows;
641  while (ep < fp)
642  sum += *ep++;
643 
644  return sum;
645 }
646 
647 ////////////////////////////////////////////////////////////////////////////////
648 /// return minimum vector element value
649 
650 template<class Element>
651 Element TVectorT<Element>::Min() const
652 {
653  R__ASSERT(IsValid());
654 
655  const Int_t index = TMath::LocMin(fNrows,fElements);
656  return fElements[index];
657 }
658 
659 ////////////////////////////////////////////////////////////////////////////////
660 /// return maximum vector element value
661 
662 template<class Element>
663 Element TVectorT<Element>::Max() const
664 {
665  R__ASSERT(IsValid());
666 
667  const Int_t index = TMath::LocMax(fNrows,fElements);
668  return fElements[index];
669 }
670 
671 ////////////////////////////////////////////////////////////////////////////////
672 /// Notice that this assignment does NOT change the ownership :
673 /// if the storage space was adopted, source is copied to
674 /// this space .
675 
676 template<class Element>
677 TVectorT<Element> &TVectorT<Element>::operator=(const TVectorT<Element> &source)
678 {
679  if (gMatrixCheck && !AreCompatible(*this,source)) {
680  Error("operator=(const TVectorT<Element> &)","vectors not compatible");
681  return *this;
682  }
683 
684  if (this->GetMatrixArray() != source.GetMatrixArray()) {
685  TObject::operator=(source);
686  memcpy(fElements,source.GetMatrixArray(),fNrows*sizeof(Element));
687  }
688  return *this;
689 }
690 
691 ////////////////////////////////////////////////////////////////////////////////
692 /// Assign a matrix row to a vector.
693 
694 template<class Element>
695 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTRow_const<Element> &mr)
696 {
697  const TMatrixTBase<Element> *mt = mr.GetMatrix();
698 
699  if (gMatrixCheck) {
700  R__ASSERT(IsValid());
701  R__ASSERT(mt->IsValid());
702  if (mt->GetColLwb() != fRowLwb || mt->GetNcols() != fNrows) {
703  Error("operator=(const TMatrixTRow_const &)","vector and row not compatible");
704  return *this;
705  }
706  }
707 
708  const Int_t inc = mr.GetInc();
709  const Element *rp = mr.GetPtr(); // Row ptr
710  Element *ep = this->GetMatrixArray(); // Vector ptr
711  const Element * const fp = ep+fNrows;
712  while (ep < fp) {
713  *ep++ = *rp;
714  rp += inc;
715  }
716 
717  R__ASSERT(rp == mr.GetPtr()+mt->GetNcols());
718 
719  return *this;
720 }
721 
722 ////////////////////////////////////////////////////////////////////////////////
723 /// Assign a matrix column to a vector.
724 
725 template<class Element>
726 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTColumn_const<Element> &mc)
727 {
728  const TMatrixTBase<Element> *mt = mc.GetMatrix();
729 
730  if (gMatrixCheck) {
731  R__ASSERT(IsValid());
732  R__ASSERT(mt->IsValid());
733  if (mt->GetRowLwb() != fRowLwb || mt->GetNrows() != fNrows) {
734  Error("operator=(const TMatrixTColumn_const &)","vector and column not compatible");
735  return *this;
736  }
737  }
738 
739  const Int_t inc = mc.GetInc();
740  const Element *cp = mc.GetPtr(); // Column ptr
741  Element *ep = this->GetMatrixArray(); // Vector ptr
742  const Element * const fp = ep+fNrows;
743  while (ep < fp) {
744  *ep++ = *cp;
745  cp += inc;
746  }
747 
748  R__ASSERT(cp == mc.GetPtr()+mt->GetNoElements());
749 
750  return *this;
751 }
752 
753 ////////////////////////////////////////////////////////////////////////////////
754 /// Assign the matrix diagonal to a vector.
755 
756 template<class Element>
757 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTDiag_const<Element> &md)
758 {
759  const TMatrixTBase<Element> *mt = md.GetMatrix();
760 
761  if (gMatrixCheck) {
762  R__ASSERT(IsValid());
763  R__ASSERT(mt->IsValid());
764  if (md.GetNdiags() != fNrows) {
765  Error("operator=(const TMatrixTDiag_const &)","vector and matrix-diagonal not compatible");
766  return *this;
767  }
768  }
769 
770  const Int_t inc = md.GetInc();
771  const Element *dp = md.GetPtr(); // Diag ptr
772  Element *ep = this->GetMatrixArray(); // Vector ptr
773  const Element * const fp = ep+fNrows;
774  while (ep < fp) {
775  *ep++ = *dp;
776  dp += inc;
777  }
778 
779  R__ASSERT(dp < md.GetPtr()+mt->GetNoElements()+inc);
780 
781  return *this;
782 }
783 
784 ////////////////////////////////////////////////////////////////////////////////
785 /// Assign a sparse matrix row to a vector. The matrix row is implicitly transposed
786 /// to allow the assignment in the strict sense.
787 
788 template<class Element>
789 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTSparseRow_const<Element> &mr)
790 {
791  const TMatrixTBase<Element> *mt = mr.GetMatrix();
792 
793  if (gMatrixCheck) {
794  R__ASSERT(IsValid());
795  R__ASSERT(mt->IsValid());
796  if (mt->GetColLwb() != fRowLwb || mt->GetNcols() != fNrows) {
797  Error("operator=(const TMatrixTSparseRow_const &)","vector and row not compatible");
798  return *this;
799  }
800  }
801 
802  const Int_t nIndex = mr.GetNindex();
803  const Element * const prData = mr.GetDataPtr(); // Row Data ptr
804  const Int_t * const prCol = mr.GetColPtr(); // Col ptr
805  Element * const pvData = this->GetMatrixArray(); // Vector ptr
806 
807  memset(pvData,0,fNrows*sizeof(Element));
808  for (Int_t index = 0; index < nIndex; index++) {
809  const Int_t icol = prCol[index];
810  pvData[icol] = prData[index];
811  }
812 
813  return *this;
814 }
815 
816 ////////////////////////////////////////////////////////////////////////////////
817 /// Assign a sparse matrix diagonal to a vector.
818 
819 template<class Element>
820 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTSparseDiag_const<Element> &md)
821 {
822  const TMatrixTBase<Element> *mt = md.GetMatrix();
823 
824  if (gMatrixCheck) {
825  R__ASSERT(IsValid());
826  R__ASSERT(mt->IsValid());
827  if (md.GetNdiags() != fNrows) {
828  Error("operator=(const TMatrixTSparseDiag_const &)","vector and matrix-diagonal not compatible");
829  return *this;
830  }
831  }
832 
833  Element * const pvData = this->GetMatrixArray();
834  for (Int_t idiag = 0; idiag < fNrows; idiag++)
835  pvData[idiag] = md(idiag);
836 
837  return *this;
838 }
839 
840 ////////////////////////////////////////////////////////////////////////////////
841 /// Assign val to every element of the vector.
842 
843 template<class Element>
844 TVectorT<Element> &TVectorT<Element>::operator=(Element val)
845 {
846  R__ASSERT(IsValid());
847 
848  Element *ep = this->GetMatrixArray();
849  const Element * const fp = ep+fNrows;
850  while (ep < fp)
851  *ep++ = val;
852 
853  return *this;
854 }
855 
856 ////////////////////////////////////////////////////////////////////////////////
857 /// Add val to every element of the vector.
858 
859 template<class Element>
860 TVectorT<Element> &TVectorT<Element>::operator+=(Element val)
861 {
862  R__ASSERT(IsValid());
863 
864  Element *ep = this->GetMatrixArray();
865  const Element * const fp = ep+fNrows;
866  while (ep < fp)
867  *ep++ += val;
868 
869  return *this;
870 }
871 
872 ////////////////////////////////////////////////////////////////////////////////
873 /// Subtract val from every element of the vector.
874 
875 template<class Element>
876 TVectorT<Element> &TVectorT<Element>::operator-=(Element val)
877 {
878  R__ASSERT(IsValid());
879 
880  Element *ep = this->GetMatrixArray();
881  const Element * const fp = ep+fNrows;
882  while (ep < fp)
883  *ep++ -= val;
884 
885  return *this;
886 }
887 
888 ////////////////////////////////////////////////////////////////////////////////
889 /// Multiply every element of the vector with val.
890 
891 template<class Element>
892 TVectorT<Element> &TVectorT<Element>::operator*=(Element val)
893 {
894  R__ASSERT(IsValid());
895 
896  Element *ep = this->GetMatrixArray();
897  const Element * const fp = ep+fNrows;
898  while (ep < fp)
899  *ep++ *= val;
900 
901  return *this;
902 }
903 
904 ////////////////////////////////////////////////////////////////////////////////
905 /// Add vector source
906 
907 template<class Element>
908 TVectorT<Element> &TVectorT<Element>::operator+=(const TVectorT<Element> &source)
909 {
910  if (gMatrixCheck && !AreCompatible(*this,source)) {
911  Error("operator+=(const TVectorT<Element> &)","vector's not compatible");
912  return *this;
913  }
914 
915  const Element *sp = source.GetMatrixArray();
916  Element *tp = this->GetMatrixArray();
917  const Element * const tp_last = tp+fNrows;
918  while (tp < tp_last)
919  *tp++ += *sp++;
920 
921  return *this;
922 }
923 
924 ////////////////////////////////////////////////////////////////////////////////
925 /// Subtract vector source
926 
927 template<class Element>
928 TVectorT<Element> &TVectorT<Element>::operator-=(const TVectorT<Element> &source)
929 {
930  if (gMatrixCheck && !AreCompatible(*this,source)) {
931  Error("operator-=(const TVectorT<Element> &)","vector's not compatible");
932  return *this;
933  }
934 
935  const Element *sp = source.GetMatrixArray();
936  Element *tp = this->GetMatrixArray();
937  const Element * const tp_last = tp+fNrows;
938  while (tp < tp_last)
939  *tp++ -= *sp++;
940 
941  return *this;
942 }
943 
944 ////////////////////////////////////////////////////////////////////////////////
945 /// "Inplace" multiplication target = A*target. A needn't be a square one
946 /// If target has to be resized, it should own the storage: fIsOwner = kTRUE
947 
948 template<class Element>
949 TVectorT<Element> &TVectorT<Element>::operator*=(const TMatrixT<Element> &a)
950 {
951  if (gMatrixCheck) {
952  R__ASSERT(IsValid());
953  R__ASSERT(a.IsValid());
954  if (a.GetNcols() != fNrows || a.GetColLwb() != fRowLwb) {
955  Error("operator*=(const TMatrixT &)","vector and matrix incompatible");
956  return *this;
957  }
958  }
959 
960  const Bool_t doResize = (fNrows != a.GetNrows() || fRowLwb != a.GetRowLwb());
961  if (doResize && !fIsOwner) {
962  Error("operator*=(const TMatrixT &)","vector has to be resized but not owner");
963  return *this;
964  }
965 
966  Element work[kWorkMax];
967  Bool_t isAllocated = kFALSE;
968  Element *elements_old = work;
969  const Int_t nrows_old = fNrows;
970  if (nrows_old > kWorkMax) {
971  isAllocated = kTRUE;
972  elements_old = new Element[nrows_old];
973  }
974  memcpy(elements_old,fElements,nrows_old*sizeof(Element));
975 
976  if (doResize) {
977  const Int_t rowlwb_new = a.GetRowLwb();
978  const Int_t nrows_new = a.GetNrows();
979  ResizeTo(rowlwb_new,rowlwb_new+nrows_new-1);
980  }
981  memset(fElements,0,fNrows*sizeof(Element));
982 
983  const Element *mp = a.GetMatrixArray(); // Matrix row ptr
984  Element *tp = this->GetMatrixArray(); // Target vector ptr
985 #ifdef CBLAS
986  if (typeid(Element) == typeid(Double_t))
987  cblas_dgemv(CblasRowMajor,CblasNoTrans,a.GetNrows(),a.GetNcols(),1.0,mp,
988  a.GetNcols(),elements_old,1,0.0,tp,1);
989  else if (typeid(Element) != typeid(Float_t))
990  cblas_sgemv(CblasRowMajor,CblasNoTrans,a.GetNrows(),a.GetNcols(),1.0,mp,
991  a.GetNcols(),elements_old,1,0.0,tp,1);
992  else
993  Error("operator*=","type %s not implemented in BLAS library",typeid(Element));
994 #else
995  const Element * const tp_last = tp+fNrows;
996  while (tp < tp_last) {
997  Element sum = 0;
998  for (const Element *sp = elements_old; sp < elements_old+nrows_old; )
999  sum += *sp++ * *mp++;
1000  *tp++ = sum;
1001  }
1002  R__ASSERT(mp == a.GetMatrixArray()+a.GetNoElements());
1003 #endif
1004 
1005  if (isAllocated)
1006  delete [] elements_old;
1007 
1008  return *this;
1009 }
1010 
1011 ////////////////////////////////////////////////////////////////////////////////
1012 /// "Inplace" multiplication target = A*target. A needn't be a square one
1013 /// If target has to be resized, it should own the storage: fIsOwner = kTRUE
1014 
1015 template<class Element>
1016 TVectorT<Element> &TVectorT<Element>::operator*=(const TMatrixTSparse<Element> &a)
1017 {
1018  if (gMatrixCheck) {
1019  R__ASSERT(IsValid());
1020  R__ASSERT(a.IsValid());
1021  if (a.GetNcols() != fNrows || a.GetColLwb() != fRowLwb) {
1022  Error("operator*=(const TMatrixTSparse &)","vector and matrix incompatible");
1023  return *this;
1024  }
1025  }
1026 
1027  const Bool_t doResize = (fNrows != a.GetNrows() || fRowLwb != a.GetRowLwb());
1028  if (doResize && !fIsOwner) {
1029  Error("operator*=(const TMatrixTSparse &)","vector has to be resized but not owner");
1030  return *this;
1031  }
1032 
1033  Element work[kWorkMax];
1034  Bool_t isAllocated = kFALSE;
1035  Element *elements_old = work;
1036  const Int_t nrows_old = fNrows;
1037  if (nrows_old > kWorkMax) {
1038  isAllocated = kTRUE;
1039  elements_old = new Element[nrows_old];
1040  }
1041  memcpy(elements_old,fElements,nrows_old*sizeof(Element));
1042 
1043  if (doResize) {
1044  const Int_t rowlwb_new = a.GetRowLwb();
1045  const Int_t nrows_new = a.GetNrows();
1046  ResizeTo(rowlwb_new,rowlwb_new+nrows_new-1);
1047  }
1048  memset(fElements,0,fNrows*sizeof(Element));
1049 
1050  const Int_t * const pRowIndex = a.GetRowIndexArray();
1051  const Int_t * const pColIndex = a.GetColIndexArray();
1052  const Element * const mp = a.GetMatrixArray(); // Matrix row ptr
1053 
1054  const Element * const sp = elements_old;
1055  Element * tp = this->GetMatrixArray(); // Target vector ptr
1056 
1057  for (Int_t irow = 0; irow < fNrows; irow++) {
1058  const Int_t sIndex = pRowIndex[irow];
1059  const Int_t eIndex = pRowIndex[irow+1];
1060  Element sum = 0.0;
1061  for (Int_t index = sIndex; index < eIndex; index++) {
1062  const Int_t icol = pColIndex[index];
1063  sum += mp[index]*sp[icol];
1064  }
1065  tp[irow] = sum;
1066  }
1067 
1068  if (isAllocated)
1069  delete [] elements_old;
1070 
1071  return *this;
1072 }
1073 
1074 ////////////////////////////////////////////////////////////////////////////////
1075 /// "Inplace" multiplication target = A*target. A is symmetric .
1076 /// vector size will not change
1077 
1078 template<class Element>
1079 TVectorT<Element> &TVectorT<Element>::operator*=(const TMatrixTSym<Element> &a)
1080 {
1081  if (gMatrixCheck) {
1082  R__ASSERT(IsValid());
1083  R__ASSERT(a.IsValid());
1084  if (a.GetNcols() != fNrows || a.GetColLwb() != fRowLwb) {
1085  Error("operator*=(const TMatrixTSym &)","vector and matrix incompatible");
1086  return *this;
1087  }
1088  }
1089 
1090  Element work[kWorkMax];
1091  Bool_t isAllocated = kFALSE;
1092  Element *elements_old = work;
1093  const Int_t nrows_old = fNrows;
1094  if (nrows_old > kWorkMax) {
1095  isAllocated = kTRUE;
1096  elements_old = new Element[nrows_old];
1097  }
1098  memcpy(elements_old,fElements,nrows_old*sizeof(Element));
1099  memset(fElements,0,fNrows*sizeof(Element));
1100 
1101  const Element *mp = a.GetMatrixArray(); // Matrix row ptr
1102  Element *tp = this->GetMatrixArray(); // Target vector ptr
1103 #ifdef CBLAS
1104  if (typeid(Element) == typeid(Double_t))
1105  cblas_dsymv(CblasRowMajor,CblasUpper,fNrows,1.0,mp,
1106  fNrows,elements_old,1,0.0,tp,1);
1107  else if (typeid(Element) != typeid(Float_t))
1108  cblas_ssymv(CblasRowMajor,CblasUpper,fNrows,1.0,mp,
1109  fNrows,elements_old,1,0.0,tp,1);
1110  else
1111  Error("operator*=","type %s not implemented in BLAS library",typeid(Element));
1112 #else
1113  const Element * const tp_last = tp+fNrows;
1114  while (tp < tp_last) {
1115  Element sum = 0;
1116  for (const Element *sp = elements_old; sp < elements_old+nrows_old; )
1117  sum += *sp++ * *mp++;
1118  *tp++ = sum;
1119  }
1120  R__ASSERT(mp == a.GetMatrixArray()+a.GetNoElements());
1121 #endif
1122 
1123  if (isAllocated)
1124  delete [] elements_old;
1125 
1126  return *this;
1127 }
1128 
1129 ////////////////////////////////////////////////////////////////////////////////
1130 /// Are all vector elements equal to val?
1131 
1132 template<class Element>
1133 Bool_t TVectorT<Element>::operator==(Element val) const
1134 {
1135  R__ASSERT(IsValid());
1136 
1137  const Element *ep = this->GetMatrixArray();
1138  const Element * const fp = ep+fNrows;
1139  while (ep < fp)
1140  if (!(*ep++ == val))
1141  return kFALSE;
1142 
1143  return kTRUE;
1144 }
1145 
1146 ////////////////////////////////////////////////////////////////////////////////
1147 /// Are all vector elements not equal to val?
1148 
1149 template<class Element>
1150 Bool_t TVectorT<Element>::operator!=(Element val) const
1151 {
1152  R__ASSERT(IsValid());
1153 
1154  const Element *ep = this->GetMatrixArray();
1155  const Element * const fp = ep+fNrows;
1156  while (ep < fp)
1157  if (!(*ep++ != val))
1158  return kFALSE;
1159 
1160  return kTRUE;
1161 }
1162 
1163 ////////////////////////////////////////////////////////////////////////////////
1164 /// Are all vector elements < val?
1165 
1166 template<class Element>
1167 Bool_t TVectorT<Element>::operator<(Element val) const
1168 {
1169  R__ASSERT(IsValid());
1170 
1171  const Element *ep = this->GetMatrixArray();
1172  const Element * const fp = ep+fNrows;
1173  while (ep < fp)
1174  if (!(*ep++ < val))
1175  return kFALSE;
1176 
1177  return kTRUE;
1178 }
1179 
1180 ////////////////////////////////////////////////////////////////////////////////
1181 /// Are all vector elements <= val?
1182 
1183 template<class Element>
1184 Bool_t TVectorT<Element>::operator<=(Element val) const
1185 {
1186  R__ASSERT(IsValid());
1187 
1188  const Element *ep = this->GetMatrixArray();
1189  const Element * const fp = ep+fNrows;
1190  while (ep < fp)
1191  if (!(*ep++ <= val))
1192  return kFALSE;
1193 
1194  return kTRUE;
1195 }
1196 
1197 ////////////////////////////////////////////////////////////////////////////////
1198 /// Are all vector elements > val?
1199 
1200 template<class Element>
1201 Bool_t TVectorT<Element>::operator>(Element val) const
1202 {
1203  R__ASSERT(IsValid());
1204 
1205  const Element *ep = this->GetMatrixArray();
1206  const Element * const fp = ep+fNrows;
1207  while (ep < fp)
1208  if (!(*ep++ > val))
1209  return kFALSE;
1210 
1211  return kTRUE;
1212 }
1213 
1214 ////////////////////////////////////////////////////////////////////////////////
1215 /// Are all vector elements >= val?
1216 
1217 template<class Element>
1218 Bool_t TVectorT<Element>::operator>=(Element val) const
1219 {
1220  R__ASSERT(IsValid());
1221 
1222  const Element *ep = this->GetMatrixArray();
1223  const Element * const fp = ep+fNrows;
1224  while (ep < fp)
1225  if (!(*ep++ >= val))
1226  return kFALSE;
1227 
1228  return kTRUE;
1229 }
1230 
1231 ////////////////////////////////////////////////////////////////////////////////
1232 /// Check if vector elements as selected through array select are non-zero
1233 
1234 template<class Element>
1235 Bool_t TVectorT<Element>::MatchesNonZeroPattern(const TVectorT<Element> &select)
1236 {
1237  if (gMatrixCheck && !AreCompatible(*this,select)) {
1238  Error("MatchesNonZeroPattern(const TVectorT&)","vector's not compatible");
1239  return kFALSE;
1240  }
1241 
1242  const Element *sp = select.GetMatrixArray();
1243  const Element *ep = this->GetMatrixArray();
1244  const Element * const fp = ep+fNrows;
1245  while (ep < fp) {
1246  if (*sp == 0.0 && *ep != 0.0)
1247  return kFALSE;
1248  sp++; ep++;
1249  }
1250 
1251  return kTRUE;
1252 }
1253 
1254 ////////////////////////////////////////////////////////////////////////////////
1255 /// Check if vector elements as selected through array select are all positive
1256 
1257 template<class Element>
1258 Bool_t TVectorT<Element>::SomePositive(const TVectorT<Element> &select)
1259 {
1260  if (gMatrixCheck && !AreCompatible(*this,select)) {
1261  Error("SomePositive(const TVectorT&)","vector's not compatible");
1262  return kFALSE;
1263  }
1264 
1265  const Element *sp = select.GetMatrixArray();
1266  const Element *ep = this->GetMatrixArray();
1267  const Element * const fp = ep+fNrows;
1268  while (ep < fp) {
1269  if (*sp != 0.0 && *ep <= 0.0)
1270  return kFALSE;
1271  sp++; ep++;
1272  }
1273 
1274  return kTRUE;
1275 }
1276 
1277 ////////////////////////////////////////////////////////////////////////////////
1278 /// Add to vector elements as selected through array select the value val
1279 
1280 template<class Element>
1281 void TVectorT<Element>::AddSomeConstant(Element val,const TVectorT<Element> &select)
1282 {
1283  if (gMatrixCheck && !AreCompatible(*this,select))
1284  Error("AddSomeConstant(Element,const TVectorT&)(const TVectorT&)","vector's not compatible");
1285 
1286  const Element *sp = select.GetMatrixArray();
1287  Element *ep = this->GetMatrixArray();
1288  const Element * const fp = ep+fNrows;
1289  while (ep < fp) {
1290  if (*sp)
1291  *ep += val;
1292  sp++; ep++;
1293  }
1294 }
1295 
1296 extern Double_t Drand(Double_t &ix);
1297 
1298 ////////////////////////////////////////////////////////////////////////////////
1299 /// randomize vector elements value
1300 
1301 template<class Element>
1302 void TVectorT<Element>::Randomize(Element alpha,Element beta,Double_t &seed)
1303 {
1304  R__ASSERT(IsValid());
1305 
1306  const Element scale = beta-alpha;
1307  const Element shift = alpha/scale;
1308 
1309  Element * ep = GetMatrixArray();
1310  const Element * const fp = ep+fNrows;
1311  while (ep < fp)
1312  *ep++ = scale*(Drand(seed)+shift);
1313 }
1314 
1315 ////////////////////////////////////////////////////////////////////////////////
1316 /// Apply action to each element of the vector.
1317 
1318 template<class Element>
1319 TVectorT<Element> &TVectorT<Element>::Apply(const TElementActionT<Element> &action)
1320 {
1321  R__ASSERT(IsValid());
1322  for (Element *ep = fElements; ep < fElements+fNrows; ep++)
1323  action.Operation(*ep);
1324  return *this;
1325 }
1326 
1327 ////////////////////////////////////////////////////////////////////////////////
1328 /// Apply action to each element of the vector. In action the location
1329 /// of the current element is known.
1330 
1331 template<class Element>
1332 TVectorT<Element> &TVectorT<Element>::Apply(const TElementPosActionT<Element> &action)
1333 {
1334  R__ASSERT(IsValid());
1335 
1336  Element *ep = fElements;
1337  for (action.fI = fRowLwb; action.fI < fRowLwb+fNrows; action.fI++)
1338  action.Operation(*ep++);
1339 
1340  R__ASSERT(ep == fElements+fNrows);
1341 
1342  return *this;
1343 }
1344 
1345 ////////////////////////////////////////////////////////////////////////////////
1346 /// Draw this vector
1347 /// The histogram is named "TVectorT" by default and no title
1348 
1349 template<class Element>
1350 void TVectorT<Element>::Draw(Option_t *option)
1351 {
1352  gROOT->ProcessLine(Form("THistPainter::PaintSpecialObjects((TObject*)0x%lx,\"%s\");",
1353  (ULong_t)this, option));
1354 }
1355 
1356 ////////////////////////////////////////////////////////////////////////////////
1357 /// Print the vector as a list of elements.
1358 
1359 template<class Element>
1360 void TVectorT<Element>::Print(Option_t *flag) const
1361 {
1362  if (!IsValid()) {
1363  Error("Print","Vector is invalid");
1364  return;
1365  }
1366 
1367  printf("\nVector (%d) %s is as follows",fNrows,flag);
1368 
1369  printf("\n\n | %6d |", 1);
1370  printf("\n%s\n", "------------------");
1371  for (Int_t i = 0; i < fNrows; i++) {
1372  printf("%4d |",i+fRowLwb);
1373  //printf("%11.4g \n",(*this)(i+fRowLwb));
1374  printf("%g \n",(*this)(i+fRowLwb));
1375  }
1376  printf("\n");
1377 }
1378 
1379 ////////////////////////////////////////////////////////////////////////////////
1380 /// Check to see if two vectors are identical.
1381 
1382 template<class Element>
1383 Bool_t operator==(const TVectorT<Element> &v1,const TVectorT<Element> &v2)
1384 {
1385  if (!AreCompatible(v1,v2)) return kFALSE;
1386  return (memcmp(v1.GetMatrixArray(),v2.GetMatrixArray(),v1.GetNrows()*sizeof(Element)) == 0);
1387 }
1388 
1389 ////////////////////////////////////////////////////////////////////////////////
1390 /// Compute the scalar product.
1391 
1392 template<class Element>
1393 Element operator*(const TVectorT<Element> &v1,const TVectorT<Element> &v2)
1394 {
1395  if (gMatrixCheck) {
1396  if (!AreCompatible(v1,v2)) {
1397  Error("operator*(const TVectorT<Element> &,const TVectorT<Element> &)","vector's are incompatible");
1398  return 0.0;
1399  }
1400  }
1401 
1402  return Dot(v1,v2);
1403 }
1404 
1405 ////////////////////////////////////////////////////////////////////////////////
1406 /// Return source1+source2
1407 
1408 template<class Element>
1409 TVectorT<Element> operator+(const TVectorT<Element> &source1,const TVectorT<Element> &source2)
1410 {
1411  TVectorT<Element> target = source1;
1412  target += source2;
1413  return target;
1414 }
1415 
1416 ////////////////////////////////////////////////////////////////////////////////
1417 /// Return source1-source2
1418 
1419 template<class Element>
1420 TVectorT<Element> operator-(const TVectorT<Element> &source1,const TVectorT<Element> &source2)
1421 {
1422  TVectorT<Element> target = source1;
1423  target -= source2;
1424  return target;
1425 }
1426 
1427 ////////////////////////////////////////////////////////////////////////////////
1428 /// return A * source
1429 
1430 template<class Element>
1431 TVectorT<Element> operator*(const TMatrixT<Element> &a,const TVectorT<Element> &source)
1432 {
1433  R__ASSERT(a.IsValid());
1434  TVectorT<Element> target(a.GetRowLwb(),a.GetRowUpb());
1435  return Add(target,Element(1.0),a,source);
1436 }
1437 
1438 ////////////////////////////////////////////////////////////////////////////////
1439 /// return A * source
1440 
1441 template<class Element>
1442 TVectorT<Element> operator*(const TMatrixTSym<Element> &a,const TVectorT<Element> &source)
1443 {
1444  R__ASSERT(a.IsValid());
1445  TVectorT<Element> target(a.GetRowLwb(),a.GetRowUpb());
1446  return Add(target,Element(1.0),a,source);
1447 }
1448 
1449 ////////////////////////////////////////////////////////////////////////////////
1450 /// return A * source
1451 
1452 template<class Element>
1453 TVectorT<Element> operator*(const TMatrixTSparse<Element> &a,const TVectorT<Element> &source)
1454 {
1455  R__ASSERT(a.IsValid());
1456  TVectorT<Element> target(a.GetRowLwb(),a.GetRowUpb());
1457  return Add(target,Element(1.0),a,source);
1458 }
1459 
1460 ////////////////////////////////////////////////////////////////////////////////
1461 /// return val * source
1462 
1463 template<class Element>
1464 TVectorT<Element> operator*(Element val,const TVectorT<Element> &source)
1465 {
1466  TVectorT<Element> target = source;
1467  target *= val;
1468  return target;
1469 }
1470 
1471 ////////////////////////////////////////////////////////////////////////////////
1472 /// return inner-produvt v1 . v2
1473 
1474 template<class Element>
1475 Element Dot(const TVectorT<Element> &v1,const TVectorT<Element> &v2)
1476 {
1477  const Element *v1p = v1.GetMatrixArray();
1478  const Element *v2p = v2.GetMatrixArray();
1479  Element sum = 0.0;
1480  const Element * const fv1p = v1p+v1.GetNrows();
1481  while (v1p < fv1p)
1482  sum += *v1p++ * *v2p++;
1483 
1484  return sum;
1485 }
1486 
1487 ////////////////////////////////////////////////////////////////////////////////
1488 /// Return the matrix M = v1 * v2'
1489 
1490 template <class Element1, class Element2>
1491 TMatrixT<Element1>
1492 OuterProduct(const TVectorT<Element1> &v1,const TVectorT<Element2> &v2)
1493 {
1494  // TMatrixD::GetSub does:
1495  // TMatrixT tmp;
1496  // Doesn't compile here, because we are outside the class?
1497  // So we'll be explicit:
1498  TMatrixT<Element1> target;
1499 
1500  return OuterProduct(target,v1,v2);
1501 }
1502 
1503 ////////////////////////////////////////////////////////////////////////////////
1504 /// Return the matrix M = v1 * v2'
1505 
1506 template <class Element1,class Element2,class Element3>
1507 TMatrixT<Element1>
1508 &OuterProduct(TMatrixT<Element1> &target,const TVectorT<Element2> &v1,const TVectorT<Element3> &v2)
1509 {
1510  target.ResizeTo(v1.GetLwb(), v1.GetUpb(), v2.GetLwb(), v2.GetUpb());
1511 
1512  Element1 * mp = target.GetMatrixArray();
1513  const Element1 * const m_last = mp + target.GetNoElements();
1514 
1515  const Element2 * v1p = v1.GetMatrixArray();
1516  const Element2 * const v1_last = v1p + v1.GetNrows();
1517 
1518  const Element3 * const v20 = v2.GetMatrixArray();
1519  const Element3 * v2p = v20;
1520  const Element3 * const v2_last = v2p + v2.GetNrows();
1521 
1522  while (v1p < v1_last) {
1523  v2p = v20;
1524  while (v2p < v2_last) {
1525  *mp++ = *v1p * *v2p++ ;
1526  }
1527  v1p++;
1528  }
1529 
1530  R__ASSERT(v1p == v1_last && mp == m_last && v2p == v2_last);
1531 
1532  return target;
1533 }
1534 
1535 ////////////////////////////////////////////////////////////////////////////////
1536 /// Perform v1 * M * v2, a scalar result
1537 
1538 template <class Element1, class Element2, class Element3>
1539 Element1 Mult(const TVectorT<Element1> &v1,const TMatrixT<Element2> &m,
1540  const TVectorT<Element3> &v2)
1541 {
1542  if (gMatrixCheck) {
1543  if (!AreCompatible(v1, m)) {
1544  ::Error("Mult", "Vector v1 and matrix m incompatible");
1545  return 0;
1546  }
1547  if (!AreCompatible(m, v2)) {
1548  ::Error("Mult", "Matrix m and vector v2 incompatible");
1549  return 0;
1550  }
1551  }
1552 
1553  const Element1 * v1p = v1.GetMatrixArray(); // first of v1
1554  const Element1 * const v1_last = v1p + v1.GetNrows(); // last of v1
1555 
1556  const Element2 * mp = m.GetMatrixArray(); // first of m
1557  const Element2 * const m_last = mp + m.GetNoElements(); // last of m
1558 
1559  const Element3 * const v20 = v2.GetMatrixArray(); // first of v2
1560  const Element3 * v2p = v20; // running v2
1561  const Element3 * const v2_last = v2p + v2.GetNrows(); // last of v2
1562 
1563  Element1 sum = 0; // scalar result accumulator
1564  Element3 dot = 0; // M_row * v2 dot product accumulator
1565 
1566  while (v1p < v1_last) {
1567  v2p = v20; // at beginning of v2
1568  while (v2p < v2_last) { // compute (M[i] * v2) dot product
1569  dot += *mp++ * *v2p++;
1570  }
1571  sum += *v1p++ * dot; // v1[i] * (M[i] * v2)
1572  dot = 0; // start next dot product
1573  }
1574 
1575  R__ASSERT(v1p == v1_last && mp == m_last && v2p == v2_last);
1576 
1577  return sum;
1578 }
1579 
1580 ////////////////////////////////////////////////////////////////////////////////
1581 /// Modify addition: target += scalar * source.
1582 
1583 template<class Element>
1584 TVectorT<Element> &Add(TVectorT<Element> &target,Element scalar,const TVectorT<Element> &source)
1585 {
1586  if (gMatrixCheck && !AreCompatible(target,source)) {
1587  Error("Add(TVectorT<Element> &,Element,const TVectorT<Element> &)","vector's are incompatible");
1588  return target;
1589  }
1590 
1591  const Element * sp = source.GetMatrixArray();
1592  Element * tp = target.GetMatrixArray();
1593  const Element * const ftp = tp+target.GetNrows();
1594  if (scalar == 1.0 ) {
1595  while ( tp < ftp )
1596  *tp++ += *sp++;
1597  } else if (scalar == -1.0) {
1598  while ( tp < ftp )
1599  *tp++ -= *sp++;
1600  } else {
1601  while ( tp < ftp )
1602  *tp++ += scalar * *sp++;
1603  }
1604 
1605  return target;
1606 }
1607 
1608 ////////////////////////////////////////////////////////////////////////////////
1609 /// Modify addition: target += scalar * A * source.
1610 /// NOTE: in case scalar=0, do target = A * source.
1611 
1612 template<class Element>
1613 TVectorT<Element> &Add(TVectorT<Element> &target,Element scalar,
1614  const TMatrixT<Element> &a,const TVectorT<Element> &source)
1615 {
1616  if (gMatrixCheck) {
1617  R__ASSERT(target.IsValid());
1618  R__ASSERT(a.IsValid());
1619  if (a.GetNrows() != target.GetNrows() || a.GetRowLwb() != target.GetLwb()) {
1620  Error("Add(TVectorT &,const TMatrixT &,const TVectorT &)","target vector and matrix are incompatible");
1621  return target;
1622  }
1623 
1624  R__ASSERT(source.IsValid());
1625  if (a.GetNcols() != source.GetNrows() || a.GetColLwb() != source.GetLwb()) {
1626  Error("Add(TVectorT &,const TMatrixT &,const TVectorT &)","source vector and matrix are incompatible");
1627  return target;
1628  }
1629  }
1630 
1631  const Element * const sp = source.GetMatrixArray(); // sources vector ptr
1632  const Element * mp = a.GetMatrixArray(); // Matrix row ptr
1633  Element * tp = target.GetMatrixArray(); // Target vector ptr
1634 #ifdef CBLAS
1635  if (typeid(Element) == typeid(Double_t))
1636  cblas_dsymv(CblasRowMajor,CblasUpper,fNrows,scalar,mp,
1637  fNrows,sp,1,0.0,tp,1);
1638  else if (typeid(Element) != typeid(Float_t))
1639  cblas_ssymv(CblasRowMajor,CblasUpper,fNrows,scalar,mp,
1640  fNrows,sp,1,0.0,tp,1);
1641  else
1642  Error("operator*=","type %s not implemented in BLAS library",typeid(Element));
1643 #else
1644  const Element * const sp_last = sp+source.GetNrows();
1645  const Element * const tp_last = tp+target.GetNrows();
1646  if (scalar == 1.0) {
1647  while (tp < tp_last) {
1648  const Element *sp1 = sp;
1649  Element sum = 0;
1650  while (sp1 < sp_last)
1651  sum += *sp1++ * *mp++;
1652  *tp++ += sum;
1653  }
1654  } else if (scalar == 0.0) {
1655  while (tp < tp_last) {
1656  const Element *sp1 = sp;
1657  Element sum = 0;
1658  while (sp1 < sp_last)
1659  sum += *sp1++ * *mp++;
1660  *tp++ = sum;
1661  }
1662  } else if (scalar == -1.0) {
1663  while (tp < tp_last) {
1664  const Element *sp1 = sp;
1665  Element sum = 0;
1666  while (sp1 < sp_last)
1667  sum += *sp1++ * *mp++;
1668  *tp++ -= sum;
1669  }
1670  } else {
1671  while (tp < tp_last) {
1672  const Element *sp1 = sp;
1673  Element sum = 0;
1674  while (sp1 < sp_last)
1675  sum += *sp1++ * *mp++;
1676  *tp++ += scalar * sum;
1677  }
1678  }
1679 
1680  if (gMatrixCheck) R__ASSERT(mp == a.GetMatrixArray()+a.GetNoElements());
1681 #endif
1682 
1683  return target;
1684 }
1685 
1686 ////////////////////////////////////////////////////////////////////////////////
1687 /// Modify addition: target += A * source.
1688 /// NOTE: in case scalar=0, do target = A * source.
1689 
1690 template<class Element>
1691 TVectorT<Element> &Add(TVectorT<Element> &target,Element scalar,
1692  const TMatrixTSym<Element> &a,const TVectorT<Element> &source)
1693 {
1694  if (gMatrixCheck) {
1695  R__ASSERT(target.IsValid());
1696  R__ASSERT(source.IsValid());
1697  R__ASSERT(a.IsValid());
1698  if (a.GetNrows() != target.GetNrows() || a.GetRowLwb() != target.GetLwb()) {
1699  Error("Add(TVectorT &,const TMatrixT &,const TVectorT &)","target vector and matrix are incompatible");
1700  return target;
1701  }
1702  }
1703 
1704  const Element * const sp = source.GetMatrixArray(); // sources vector ptr
1705  const Element * mp = a.GetMatrixArray(); // Matrix row ptr
1706  Element * tp = target.GetMatrixArray(); // Target vector ptr
1707 #ifdef CBLAS
1708  if (typeid(Element) == typeid(Double_t))
1709  cblas_dsymv(CblasRowMajor,CblasUpper,fNrows,1.0,mp,
1710  fNrows,sp,1,0.0,tp,1);
1711  else if (typeid(Element) != typeid(Float_t))
1712  cblas_ssymv(CblasRowMajor,CblasUpper,fNrows,1.0,mp,
1713  fNrows,sp,1,0.0,tp,1);
1714  else
1715  Error("operator*=","type %s not implemented in BLAS library",typeid(Element));
1716 #else
1717  const Element * const sp_last = sp+source.GetNrows();
1718  const Element * const tp_last = tp+target.GetNrows();
1719  if (scalar == 1.0) {
1720  while (tp < tp_last) {
1721  const Element *sp1 = sp;
1722  Element sum = 0;
1723  while (sp1 < sp_last)
1724  sum += *sp1++ * *mp++;
1725  *tp++ += sum;
1726  }
1727  } else if (scalar == 0.0) {
1728  while (tp < tp_last) {
1729  const Element *sp1 = sp;
1730  Element sum = 0;
1731  while (sp1 < sp_last)
1732  sum += *sp1++ * *mp++;
1733  *tp++ = sum;
1734  }
1735  } else if (scalar == -1.0) {
1736  while (tp < tp_last) {
1737  const Element *sp1 = sp;
1738  Element sum = 0;
1739  while (sp1 < sp_last)
1740  sum += *sp1++ * *mp++;
1741  *tp++ -= sum;
1742  }
1743  } else {
1744  while (tp < tp_last) {
1745  const Element *sp1 = sp;
1746  Element sum = 0;
1747  while (sp1 < sp_last)
1748  sum += *sp1++ * *mp++;
1749  *tp++ += scalar * sum;
1750  }
1751  }
1752  R__ASSERT(mp == a.GetMatrixArray()+a.GetNoElements());
1753 #endif
1754 
1755  return target;
1756 }
1757 
1758 ////////////////////////////////////////////////////////////////////////////////
1759 /// Modify addition: target += A * source.
1760 /// NOTE: in case scalar=0, do target = A * source.
1761 
1762 template<class Element>
1763 TVectorT<Element> &Add(TVectorT<Element> &target,Element scalar,
1764  const TMatrixTSparse<Element> &a,const TVectorT<Element> &source)
1765 {
1766  if (gMatrixCheck) {
1767  R__ASSERT(target.IsValid());
1768  R__ASSERT(a.IsValid());
1769  if (a.GetNrows() != target.GetNrows() || a.GetRowLwb() != target.GetLwb()) {
1770  Error("Add(TVectorT &,const TMatrixT &,const TVectorT &)","target vector and matrix are incompatible");
1771  return target;
1772  }
1773 
1774  R__ASSERT(source.IsValid());
1775  if (a.GetNcols() != source.GetNrows() || a.GetColLwb() != source.GetLwb()) {
1776  Error("Add(TVectorT &,const TMatrixT &,const TVectorT &)","source vector and matrix are incompatible");
1777  return target;
1778  }
1779  }
1780 
1781  const Int_t * const pRowIndex = a.GetRowIndexArray();
1782  const Int_t * const pColIndex = a.GetColIndexArray();
1783  const Element * const mp = a.GetMatrixArray(); // Matrix row ptr
1784 
1785  const Element * const sp = source.GetMatrixArray(); // Source vector ptr
1786  Element * tp = target.GetMatrixArray(); // Target vector ptr
1787 
1788  if (scalar == 1.0) {
1789  for (Int_t irow = 0; irow < a.GetNrows(); irow++) {
1790  const Int_t sIndex = pRowIndex[irow];
1791  const Int_t eIndex = pRowIndex[irow+1];
1792  Element sum = 0.0;
1793  for (Int_t index = sIndex; index < eIndex; index++) {
1794  const Int_t icol = pColIndex[index];
1795  sum += mp[index]*sp[icol];
1796  }
1797  tp[irow] += sum;
1798  }
1799  } else if (scalar == 0.0) {
1800  for (Int_t irow = 0; irow < a.GetNrows(); irow++) {
1801  const Int_t sIndex = pRowIndex[irow];
1802  const Int_t eIndex = pRowIndex[irow+1];
1803  Element sum = 0.0;
1804  for (Int_t index = sIndex; index < eIndex; index++) {
1805  const Int_t icol = pColIndex[index];
1806  sum += mp[index]*sp[icol];
1807  }
1808  tp[irow] = sum;
1809  }
1810  } else if (scalar == -1.0) {
1811  for (Int_t irow = 0; irow < a.GetNrows(); irow++) {
1812  const Int_t sIndex = pRowIndex[irow];
1813  const Int_t eIndex = pRowIndex[irow+1];
1814  Element sum = 0.0;
1815  for (Int_t index = sIndex; index < eIndex; index++) {
1816  const Int_t icol = pColIndex[index];
1817  sum += mp[index]*sp[icol];
1818  }
1819  tp[irow] -= sum;
1820  }
1821  } else {
1822  for (Int_t irow = 0; irow < a.GetNrows(); irow++) {
1823  const Int_t sIndex = pRowIndex[irow];
1824  const Int_t eIndex = pRowIndex[irow+1];
1825  Element sum = 0.0;
1826  for (Int_t index = sIndex; index < eIndex; index++) {
1827  const Int_t icol = pColIndex[index];
1828  sum += mp[index]*sp[icol];
1829  }
1830  tp[irow] += scalar * sum;
1831  }
1832  }
1833 
1834  return target;
1835 }
1836 
1837 ////////////////////////////////////////////////////////////////////////////////
1838 /// Modify addition: target += scalar * ElementMult(source1,source2) .
1839 
1840 template<class Element>
1841 TVectorT<Element> &AddElemMult(TVectorT<Element> &target,Element scalar,
1842  const TVectorT<Element> &source1,const TVectorT<Element> &source2)
1843 {
1844  if (gMatrixCheck && !(AreCompatible(target,source1) && AreCompatible(target,source1))) {
1845  Error("AddElemMult(TVectorT<Element> &,Element,const TVectorT<Element> &,const TVectorT<Element> &)",
1846  "vector's are incompatible");
1847  return target;
1848  }
1849 
1850  const Element * sp1 = source1.GetMatrixArray();
1851  const Element * sp2 = source2.GetMatrixArray();
1852  Element * tp = target.GetMatrixArray();
1853  const Element * const ftp = tp+target.GetNrows();
1854 
1855  if (scalar == 1.0 ) {
1856  while ( tp < ftp )
1857  *tp++ += *sp1++ * *sp2++;
1858  } else if (scalar == -1.0) {
1859  while ( tp < ftp )
1860  *tp++ -= *sp1++ * *sp2++;
1861  } else {
1862  while ( tp < ftp )
1863  *tp++ += scalar * *sp1++ * *sp2++;
1864  }
1865 
1866  return target;
1867 }
1868 
1869 ////////////////////////////////////////////////////////////////////////////////
1870 /// Modify addition: target += scalar * ElementMult(source1,source2) only for those elements
1871 /// where select[i] != 0.0
1872 
1873 template<class Element>
1874 TVectorT<Element> &AddElemMult(TVectorT<Element> &target,Element scalar,
1875  const TVectorT<Element> &source1,const TVectorT<Element> &source2,const TVectorT<Element> &select)
1876 {
1877  if (gMatrixCheck && !( AreCompatible(target,source1) && AreCompatible(target,source1) &&
1878  AreCompatible(target,select) )) {
1879  Error("AddElemMult(TVectorT<Element> &,Element,const TVectorT<Element> &,const TVectorT<Element> &,onst TVectorT<Element> &)",
1880  "vector's are incompatible");
1881  return target;
1882  }
1883 
1884  const Element * sp1 = source1.GetMatrixArray();
1885  const Element * sp2 = source2.GetMatrixArray();
1886  const Element * mp = select.GetMatrixArray();
1887  Element * tp = target.GetMatrixArray();
1888  const Element * const ftp = tp+target.GetNrows();
1889 
1890  if (scalar == 1.0 ) {
1891  while ( tp < ftp ) {
1892  if (*mp) *tp += *sp1 * *sp2;
1893  mp++; tp++; sp1++; sp2++;
1894  }
1895  } else if (scalar == -1.0) {
1896  while ( tp < ftp ) {
1897  if (*mp) *tp -= *sp1 * *sp2;
1898  mp++; tp++; sp1++; sp2++;
1899  }
1900  } else {
1901  while ( tp < ftp ) {
1902  if (*mp) *tp += scalar * *sp1 * *sp2;
1903  mp++; tp++; sp1++; sp2++;
1904  }
1905  }
1906 
1907  return target;
1908 }
1909 
1910 ////////////////////////////////////////////////////////////////////////////////
1911 /// Modify addition: target += scalar * ElementDiv(source1,source2) .
1912 
1913 template<class Element>
1914 TVectorT<Element> &AddElemDiv(TVectorT<Element> &target,Element scalar,
1915  const TVectorT<Element> &source1,const TVectorT<Element> &source2)
1916 {
1917  if (gMatrixCheck && !(AreCompatible(target,source1) && AreCompatible(target,source1))) {
1918  Error("AddElemDiv(TVectorT<Element> &,Element,const TVectorT<Element> &,const TVectorT<Element> &)",
1919  "vector's are incompatible");
1920  return target;
1921  }
1922 
1923  const Element * sp1 = source1.GetMatrixArray();
1924  const Element * sp2 = source2.GetMatrixArray();
1925  Element * tp = target.GetMatrixArray();
1926  const Element * const ftp = tp+target.GetNrows();
1927 
1928  if (scalar == 1.0 ) {
1929  while ( tp < ftp ) {
1930  if (*sp2 != 0.0)
1931  *tp += *sp1 / *sp2;
1932  else {
1933  const Int_t irow = (sp2-source2.GetMatrixArray())/source2.GetNrows();
1934  Error("AddElemDiv","source2 (%d) is zero",irow);
1935  }
1936  tp++; sp1++; sp2++;
1937  }
1938  } else if (scalar == -1.0) {
1939  while ( tp < ftp ) {
1940  if (*sp2 != 0.0)
1941  *tp -= *sp1 / *sp2;
1942  else {
1943  const Int_t irow = (sp2-source2.GetMatrixArray())/source2.GetNrows();
1944  Error("AddElemDiv","source2 (%d) is zero",irow);
1945  }
1946  tp++; sp1++; sp2++;
1947  }
1948  } else {
1949  while ( tp < ftp ) {
1950  if (*sp2 != 0.0)
1951  *tp += scalar * *sp1 / *sp2;
1952  else {
1953  const Int_t irow = (sp2-source2.GetMatrixArray())/source2.GetNrows();
1954  Error("AddElemDiv","source2 (%d) is zero",irow);
1955  }
1956  tp++; sp1++; sp2++;
1957  }
1958  }
1959 
1960  return target;
1961 }
1962 
1963 ////////////////////////////////////////////////////////////////////////////////
1964 /// Modify addition: target += scalar * ElementDiv(source1,source2) only for those elements
1965 /// where select[i] != 0.0
1966 
1967 template<class Element>
1968 TVectorT<Element> &AddElemDiv(TVectorT<Element> &target,Element scalar,
1969  const TVectorT<Element> &source1,const TVectorT<Element> &source2,const TVectorT<Element> &select)
1970 {
1971  if (gMatrixCheck && !( AreCompatible(target,source1) && AreCompatible(target,source1) &&
1972  AreCompatible(target,select) )) {
1973  Error("AddElemDiv(TVectorT<Element> &,Element,const TVectorT<Element> &,const TVectorT<Element> &,onst TVectorT<Element> &)",
1974  "vector's are incompatible");
1975  return target;
1976  }
1977 
1978  const Element * sp1 = source1.GetMatrixArray();
1979  const Element * sp2 = source2.GetMatrixArray();
1980  const Element * mp = select.GetMatrixArray();
1981  Element * tp = target.GetMatrixArray();
1982  const Element * const ftp = tp+target.GetNrows();
1983 
1984  if (scalar == 1.0 ) {
1985  while ( tp < ftp ) {
1986  if (*mp) {
1987  if (*sp2 != 0.0)
1988  *tp += *sp1 / *sp2;
1989  else {
1990  const Int_t irow = (sp2-source2.GetMatrixArray())/source2.GetNrows();
1991  Error("AddElemDiv","source2 (%d) is zero",irow);
1992  }
1993  }
1994  mp++; tp++; sp1++; sp2++;
1995  }
1996  } else if (scalar == -1.0) {
1997  while ( tp < ftp ) {
1998  if (*mp) {
1999  if (*sp2 != 0.0)
2000  *tp -= *sp1 / *sp2;
2001  else {
2002  const Int_t irow = (sp2-source2.GetMatrixArray())/source2.GetNrows();
2003  Error("AddElemDiv","source2 (%d) is zero",irow);
2004  }
2005  }
2006  mp++; tp++; sp1++; sp2++;
2007  }
2008  } else {
2009  while ( tp < ftp ) {
2010  if (*mp) {
2011  if (*sp2 != 0.0)
2012  *tp += scalar * *sp1 / *sp2;
2013  else {
2014  const Int_t irow = (sp2-source2.GetMatrixArray())/source2.GetNrows();
2015  Error("AddElemDiv","source2 (%d) is zero",irow);
2016  }
2017  }
2018  mp++; tp++; sp1++; sp2++;
2019  }
2020  }
2021 
2022  return target;
2023 }
2024 
2025 ////////////////////////////////////////////////////////////////////////////////
2026 /// Multiply target by the source, element-by-element.
2027 
2028 template<class Element>
2029 TVectorT<Element> &ElementMult(TVectorT<Element> &target,const TVectorT<Element> &source)
2030 {
2031  if (gMatrixCheck && !AreCompatible(target,source)) {
2032  Error("ElementMult(TVectorT<Element> &,const TVectorT<Element> &)","vector's are incompatible");
2033  return target;
2034  }
2035 
2036  const Element * sp = source.GetMatrixArray();
2037  Element * tp = target.GetMatrixArray();
2038  const Element * const ftp = tp+target.GetNrows();
2039  while ( tp < ftp )
2040  *tp++ *= *sp++;
2041 
2042  return target;
2043 }
2044 
2045 ////////////////////////////////////////////////////////////////////////////////
2046 /// Multiply target by the source, element-by-element only where select[i] != 0.0
2047 
2048 template<class Element>
2049 TVectorT<Element> &ElementMult(TVectorT<Element> &target,const TVectorT<Element> &source,const TVectorT<Element> &select)
2050 {
2051  if (gMatrixCheck && !(AreCompatible(target,source) && AreCompatible(target,select))) {
2052  Error("ElementMult(TVectorT<Element> &,const TVectorT<Element> &,const TVectorT<Element> &)","vector's are incompatible");
2053  return target;
2054  }
2055 
2056  const Element * sp = source.GetMatrixArray();
2057  const Element * mp = select.GetMatrixArray();
2058  Element * tp = target.GetMatrixArray();
2059  const Element * const ftp = tp+target.GetNrows();
2060  while ( tp < ftp ) {
2061  if (*mp) *tp *= *sp;
2062  mp++; tp++; sp++;
2063  }
2064 
2065  return target;
2066 }
2067 
2068 ////////////////////////////////////////////////////////////////////////////////
2069 /// Divide target by the source, element-by-element.
2070 
2071 template<class Element>
2072 TVectorT<Element> &ElementDiv(TVectorT<Element> &target,const TVectorT<Element> &source)
2073 {
2074  if (gMatrixCheck && !AreCompatible(target,source)) {
2075  Error("ElementDiv(TVectorT<Element> &,const TVectorT<Element> &)","vector's are incompatible");
2076  return target;
2077  }
2078 
2079  const Element * sp = source.GetMatrixArray();
2080  Element * tp = target.GetMatrixArray();
2081  const Element * const ftp = tp+target.GetNrows();
2082  while ( tp < ftp ) {
2083  if (*sp != 0.0)
2084  *tp++ /= *sp++;
2085  else {
2086  const Int_t irow = (sp-source.GetMatrixArray())/source.GetNrows();
2087  Error("ElementDiv","source (%d) is zero",irow);
2088  }
2089  }
2090 
2091  return target;
2092 }
2093 
2094 ////////////////////////////////////////////////////////////////////////////////
2095 /// Divide target by the source, element-by-element only where select[i] != 0.0
2096 
2097 template<class Element>
2098 TVectorT<Element> &ElementDiv(TVectorT<Element> &target,const TVectorT<Element> &source,const TVectorT<Element> &select)
2099 {
2100  if (gMatrixCheck && !AreCompatible(target,source)) {
2101  Error("ElementDiv(TVectorT<Element> &,const TVectorT<Element> &,const TVectorT<Element> &)","vector's are incompatible");
2102  return target;
2103  }
2104 
2105  const Element * sp = source.GetMatrixArray();
2106  const Element * mp = select.GetMatrixArray();
2107  Element * tp = target.GetMatrixArray();
2108  const Element * const ftp = tp+target.GetNrows();
2109  while ( tp < ftp ) {
2110  if (*mp) {
2111  if (*sp != 0.0)
2112  *tp /= *sp;
2113  else {
2114  const Int_t irow = (sp-source.GetMatrixArray())/source.GetNrows();
2115  Error("ElementDiv","source (%d) is zero",irow);
2116  }
2117  }
2118  mp++; tp++; sp++;
2119  }
2120 
2121  return target;
2122 }
2123 
2124 ////////////////////////////////////////////////////////////////////////////////
2125 /// Check if v1 and v2 are both valid and have the same shape
2126 
2127 template<class Element1,class Element2>
2128 Bool_t AreCompatible(const TVectorT<Element1> &v1,const TVectorT<Element2> &v2,Int_t verbose)
2129 {
2130  if (!v1.IsValid()) {
2131  if (verbose)
2132  ::Error("AreCompatible", "vector 1 not valid");
2133  return kFALSE;
2134  }
2135  if (!v2.IsValid()) {
2136  if (verbose)
2137  ::Error("AreCompatible", "vector 2 not valid");
2138  return kFALSE;
2139  }
2140 
2141  if (v1.GetNrows() != v2.GetNrows() || v1.GetLwb() != v2.GetLwb()) {
2142  if (verbose)
2143  ::Error("AreCompatible", "matrices 1 and 2 not compatible");
2144  return kFALSE;
2145  }
2146 
2147  return kTRUE;
2148 }
2149 
2150 ////////////////////////////////////////////////////////////////////////////////
2151 /// Check if m and v are both valid and have compatible shapes for M * v
2152 
2153 template<class Element1, class Element2>
2154 Bool_t AreCompatible(const TMatrixT<Element1> &m,const TVectorT<Element2> &v,Int_t verbose)
2155 {
2156  if (!m.IsValid()) {
2157  if (verbose)
2158  ::Error("AreCompatible", "Matrix not valid");
2159  return kFALSE;
2160  }
2161  if (!v.IsValid()) {
2162  if (verbose)
2163  ::Error("AreCompatible", "vector not valid");
2164  return kFALSE;
2165  }
2166 
2167  if (m.GetNcols() != v.GetNrows() ) {
2168  if (verbose)
2169  ::Error("AreCompatible", "matrix and vector not compatible");
2170  return kFALSE;
2171  }
2172 
2173  return kTRUE;
2174 }
2175 
2176 ////////////////////////////////////////////////////////////////////////////////
2177 /// Check if m and v are both valid and have compatible shapes for v * M
2178 
2179 template<class Element1, class Element2>
2180 Bool_t AreCompatible(const TVectorT<Element1> &v,const TMatrixT<Element2> &m,Int_t verbose)
2181 {
2182  if (!m.IsValid()) {
2183  if (verbose)
2184  ::Error("AreCompatible", "Matrix not valid");
2185  return kFALSE;
2186  }
2187  if (!v.IsValid()) {
2188  if (verbose)
2189  ::Error("AreCompatible", "vector not valid");
2190  return kFALSE;
2191  }
2192 
2193  if (v.GetNrows() != m.GetNrows() ) {
2194  if (verbose)
2195  ::Error("AreCompatible", "vector and matrix not compatible");
2196  return kFALSE;
2197  }
2198 
2199  return kTRUE;
2200 }
2201 
2202 ////////////////////////////////////////////////////////////////////////////////
2203 /// Compare two vectors and print out the result of the comparison.
2204 
2205 template<class Element>
2206 void Compare(const TVectorT<Element> &v1,const TVectorT<Element> &v2)
2207 {
2208  if (!AreCompatible(v1,v2)) {
2209  Error("Compare(const TVectorT<Element> &,const TVectorT<Element> &)","vectors are incompatible");
2210  return;
2211  }
2212 
2213  printf("\n\nComparison of two TVectorTs:\n");
2214 
2215  Element norm1 = 0; // Norm of the Matrices
2216  Element norm2 = 0; // Norm of the Matrices
2217  Element ndiff = 0; // Norm of the difference
2218  Int_t imax = 0; // For the elements that differ most
2219  Element difmax = -1;
2220  const Element *mp1 = v1.GetMatrixArray(); // Vector element pointers
2221  const Element *mp2 = v2.GetMatrixArray();
2222 
2223  for (Int_t i = 0; i < v1.GetNrows(); i++) {
2224  const Element mv1 = *mp1++;
2225  const Element mv2 = *mp2++;
2226  const Element diff = TMath::Abs(mv1-mv2);
2227 
2228  if (diff > difmax) {
2229  difmax = diff;
2230  imax = i;
2231  }
2232  norm1 += TMath::Abs(mv1);
2233  norm2 += TMath::Abs(mv2);
2234  ndiff += TMath::Abs(diff);
2235  }
2236 
2237  imax += v1.GetLwb();
2238  printf("\nMaximal discrepancy \t\t%g",difmax);
2239  printf("\n occured at the point\t\t(%d)",imax);
2240  const Element mv1 = v1(imax);
2241  const Element mv2 = v2(imax);
2242  printf("\n Vector 1 element is \t\t%g",mv1);
2243  printf("\n Vector 2 element is \t\t%g",mv2);
2244  printf("\n Absolute error v2[i]-v1[i]\t\t%g",mv2-mv1);
2245  printf("\n Relative error\t\t\t\t%g\n",
2246  (mv2-mv1)/TMath::Max(TMath::Abs(mv2+mv1)/2,(Element)1e-7));
2247 
2248  printf("\n||Vector 1|| \t\t\t%g",norm1);
2249  printf("\n||Vector 2|| \t\t\t%g",norm2);
2250  printf("\n||Vector1-Vector2||\t\t\t\t%g",ndiff);
2251  printf("\n||Vector1-Vector2||/sqrt(||Vector1|| ||Vector2||)\t%g\n\n",
2252  ndiff/TMath::Max(TMath::Sqrt(norm1*norm2),1e-7));
2253 }
2254 
2255 ////////////////////////////////////////////////////////////////////////////////
2256 /// Validate that all elements of vector have value val within maxDevAllow .
2257 
2258 template<class Element>
2259 Bool_t VerifyVectorValue(const TVectorT<Element> &v,Element val,
2260  Int_t verbose,Element maxDevAllow)
2261 {
2262  Int_t imax = 0;
2263  Element maxDevObs = 0;
2264 
2265  if (TMath::Abs(maxDevAllow) <= 0.0)
2266  maxDevAllow = std::numeric_limits<Element>::epsilon();
2267 
2268  for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++) {
2269  const Element dev = TMath::Abs(v(i)-val);
2270  if (dev > maxDevObs) {
2271  imax = i;
2272  maxDevObs = dev;
2273  }
2274  }
2275 
2276  if (maxDevObs == 0)
2277  return kTRUE;
2278 
2279  if (verbose) {
2280  printf("Largest dev for (%d); dev = |%g - %g| = %g\n",imax,v(imax),val,maxDevObs);
2281  if (maxDevObs > maxDevAllow)
2282  Error("VerifyVectorValue","Deviation > %g\n",maxDevAllow);
2283  }
2284 
2285  if (maxDevObs > maxDevAllow)
2286  return kFALSE;
2287  return kTRUE;
2288 }
2289 
2290 ////////////////////////////////////////////////////////////////////////////////
2291 /// Verify that elements of the two vectors are equal within maxDevAllow .
2292 
2293 template<class Element>
2294 Bool_t VerifyVectorIdentity(const TVectorT<Element> &v1,const TVectorT<Element> &v2,
2295  Int_t verbose, Element maxDevAllow)
2296 {
2297  Int_t imax = 0;
2298  Element maxDevObs = 0;
2299 
2300  if (!AreCompatible(v1,v2))
2301  return kFALSE;
2302 
2303  if (TMath::Abs(maxDevAllow) <= 0.0)
2304  maxDevAllow = std::numeric_limits<Element>::epsilon();
2305 
2306  for (Int_t i = v1.GetLwb(); i <= v1.GetUpb(); i++) {
2307  const Element dev = TMath::Abs(v1(i)-v2(i));
2308  if (dev > maxDevObs) {
2309  imax = i;
2310  maxDevObs = dev;
2311  }
2312  }
2313 
2314  if (maxDevObs == 0)
2315  return kTRUE;
2316 
2317  if (verbose) {
2318  printf("Largest dev for (%d); dev = |%g - %g| = %g\n",imax,v1(imax),v2(imax),maxDevObs);
2319  if (maxDevObs > maxDevAllow)
2320  Error("VerifyVectorIdentity","Deviation > %g\n",maxDevAllow);
2321  }
2322 
2323  if (maxDevObs > maxDevAllow) {
2324  return kFALSE;
2325  }
2326  return kTRUE;
2327 }
2328 
2329 ////////////////////////////////////////////////////////////////////////////////
2330 /// Stream an object of class TVectorT.
2331 
2332 template<class Element>
2333 void TVectorT<Element>::Streamer(TBuffer &R__b)
2334 {
2335  if (R__b.IsReading()) {
2336  UInt_t R__s, R__c;
2337  Version_t R__v = R__b.ReadVersion(&R__s,&R__c);
2338  if (R__v > 1) {
2339  Clear();
2340  R__b.ReadClassBuffer(TVectorT<Element>::Class(),this,R__v,R__s,R__c);
2341  } else { //====process old versions before automatic schema evolution
2342  TObject::Streamer(R__b);
2343  R__b >> fRowLwb;
2344  fNrows = R__b.ReadArray(fElements);
2345  R__b.CheckByteCount(R__s, R__c, TVectorT<Element>::IsA());
2346  }
2347  if (fNrows > 0 && fNrows <= kSizeMax) {
2348  memcpy(fDataStack,fElements,fNrows*sizeof(Element));
2349  delete [] fElements;
2350  fElements = fDataStack;
2351  } else if (fNrows < 0)
2352  Invalidate();
2353 
2354  if (R__v < 3)
2355  MakeValid();
2356  } else {
2357  R__b.WriteClassBuffer(TVectorT<Element>::Class(),this);
2358  }
2359 }
2360 
2361 #ifndef ROOT_TMatrixFfwd
2362 #include "TMatrixFfwd.h"
2363 #endif
2364 #ifndef ROOT_TMatrixFSymfwd
2365 #include "TMatrixFSymfwd.h"
2366 #endif
2367 #ifndef ROOT_TMatrixFSparsefwd
2368 #include "TMatrixFSparsefwd.h"
2369 #endif
2370 
2371 template class TVectorT<Float_t>;
2372 
2373 template Bool_t operator== <Float_t>(const TVectorF &source1,const TVectorF &source2);
2374 template TVectorF operator+ <Float_t>(const TVectorF &source1,const TVectorF &source2);
2375 template TVectorF operator- <Float_t>(const TVectorF &source1,const TVectorF &source2);
2376 template Float_t operator* <Float_t>(const TVectorF &source1,const TVectorF &source2);
2377 template TVectorF operator* <Float_t>(const TMatrixF &a, const TVectorF &source);
2378 template TVectorF operator* <Float_t>(const TMatrixFSym &a, const TVectorF &source);
2379 template TVectorF operator* <Float_t>(const TMatrixFSparse &a, const TVectorF &source);
2380 template TVectorF operator* <Float_t>( Float_t val, const TVectorF &source);
2381 
2382 template Float_t Dot <Float_t>(const TVectorF &v1, const TVectorF &v2);
2383 template TMatrixF OuterProduct <Float_t,Float_t>
2384  (const TVectorF &v1, const TVectorF &v2);
2385 template TMatrixF &OuterProduct <Float_t,Float_t,Float_t>
2386  ( TMatrixF &target, const TVectorF &v1, const TVectorF &v2);
2387 template Float_t Mult <Float_t,Float_t,Float_t>
2388  (const TVectorF &v1, const TMatrixF &m, const TVectorF &v2);
2389 
2390 template TVectorF &Add <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source);
2391 template TVectorF &Add <Float_t>( TVectorF &target, Float_t scalar, const TMatrixF &a,
2392  const TVectorF &source);
2393 template TVectorF &Add <Float_t>( TVectorF &target, Float_t scalar, const TMatrixFSym &a,
2394  const TVectorF &source);
2395 template TVectorF &Add <Float_t>( TVectorF &target, Float_t scalar, const TMatrixFSparse &a,
2396  const TVectorF &source);
2397 template TVectorF &AddElemMult <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source1,
2398  const TVectorF &source2);
2399 template TVectorF &AddElemMult <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source1,
2400  const TVectorF &source2,const TVectorF &select);
2401 template TVectorF &AddElemDiv <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source1,
2402  const TVectorF &source2);
2403 template TVectorF &AddElemDiv <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source1,
2404  const TVectorF &source2,const TVectorF &select);
2405 template TVectorF &ElementMult <Float_t>( TVectorF &target, const TVectorF &source);
2406 template TVectorF &ElementMult <Float_t>( TVectorF &target, const TVectorF &source, const TVectorF &select);
2407 template TVectorF &ElementDiv <Float_t>( TVectorF &target, const TVectorF &source);
2408 template TVectorF &ElementDiv <Float_t>( TVectorF &target, const TVectorF &source, const TVectorF &select);
2409 
2410 template Bool_t AreCompatible <Float_t,Float_t> (const TVectorF &v1,const TVectorF &v2,Int_t verbose);
2411 template Bool_t AreCompatible <Float_t,Double_t>(const TVectorF &v1,const TVectorD &v2,Int_t verbose);
2412 template Bool_t AreCompatible <Float_t,Float_t> (const TMatrixF &m, const TVectorF &v, Int_t verbose);
2413 template Bool_t AreCompatible <Float_t,Float_t> (const TVectorF &v, const TMatrixF &m, Int_t verbose);
2414 
2415 template void Compare <Float_t> (const TVectorF &v1,const TVectorF &v2);
2416 template Bool_t VerifyVectorValue <Float_t> (const TVectorF &m, Float_t val,Int_t verbose,Float_t maxDevAllow);
2417 template Bool_t VerifyVectorIdentity<Float_t> (const TVectorF &m1,const TVectorF &m2, Int_t verbose,Float_t maxDevAllow);
2418 
2419 #ifndef ROOT_TMatrixDfwd
2420 #include "TMatrixDfwd.h"
2421 #endif
2422 #ifndef ROOT_TMatrixDSymfwd
2423 #include "TMatrixDSymfwd.h"
2424 #endif
2425 #ifndef ROOT_TMatrixDSparsefwd
2426 #include "TMatrixDSparsefwd.h"
2427 #endif
2428 
2429 template class TVectorT<Double_t>;
2430 
2431 template Bool_t operator== <Double_t>(const TVectorD &source1,const TVectorD &source2);
2432 template TVectorD operator+ <Double_t>(const TVectorD &source1,const TVectorD &source2);
2433 template TVectorD operator- <Double_t>(const TVectorD &source1,const TVectorD &source2);
2434 template Double_t operator* <Double_t>(const TVectorD &source1,const TVectorD &source2);
2435 template TVectorD operator* <Double_t>(const TMatrixD &a, const TVectorD &source);
2436 template TVectorD operator* <Double_t>(const TMatrixDSym &a, const TVectorD &source);
2437 template TVectorD operator* <Double_t>(const TMatrixDSparse &a, const TVectorD &source);
2438 template TVectorD operator* <Double_t>( Double_t val, const TVectorD &source);
2439 
2440 template Double_t Dot <Double_t>(const TVectorD &v1, const TVectorD &v2);
2441 template TMatrixD OuterProduct <Double_t,Double_t>
2442  (const TVectorD &v1, const TVectorD &v2);
2443 template TMatrixD &OuterProduct <Double_t,Double_t,Double_t>
2444  ( TMatrixD &target, const TVectorD &v1, const TVectorD &v2);
2445 template Double_t Mult <Double_t,Double_t,Double_t>
2446  (const TVectorD &v1, const TMatrixD &m, const TVectorD &v2);
2447 
2448 template TVectorD &Add <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source);
2449 template TVectorD &Add <Double_t>( TVectorD &target, Double_t scalar, const TMatrixD &a,
2450  const TVectorD &source);
2451 template TVectorD &Add <Double_t>( TVectorD &target, Double_t scalar, const TMatrixDSym &a
2452  , const TVectorD &source);
2453 template TVectorD &Add <Double_t>( TVectorD &target, Double_t scalar, const TMatrixDSparse &a
2454  , const TVectorD &source);
2455 template TVectorD &AddElemMult <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source1,
2456  const TVectorD &source2);
2457 template TVectorD &AddElemMult <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source1,
2458  const TVectorD &source2,const TVectorD &select);
2459 template TVectorD &AddElemDiv <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source1,
2460  const TVectorD &source2);
2461 template TVectorD &AddElemDiv <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source1,
2462  const TVectorD &source2,const TVectorD &select);
2463 template TVectorD &ElementMult <Double_t>( TVectorD &target, const TVectorD &source);
2464 template TVectorD &ElementMult <Double_t>( TVectorD &target, const TVectorD &source, const TVectorD &select);
2465 template TVectorD &ElementDiv <Double_t>( TVectorD &target, const TVectorD &source);
2466 template TVectorD &ElementDiv <Double_t>( TVectorD &target, const TVectorD &source, const TVectorD &select);
2467 
2468 template Bool_t AreCompatible <Double_t,Double_t>(const TVectorD &v1,const TVectorD &v2,Int_t verbose);
2469 template Bool_t AreCompatible <Double_t,Float_t> (const TVectorD &v1,const TVectorF &v2,Int_t verbose);
2470 template Bool_t AreCompatible <Double_t,Double_t>(const TMatrixD &m, const TVectorD &v, Int_t verbose);
2471 template Bool_t AreCompatible <Double_t,Double_t>(const TVectorD &v, const TMatrixD &m, Int_t verbose);
2472 
2473 template void Compare <Double_t> (const TVectorD &v1,const TVectorD &v2);
2474 template Bool_t VerifyVectorValue <Double_t> (const TVectorD &m, Double_t val,Int_t verbose,Double_t maxDevAllow);
2475 template Bool_t VerifyVectorIdentity<Double_t> (const TVectorD &m1,const TVectorD &m2, Int_t verbose,Double_t maxDevAllow);
TMatrixTBase::GetRowLwb
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:132
TVectorT::Sum
Element Sum() const
Compute sum of elements.
Definition: TVectorT.cxx:634
TVectorT::ResizeTo
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:291
TBuffer::WriteClassBuffer
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
TVectorT::GetUpb
Int_t GetUpb() const
Definition: TVectorT.h:80
template
Small helper to encapsulate whether to return the value pointed to by the iterator or its address...
TMath::LocMax
Long64_t LocMax(Long64_t n, const T *a)
Definition: TMath.h:724
v1
const Double_t * v1
Definition: TArcBall.cxx:33
TVectorT::operator-=
TVectorT< Element > & operator-=(Element val)
Subtract val from every element of the vector.
Definition: TVectorT.cxx:876
TBuffer::IsReading
Bool_t IsReading() const
Definition: TBuffer.h:81
Version_t
short Version_t
Definition: RtypesCore.h:61
TVectorT::Max
Element Max() const
return maximum vector element value
Definition: TVectorT.cxx:663
TVectorT::IsValid
Bool_t IsValid() const
Definition: TVectorT.h:89
Float_t
float Float_t
Definition: RtypesCore.h:53
TMatrixTSparse::GetMatrixArray
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSparse.h:208
TMatrixTBase::GetRowUpb
Int_t GetRowUpb() const
Definition: TMatrixTBase.h:133
TVectorT::Add
void Add(const TVectorT< Element > &v)
Add vector v to this vector.
Definition: TVectorT.cxx:81
Option_t
const char Option_t
Definition: RtypesCore.h:62
TMatrixTColumn_const
Definition: TMatrixDUtilsfwd.h:38
Dot< Double_t >
template Double_t Dot< Double_t >(const TVectorD &v1, const TVectorD &v2)
TVectorT::Draw
void Draw(Option_t *option="")
Draw this vector The histogram is named "TVectorT" by default and no title.
Definition: TVectorT.cxx:1350
Compare< Double_t >
template void Compare< Double_t >(const TVectorD &v1, const TVectorD &v2)
TVectorT::GetNrows
Int_t GetNrows() const
Definition: TVectorT.h:81
AddElemDiv
TVectorT< Element > & AddElemDiv(TVectorT< Element > &target, Element scalar, const TVectorT< Element > &source1, const TVectorT< Element > &source2)
Modify addition: target += scalar * ElementDiv(source1,source2) .
Definition: TVectorT.cxx:1914
OuterProduct
TMatrixT< Element1 > OuterProduct(const TVectorT< Element1 > &v1, const TVectorT< Element2 > &v2)
Return the matrix M = v1 * v2'.
Definition: TVectorT.cxx:1492
TString::ToUpper
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1101
TBuffer
Buffer base class used for serializing objects.
Definition: TBuffer.h:40
AddElemMult< Double_t >
template TVectorD & AddElemMult< Double_t >(TVectorD &target, Double_t scalar, const TVectorD &source1, const TVectorD &source2)
R__ASSERT
#define R__ASSERT(e)
Definition: TError.h:98
TMatrixTBase::GetNoElements
Int_t GetNoElements() const
Definition: TMatrixTBase.h:138
TMatrixTSparseDiag_const::GetMatrix
const TMatrixTBase< Element > * GetMatrix() const
Definition: TMatrixTUtils.h:663
TBuffer::CheckByteCount
virtual Int_t CheckByteCount(UInt_t startpos, UInt_t bcnt, const TClass *clss)=0
gROOT
#define gROOT
Definition: TROOT.h:340
TVectorT::Invert
TVectorT< Element > & Invert()
v[i] = 1/v[i]
Definition: TVectorT.cxx:519
TVectorT::GetSub
TVectorT< Element > & GetSub(Int_t row_lwb, Int_t row_upb, TVectorT< Element > &target, Option_t *option="S") const
Get subvector [row_lwb..row_upb]; The indexing range of the returned vector depends on the argument o...
Definition: TVectorT.cxx:370
TString
Basic string class.
Definition: TString.h:137
ElementDiv< Float_t >
template TVectorF & ElementDiv< Float_t >(TVectorF &target, const TVectorF &source)
Mult< Float_t, Float_t, Float_t >
template Float_t Mult< Float_t, Float_t, Float_t >(const TVectorF &v1, const TMatrixF &m, const TVectorF &v2)
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:170
Int_t
int Int_t
Definition: RtypesCore.h:41
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
a
TArc * a
Definition: textangle.C:12
kFALSE
const Bool_t kFALSE
Definition: Rtypes.h:92
operator-
TVectorT< Element > operator-(const TVectorT< Element > &source1, const TVectorT< Element > &source2)
Return source1-source2.
Definition: TVectorT.cxx:1420
VerifyVectorIdentity
Bool_t VerifyVectorIdentity(const TVectorT< Element > &v1, const TVectorT< Element > &v2, Int_t verbose, Element maxDevAllow)
Verify that elements of the two vectors are equal within maxDevAllow .
Definition: TVectorT.cxx:2294
AreCompatible< Float_t, Double_t >
template Bool_t AreCompatible< Float_t, Double_t >(const TVectorF &v1, const TVectorD &v2, Int_t verbose)
ElementDiv< Double_t >
template TVectorD & ElementDiv< Double_t >(TVectorD &target, const TVectorD &source)
AddElemDiv< Float_t >
template TVectorF & AddElemDiv< Float_t >(TVectorF &target, Float_t scalar, const TVectorF &source1, const TVectorF &source2)
TMatrixT::ResizeTo
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1201
TVectorT::operator>
Bool_t operator>(Element val) const
Are all vector elements > val?
Definition: TVectorT.cxx:1201
TVectorT::Sqr
TVectorT< Element > & Sqr()
Square each element of the vector.
Definition: TVectorT.cxx:479
TMatrixTSparseRow_const::GetNindex
Int_t GetNindex() const
Definition: TMatrixTUtils.h:603
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
TMatrixTRow_const::GetMatrix
const TMatrixTBase< Element > * GetMatrix() const
Definition: TMatrixTUtils.h:132
ROOT::Math::beta
double beta(double x, double y)
Calculates the beta function.
Definition: SpecFuncMathCore.cxx:111
TVectorT::Randomize
void Randomize(Element alpha, Element beta, Double_t &seed)
randomize vector elements value
Definition: TVectorT.cxx:1302
Mult
Element1 Mult(const TVectorT< Element1 > &v1, const TMatrixT< Element2 > &m, const TVectorT< Element3 > &v2)
Perform v1 * M * v2, a scalar result.
Definition: TVectorT.cxx:1539
TMatrixTBase::IsValid
Bool_t IsValid() const
Definition: TMatrixTBase.h:157
OuterProduct< Float_t, Float_t >
template TMatrixF OuterProduct< Float_t, Float_t >(const TVectorF &v1, const TVectorF &v2)
TVectorT::Use
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:346
RootCsg::dot
Double_t dot(const TVector2 &v1, const TVector2 &v2)
Definition: CsgOps.cxx:333
Add< Double_t >
template TVectorD & Add< Double_t >(TVectorD &target, Double_t scalar, const TVectorD &source)
TElementActionT::Operation
virtual void Operation(Element &element) const =0
AddElemMult< Float_t >
template TVectorF & AddElemMult< Float_t >(TVectorF &target, Float_t scalar, const TVectorF &source1, const TVectorF &source2)
TVectorT::SetSub
TVectorT< Element > & SetSub(Int_t row_lwb, const TVectorT< Element > &source)
Insert vector source starting at [row_lwb], thereby overwriting the part [row_lwb..row_lwb+nrows_source];.
Definition: TVectorT.cxx:419
Compare< Float_t >
template void Compare< Float_t >(const TVectorF &v1, const TVectorF &v2)
operator*
Element operator*(const TVectorT< Element > &v1, const TVectorT< Element > &v2)
Compute the scalar product.
Definition: TVectorT.cxx:1393
TMatrixTSparseRow_const
Definition: TMatrixDUtilsfwd.h:42
Dot
Element Dot(const TVectorT< Element > &v1, const TVectorT< Element > &v2)
return inner-produvt v1 . v2
Definition: TVectorT.cxx:1475
va_
#define va_(arg)
Definition: Varargs.h:41
TMatrixTRow_const::GetInc
Int_t GetInc() const
Definition: TMatrixTUtils.h:134
TObject::operator=
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.cxx:102
Add< Float_t >
template TVectorF & Add< Float_t >(TVectorF &target, Float_t scalar, const TVectorF &source)
TVectorT::Sqrt
TVectorT< Element > & Sqrt()
Take square root of all elements.
Definition: TVectorT.cxx:497
OuterProduct< Float_t, Float_t, Float_t >
template TMatrixF & OuterProduct< Float_t, Float_t, Float_t >(TMatrixF &target, const TVectorF &v1, const TVectorF &v2)
Error
void Error(const char *location, const char *msgfmt,...)
TMatrixTBase::GetColUpb
Int_t GetColUpb() const
Definition: TMatrixTBase.h:136
TMatrixTSym::GetMatrixArray
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSym.h:183
ElementMult
TVectorT< Element > & ElementMult(TVectorT< Element > &target, const TVectorT< Element > &source)
Multiply target by the source, element-by-element.
Definition: TVectorT.cxx:2029
TVectorT::Allocate
void Allocate(Int_t nrows, Int_t row_lwb=0, Int_t init=0)
Allocate new vector.
Definition: TVectorT.cxx:148
TMatrixTBase::GetNrows
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
TVectorT::GetMatrixArray
Element * GetMatrixArray()
Definition: TVectorT.h:84
TMatrixTSparseRow_const::GetDataPtr
const Element * GetDataPtr() const
Definition: TMatrixTUtils.h:600
TVectorT::Print
void Print(Option_t *option="") const
Print the vector as a list of elements.
Definition: TVectorT.cxx:1360
AreCompatible< Double_t, Double_t >
template Bool_t AreCompatible< Double_t, Double_t >(const TVectorD &v1, const TVectorD &v2, Int_t verbose)
TVectorT::operator==
Bool_t operator==(Element val) const
Are all vector elements equal to val?
Definition: TVectorT.cxx:1133
VerifyVectorIdentity< Float_t >
template Bool_t VerifyVectorIdentity< Float_t >(const TVectorF &m1, const TVectorF &m2, Int_t verbose, Float_t maxDevAllow)
TVectorT::operator>=
Bool_t operator>=(Element val) const
Are all vector elements >= val?
Definition: TVectorT.cxx:1218
TVectorT::NonZeros
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TVectorT.cxx:617
Add
TVectorT< Element > & Add(TVectorT< Element > &target, Element scalar, const TVectorT< Element > &source)
Modify addition: target += scalar * source.
Definition: TVectorT.cxx:1584
VerifyVectorIdentity< Double_t >
template Bool_t VerifyVectorIdentity< Double_t >(const TVectorD &m1, const TVectorD &m2, Int_t verbose, Double_t maxDevAllow)
TVectorT::Min
Element Min() const
return minimum vector element value
Definition: TVectorT.cxx:651
TMatrixTSparse::GetRowIndexArray
virtual const Int_t * GetRowIndexArray() const
Definition: TMatrixTSparse.h:210
operator==
Bool_t operator==(const TVectorT< Element > &v1, const TVectorT< Element > &v2)
Check to see if two vectors are identical.
Definition: TVectorT.cxx:1383
v
SVector< double, 2 > v
Definition: Dict.h:5
TVectorT::NormInf
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition: TVectorT.cxx:600
TVectorT::operator!=
Bool_t operator!=(Element val) const
Are all vector elements not equal to val?
Definition: TVectorT.cxx:1150
TMatrixTDiag_const::GetMatrix
const TMatrixTBase< Element > * GetMatrix() const
Definition: TMatrixTUtils.h:330
TVectorT::operator=
TVectorT< Element > & operator=(const TVectorT< Element > &source)
Notice that this assignment does NOT change the ownership : if the storage space was adopted...
Definition: TVectorT.cxx:677
TMatrixTSparseRow_const::GetColPtr
const Int_t * GetColPtr() const
Definition: TMatrixTUtils.h:601
TVectorT::operator*=
TVectorT< Element > & operator*=(Element val)
Multiply every element of the vector with val.
Definition: TVectorT.cxx:892
TVectorT::Zero
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:450
UInt_t
unsigned int UInt_t
Definition: RtypesCore.h:42
m
TMarker * m
Definition: textangle.C:8
verbose
bool verbose
Definition: binarySearchTime.cxx:27
Form
char * Form(const char *fmt,...)
TMatrixTDiag_const
Definition: TMatrixDUtilsfwd.h:39
TMatrixT::GetMatrixArray
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:211
OuterProduct< Double_t, Double_t, Double_t >
template TMatrixD & OuterProduct< Double_t, Double_t, Double_t >(TMatrixD &target, const TVectorD &v1, const TVectorD &v2)
TMatrixTColumn_const::GetMatrix
const TMatrixTBase< Element > * GetMatrix() const
Definition: TMatrixTUtils.h:231
gMatrixCheck
R__EXTERN Int_t gMatrixCheck
Definition: TMatrixTBase.h:91
TVectorT::Norm2Sqr
Element Norm2Sqr() const
Compute the square of the 2-norm SUM{ v[i]^2 }.
Definition: TVectorT.cxx:581
operator+
TVectorT< Element > operator+(const TVectorT< Element > &source1, const TVectorT< Element > &source2)
Return source1+source2.
Definition: TVectorT.cxx:1409
TMatrixTBase::GetColLwb
Int_t GetColLwb() const
Definition: TMatrixTBase.h:135
epsilon
REAL epsilon
Definition: triangle.c:617
TVectorT::TVectorT
TVectorT()
Definition: TVectorT.h:59
AreCompatible< Double_t, Float_t >
template Bool_t AreCompatible< Double_t, Float_t >(const TVectorD &v1, const TVectorF &v2, Int_t verbose)
TVectorT::SelectNonZeros
TVectorT< Element > & SelectNonZeros(const TVectorT< Element > &select)
Keep only element as selected through array select non-zero.
Definition: TVectorT.cxx:541
Long_t
long Long_t
Definition: RtypesCore.h:50
TBuffer::ReadClassBuffer
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
TMatrixTRow_const::GetPtr
const Element * GetPtr() const
Definition: TMatrixTUtils.h:135
TVectorT::GetLwb
Int_t GetLwb() const
Definition: TVectorT.h:79
init
static Int_t init()
Definition: RooClassFactory.cxx:53
Double_t
double Double_t
Definition: RtypesCore.h:55
ElementDiv
TVectorT< Element > & ElementDiv(TVectorT< Element > &target, const TVectorT< Element > &source)
Divide target by the source, element-by-element.
Definition: TVectorT.cxx:2072
TVectorT::Norm1
Element Norm1() const
Compute the 1-norm of the vector SUM{ |v[i]| }.
Definition: TVectorT.cxx:564
printf
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
VerifyVectorValue
Bool_t VerifyVectorValue(const TVectorT< Element > &v, Element val, Int_t verbose, Element maxDevAllow)
Validate that all elements of vector have value val within maxDevAllow .
Definition: TVectorT.cxx:2259
TMatrixTBase::GetNcols
Int_t GetNcols() const
Definition: TMatrixTBase.h:137
ULong_t
unsigned long ULong_t
Definition: RtypesCore.h:51
AddElemMult
TVectorT< Element > & AddElemMult(TVectorT< Element > &target, Element scalar, const TVectorT< Element > &source1, const TVectorT< Element > &source2)
Modify addition: target += scalar * ElementMult(source1,source2) .
Definition: TVectorT.cxx:1841
TVectorT::operator<=
Bool_t operator<=(Element val) const
Are all vector elements <= val?
Definition: TVectorT.cxx:1184
ElementMult< Double_t >
template TVectorD & ElementMult< Double_t >(TVectorD &target, const TVectorD &source)
TElementPosActionT::Operation
virtual void Operation(Element &element) const =0
Cppyy::Allocate
TCppObject_t Allocate(TCppType_t type)
Definition: Cppyy.cxx:228
TObject
Mother of all ROOT objects.
Definition: TObject.h:58
TMatrixTSparse::GetColIndexArray
virtual const Int_t * GetColIndexArray() const
Definition: TMatrixTSparse.h:212
TMatrixTColumn_const::GetPtr
const Element * GetPtr() const
Definition: TMatrixTUtils.h:234
TVectorT::Memcpy_m
Int_t Memcpy_m(Element *newp, const Element *oldp, Int_t copySize, Int_t newSize, Int_t oldSize)
Copy copySize doubles from *oldp to *newp .
Definition: TVectorT.cxx:121
TMatrixTSparseRow_const::GetMatrix
const TMatrixTBase< Element > * GetMatrix() const
Definition: TMatrixTUtils.h:599
TMatrixTDiag_const::GetPtr
const Element * GetPtr() const
Definition: TMatrixTUtils.h:331
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
TMatrixTColumn_const::GetInc
Int_t GetInc() const
Definition: TMatrixTUtils.h:233
TVectorT::Abs
TVectorT< Element > & Abs()
Take an absolute value of a vector, i.e. apply Abs() to each element.
Definition: TVectorT.cxx:461
TVectorT::MatchesNonZeroPattern
Bool_t MatchesNonZeroPattern(const TVectorT< Element > &select)
Check if vector elements as selected through array select are non-zero.
Definition: TVectorT.cxx:1235
TMatrixTSparseDiag_const::GetNdiags
Int_t GetNdiags() const
Definition: TMatrixTUtils.h:665
AreCompatible< Float_t, Float_t >
template Bool_t AreCompatible< Float_t, Float_t >(const TVectorF &v1, const TVectorF &v2, Int_t verbose)
Mult< Double_t, Double_t, Double_t >
template Double_t Mult< Double_t, Double_t, Double_t >(const TVectorD &v1, const TMatrixD &m, const TVectorD &v2)
TMatrixTRow_const
Definition: TMatrixDUtilsfwd.h:37
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
Compare
void Compare(const TVectorT< Element > &v1, const TVectorT< Element > &v2)
Compare two vectors and print out the result of the comparison.
Definition: TVectorT.cxx:2206
TMatrixTSparseDiag_const
Definition: TMatrixDUtilsfwd.h:43
Dot< Float_t >
template Float_t Dot< Float_t >(const TVectorF &v1, const TVectorF &v2)
AreCompatible
Bool_t AreCompatible(const TVectorT< Element1 > &v1, const TVectorT< Element2 > &v2, Int_t verbose)
Check if v1 and v2 are both valid and have the same shape.
Definition: TVectorT.cxx:2128
TMath::LocMin
Long64_t LocMin(Long64_t n, const T *a)
Definition: TMath.h:695
VerifyVectorValue< Float_t >
template Bool_t VerifyVectorValue< Float_t >(const TVectorF &m, Float_t val, Int_t verbose, Float_t maxDevAllow)
kTRUE
const Bool_t kTRUE
Definition: Rtypes.h:91
norm
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
TVectorT::New_m
Element * New_m(Int_t size)
default kTRUE, when Use array kFALSE
Definition: TVectorT.cxx:64
VerifyVectorValue< Double_t >
template Bool_t VerifyVectorValue< Double_t >(const TVectorD &m, Double_t val, Int_t verbose, Double_t maxDevAllow)
TVectorT::AddSomeConstant
void AddSomeConstant(Element val, const TVectorT< Element > &select)
Add to vector elements as selected through array select the value val.
Definition: TVectorT.cxx:1281
n
const Int_t n
Definition: legend1.C:16
TVectorT::operator<
Bool_t operator<(Element val) const
Are all vector elements < val?
Definition: TVectorT.cxx:1167
templateClassImp
templateClassImp(TVectorT) template< class Element > void TVectorT< Element >
Delete data pointer m, if it was assigned on the heap.
Definition: TVectorT.cxx:43
TMatrixTDiag_const::GetInc
Int_t GetInc() const
Definition: TMatrixTUtils.h:332
ElementMult< Float_t >
template TVectorF & ElementMult< Float_t >(TVectorF &target, const TVectorF &source)
TVectorT::SomePositive
Bool_t SomePositive(const TVectorT< Element > &select)
Check if vector elements as selected through array select are all positive.
Definition: TVectorT.cxx:1258
AddElemDiv< Double_t >
template TVectorD & AddElemDiv< Double_t >(TVectorD &target, Double_t scalar, const TVectorD &source1, const TVectorD &source2)
TBuffer::ReadVersion
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
TVectorT::operator+=
TVectorT< Element > & operator+=(Element val)
Add val to every element of the vector.
Definition: TVectorT.cxx:860
TVectorT::Apply
TVectorT< Element > & Apply(const TElementActionT< Element > &action)
Apply action to each element of the vector.
Definition: TVectorT.cxx:1319
TBuffer::ReadArray
virtual Int_t ReadArray(Bool_t *&b)=0
OuterProduct< Double_t, Double_t >
template TMatrixD OuterProduct< Double_t, Double_t >(const TVectorD &v1, const TVectorD &v2)
TMatrixTDiag_const::GetNdiags
Int_t GetNdiags() const
Definition: TMatrixTUtils.h:344
Drand
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
Definition: TMatrixTUtils.cxx:1852