Logo ROOT  

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