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