Logo ROOT  

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