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