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