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