Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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  Unless otherwise specified, vector indices always start with 0,
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 
44 templateClassImp(TVectorT);
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>
65 Element* TVectorT<Element>::New_m(Int_t size)
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>
82 void TVectorT<Element>::Add(const TVectorT<Element> &v)
83 {
84  if (gMatrixCheck && !AreCompatible(*this,v)) {
85  Error("Add(TVectorT<Element> &)","vector's not compatible");
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>
100 void TVectorT<Element>::Add(const TVectorT<Element> &v1,const TVectorT<Element> &v2)
101 {
102  if (gMatrixCheck) {
103  if ( !AreCompatible(*this,v1) && !AreCompatible(*this,v2)) {
104  Error("Add(TVectorT<Element> &)","vectors not compatible");
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>
149 void TVectorT<Element>::Allocate(Int_t nrows,Int_t row_lwb,Int_t init)
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>
175 TVectorT<Element>::TVectorT(Int_t n)
176 {
177  Allocate(n,0,1);
178 }
179 
180 ////////////////////////////////////////////////////////////////////////////////
181 /// Constructor [lwb..upb]-vector
182 
183 template<class Element>
184 TVectorT<Element>::TVectorT(Int_t lwb,Int_t upb)
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>
224 TVectorT<Element>::TVectorT(const TMatrixTRow_const<Element> &mr) : TObject()
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>
236 TVectorT<Element>::TVectorT(const TMatrixTColumn_const<Element> &mc) : TObject()
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>
248 TVectorT<Element>::TVectorT(const TMatrixTDiag_const<Element> &md) : TObject()
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>
263 TVectorT<Element>::TVectorT(Int_t lwb,Int_t upb,Double_t iv1, ...)
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>
292 TVectorT<Element> &TVectorT<Element>::ResizeTo(Int_t lwb,Int_t upb)
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>
347 TVectorT<Element> &TVectorT<Element>::Use(Int_t lwb,Int_t upb,Element *data)
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>
371 TVectorT<Element> &TVectorT<Element>::GetSub(Int_t row_lwb,Int_t row_upb,TVectorT<Element> &target,Option_t *option) const
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>
420 TVectorT<Element> &TVectorT<Element>::SetSub(Int_t row_lwb,const TVectorT<Element> &source)
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>
451 TVectorT<Element> &TVectorT<Element>::Zero()
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>
462 TVectorT<Element> &TVectorT<Element>::Abs()
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>
480 TVectorT<Element> &TVectorT<Element>::Sqr()
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>
498 TVectorT<Element> &TVectorT<Element>::Sqrt()
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>
520 TVectorT<Element> &TVectorT<Element>::Invert()
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>
542 TVectorT<Element> &TVectorT<Element>::SelectNonZeros(const TVectorT<Element> &select)
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>
565 Element TVectorT<Element>::Norm1() const
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>
582 Element TVectorT<Element>::Norm2Sqr() const
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>
601 Element TVectorT<Element>::NormInf() const
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>
618 Int_t TVectorT<Element>::NonZeros() const
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>
678 TVectorT<Element> &TVectorT<Element>::operator=(const TVectorT<Element> &source)
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>
696 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTRow_const<Element> &mr)
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>
727 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTColumn_const<Element> &mc)
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>
758 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTDiag_const<Element> &md)
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>
790 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTSparseRow_const<Element> &mr)
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>
821 TVectorT<Element> &TVectorT<Element>::operator=(const TMatrixTSparseDiag_const<Element> &md)
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>
845 TVectorT<Element> &TVectorT<Element>::operator=(Element val)
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>
861 TVectorT<Element> &TVectorT<Element>::operator+=(Element val)
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>
877 TVectorT<Element> &TVectorT<Element>::operator-=(Element val)
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>
893 TVectorT<Element> &TVectorT<Element>::operator*=(Element val)
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 ////////////////////////////////////////////////////////////////////////////////
906 /// Add vector source
907 
908 template<class Element>
909 TVectorT<Element> &TVectorT<Element>::operator+=(const TVectorT<Element> &source)
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>
929 TVectorT<Element> &TVectorT<Element>::operator-=(const TVectorT<Element> &source)
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>
950 TVectorT<Element> &TVectorT<Element>::operator*=(const TMatrixT<Element> &a)
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>
1017 TVectorT<Element> &TVectorT<Element>::operator*=(const TMatrixTSparse<Element> &a)
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>
1080 TVectorT<Element> &TVectorT<Element>::operator*=(const TMatrixTSym<Element> &a)
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>
1134 Bool_t TVectorT<Element>::operator==(Element val) const
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>
1151 Bool_t TVectorT<Element>::operator!=(Element val) const
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>
1168 Bool_t TVectorT<Element>::operator<(Element val) const
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>
1185 Bool_t TVectorT<Element>::operator<=(Element val) const
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>
1202 Bool_t TVectorT<Element>::operator>(Element val) const
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>
1219 Bool_t TVectorT<Element>::operator>=(Element val) const
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>
1236 Bool_t TVectorT<Element>::MatchesNonZeroPattern(const TVectorT<Element> &select)
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>
1259 Bool_t TVectorT<Element>::SomePositive(const TVectorT<Element> &select)
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>
1282 void TVectorT<Element>::AddSomeConstant(Element val,const TVectorT<Element> &select)
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>
1320 TVectorT<Element> &TVectorT<Element>::Apply(const TElementActionT<Element> &action)
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>
1333 TVectorT<Element> &TVectorT<Element>::Apply(const TElementPosActionT<Element> &action)
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>
1351 void TVectorT<Element>::Draw(Option_t *option)
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>
1361 void TVectorT<Element>::Print(Option_t *flag) const
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>
1384 Bool_t operator==(const TVectorT<Element> &v1,const TVectorT<Element> &v2)
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>
1410 TVectorT<Element> operator+(const TVectorT<Element> &source1,const TVectorT<Element> &source2)
1411 {
1412  TVectorT<Element> target = source1;
1413  target += source2;
1414  return target;
1415 }
1416 
1417 ////////////////////////////////////////////////////////////////////////////////
1418 /// Return source1-source2
1419 
1420 template<class Element>
1421 TVectorT<Element> operator-(const TVectorT<Element> &source1,const TVectorT<Element> &source2)
1422 {
1423  TVectorT<Element> target = source1;
1424  target -= source2;
1425  return target;
1426 }
1427 
1428 ////////////////////////////////////////////////////////////////////////////////
1429 /// return A * source
1430 
1431 template<class Element>
1432 TVectorT<Element> operator*(const TMatrixT<Element> &a,const TVectorT<Element> &source)
1433 {
1434  R__ASSERT(a.IsValid());
1435  TVectorT<Element> target(a.GetRowLwb(),a.GetRowUpb());
1436  return Add(target,Element(1.0),a,source);
1437 }
1438 
1439 ////////////////////////////////////////////////////////////////////////////////
1440 /// return A * source
1441 
1442 template<class Element>
1443 TVectorT<Element> operator*(const TMatrixTSym<Element> &a,const TVectorT<Element> &source)
1444 {
1445  R__ASSERT(a.IsValid());
1446  TVectorT<Element> target(a.GetRowLwb(),a.GetRowUpb());
1447  return Add(target,Element(1.0),a,source);
1448 }
1449 
1450 ////////////////////////////////////////////////////////////////////////////////
1451 /// return A * source
1452 
1453 template<class Element>
1454 TVectorT<Element> operator*(const TMatrixTSparse<Element> &a,const TVectorT<Element> &source)
1455 {
1456  R__ASSERT(a.IsValid());
1457  TVectorT<Element> target(a.GetRowLwb(),a.GetRowUpb());
1458  return Add(target,Element(1.0),a,source);
1459 }
1460 
1461 ////////////////////////////////////////////////////////////////////////////////
1462 /// return val * source
1463 
1464 template<class Element>
1465 TVectorT<Element> operator*(Element val,const TVectorT<Element> &source)
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>
1492 TMatrixT<Element1>
1493 OuterProduct(const TVectorT<Element1> &v1,const TVectorT<Element2> &v2)
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>
1508 TMatrixT<Element1>
1509 &OuterProduct(TMatrixT<Element1> &target,const TVectorT<Element2> &v1,const TVectorT<Element3> &v2)
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>
1614 TVectorT<Element> &Add(TVectorT<Element> &target,Element scalar,
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>
1692 TVectorT<Element> &Add(TVectorT<Element> &target,Element scalar,
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>
1764 TVectorT<Element> &Add(TVectorT<Element> &target,Element scalar,
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>
1842 TVectorT<Element> &AddElemMult(TVectorT<Element> &target,Element scalar,
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>
1875 TVectorT<Element> &AddElemMult(TVectorT<Element> &target,Element scalar,
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>
1915 TVectorT<Element> &AddElemDiv(TVectorT<Element> &target,Element scalar,
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();
1935  Error("AddElemDiv","source2 (%d) is zero",irow);
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();
1945  Error("AddElemDiv","source2 (%d) is zero",irow);
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();
1955  Error("AddElemDiv","source2 (%d) is zero",irow);
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>
1969 TVectorT<Element> &AddElemDiv(TVectorT<Element> &target,Element scalar,
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();
1992  Error("AddElemDiv","source2 (%d) is zero",irow);
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();
2004  Error("AddElemDiv","source2 (%d) is zero",irow);
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();
2016  Error("AddElemDiv","source2 (%d) is zero",irow);
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>
2030 TVectorT<Element> &ElementMult(TVectorT<Element> &target,const TVectorT<Element> &source)
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>
2050 TVectorT<Element> &ElementMult(TVectorT<Element> &target,const TVectorT<Element> &source,const TVectorT<Element> &select)
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>
2073 TVectorT<Element> &ElementDiv(TVectorT<Element> &target,const TVectorT<Element> &source)
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>
2099 TVectorT<Element> &ElementDiv(TVectorT<Element> &target,const TVectorT<Element> &source,const TVectorT<Element> &select)
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>
2129 Bool_t AreCompatible(const TVectorT<Element1> &v1,const TVectorT<Element2> &v2,Int_t verbose)
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>
2155 Bool_t AreCompatible(const TMatrixT<Element1> &m,const TVectorT<Element2> &v,Int_t verbose)
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>
2181 Bool_t AreCompatible(const TVectorT<Element1> &v,const TMatrixT<Element2> &m,Int_t verbose)
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>
2207 void Compare(const TVectorT<Element> &v1,const TVectorT<Element> &v2)
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>
2260 Bool_t VerifyVectorValue(const TVectorT<Element> &v,Element val,
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>
2295 Bool_t VerifyVectorIdentity(const TVectorT<Element> &v1,const TVectorT<Element> &v2,
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>
2334 void TVectorT<Element>::Streamer(TBuffer &R__b)
2335 {
2336  if (R__b.IsReading()) {
2337  UInt_t R__s, R__c;
2338  Version_t R__v = R__b.ReadVersion(&R__s,&R__c);
2339  if (R__v > 1) {
2340  Clear();
2341  R__b.ReadClassBuffer(TVectorT<Element>::Class(),this,R__v,R__s,R__c);
2342  } else { //====process old versions before automatic schema evolution
2343  TObject::Streamer(R__b);
2344  R__b >> fRowLwb;
2345  fNrows = R__b.ReadArray(fElements);
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;
2351  fElements = fDataStack;
2352  } else if (fNrows < 0)
2353  Invalidate();
2354 
2355  if (R__v < 3)
2356  MakeValid();
2357  } else {
2358  R__b.WriteClassBuffer(TVectorT<Element>::Class(),this);
2359  }
2360 }
2361 
2362 #include "TMatrixFfwd.h"
2363 #include "TMatrixFSymfwd.h"
2364 #include "TMatrixFSparsefwd.h"
2365 
2366 template class TVectorT<Float_t>;
2367 
2368 template Bool_t operator== <Float_t>(const TVectorF &source1,const TVectorF &source2);
2369 template TVectorF operator+ <Float_t>(const TVectorF &source1,const TVectorF &source2);
2370 template TVectorF operator- <Float_t>(const TVectorF &source1,const TVectorF &source2);
2371 template Float_t operator* <Float_t>(const TVectorF &source1,const TVectorF &source2);
2372 template TVectorF operator* <Float_t>(const TMatrixF &a, const TVectorF &source);
2373 template TVectorF operator* <Float_t>(const TMatrixFSym &a, const TVectorF &source);
2374 template TVectorF operator* <Float_t>(const TMatrixFSparse &a, const TVectorF &source);
2375 template TVectorF operator* <Float_t>( Float_t val, const TVectorF &source);
2376 
2377 template Float_t Dot <Float_t>(const TVectorF &v1, const TVectorF &v2);
2378 template TMatrixF OuterProduct <Float_t,Float_t>
2379  (const TVectorF &v1, const TVectorF &v2);
2380 template TMatrixF &OuterProduct <Float_t,Float_t,Float_t>
2381  ( TMatrixF &target, const TVectorF &v1, const TVectorF &v2);
2382 template Float_t Mult <Float_t,Float_t,Float_t>
2383  (const TVectorF &v1, const TMatrixF &m, const TVectorF &v2);
2384 
2385 template TVectorF &Add <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source);
2386 template TVectorF &Add <Float_t>( TVectorF &target, Float_t scalar, const TMatrixF &a,
2387  const TVectorF &source);
2388 template TVectorF &Add <Float_t>( TVectorF &target, Float_t scalar, const TMatrixFSym &a,
2389  const TVectorF &source);
2390 template TVectorF &Add <Float_t>( TVectorF &target, Float_t scalar, const TMatrixFSparse &a,
2391  const TVectorF &source);
2392 template TVectorF &AddElemMult <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source1,
2393  const TVectorF &source2);
2394 template TVectorF &AddElemMult <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source1,
2395  const TVectorF &source2,const TVectorF &select);
2396 template TVectorF &AddElemDiv <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source1,
2397  const TVectorF &source2);
2398 template TVectorF &AddElemDiv <Float_t>( TVectorF &target, Float_t scalar, const TVectorF &source1,
2399  const TVectorF &source2,const TVectorF &select);
2400 template TVectorF &ElementMult <Float_t>( TVectorF &target, const TVectorF &source);
2401 template TVectorF &ElementMult <Float_t>( TVectorF &target, const TVectorF &source, const TVectorF &select);
2402 template TVectorF &ElementDiv <Float_t>( TVectorF &target, const TVectorF &source);
2403 template TVectorF &ElementDiv <Float_t>( TVectorF &target, const TVectorF &source, const TVectorF &select);
2404 
2405 template Bool_t AreCompatible <Float_t,Float_t> (const TVectorF &v1,const TVectorF &v2,Int_t verbose);
2406 template Bool_t AreCompatible <Float_t,Double_t>(const TVectorF &v1,const TVectorD &v2,Int_t verbose);
2407 template Bool_t AreCompatible <Float_t,Float_t> (const TMatrixF &m, const TVectorF &v, Int_t verbose);
2408 template Bool_t AreCompatible <Float_t,Float_t> (const TVectorF &v, const TMatrixF &m, Int_t verbose);
2409 
2410 template void Compare <Float_t> (const TVectorF &v1,const TVectorF &v2);
2411 template Bool_t VerifyVectorValue <Float_t> (const TVectorF &m, Float_t val,Int_t verbose,Float_t maxDevAllow);
2412 template Bool_t VerifyVectorIdentity<Float_t> (const TVectorF &m1,const TVectorF &m2, Int_t verbose,Float_t maxDevAllow);
2413 
2414 #include "TMatrixDfwd.h"
2415 #include "TMatrixDSymfwd.h"
2416 #include "TMatrixDSparsefwd.h"
2417 
2418 template class TVectorT<Double_t>;
2419 
2420 template Bool_t operator== <Double_t>(const TVectorD &source1,const TVectorD &source2);
2421 template TVectorD operator+ <Double_t>(const TVectorD &source1,const TVectorD &source2);
2422 template TVectorD operator- <Double_t>(const TVectorD &source1,const TVectorD &source2);
2423 template Double_t operator* <Double_t>(const TVectorD &source1,const TVectorD &source2);
2424 template TVectorD operator* <Double_t>(const TMatrixD &a, const TVectorD &source);
2425 template TVectorD operator* <Double_t>(const TMatrixDSym &a, const TVectorD &source);
2426 template TVectorD operator* <Double_t>(const TMatrixDSparse &a, const TVectorD &source);
2427 template TVectorD operator* <Double_t>( Double_t val, const TVectorD &source);
2428 
2429 template Double_t Dot <Double_t>(const TVectorD &v1, const TVectorD &v2);
2430 template TMatrixD OuterProduct <Double_t,Double_t>
2431  (const TVectorD &v1, const TVectorD &v2);
2432 template TMatrixD &OuterProduct <Double_t,Double_t,Double_t>
2433  ( TMatrixD &target, const TVectorD &v1, const TVectorD &v2);
2434 template Double_t Mult <Double_t,Double_t,Double_t>
2435  (const TVectorD &v1, const TMatrixD &m, const TVectorD &v2);
2436 
2437 template TVectorD &Add <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source);
2438 template TVectorD &Add <Double_t>( TVectorD &target, Double_t scalar, const TMatrixD &a,
2439  const TVectorD &source);
2440 template TVectorD &Add <Double_t>( TVectorD &target, Double_t scalar, const TMatrixDSym &a
2441  , const TVectorD &source);
2442 template TVectorD &Add <Double_t>( TVectorD &target, Double_t scalar, const TMatrixDSparse &a
2443  , const TVectorD &source);
2444 template TVectorD &AddElemMult <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source1,
2445  const TVectorD &source2);
2446 template TVectorD &AddElemMult <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source1,
2447  const TVectorD &source2,const TVectorD &select);
2448 template TVectorD &AddElemDiv <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source1,
2449  const TVectorD &source2);
2450 template TVectorD &AddElemDiv <Double_t>( TVectorD &target, Double_t scalar, const TVectorD &source1,
2451  const TVectorD &source2,const TVectorD &select);
2452 template TVectorD &ElementMult <Double_t>( TVectorD &target, const TVectorD &source);
2453 template TVectorD &ElementMult <Double_t>( TVectorD &target, const TVectorD &source, const TVectorD &select);
2454 template TVectorD &ElementDiv <Double_t>( TVectorD &target, const TVectorD &source);
2455 template TVectorD &ElementDiv <Double_t>( TVectorD &target, const TVectorD &source, const TVectorD &select);
2456 
2457 template Bool_t AreCompatible <Double_t,Double_t>(const TVectorD &v1,const TVectorD &v2,Int_t verbose);
2458 template Bool_t AreCompatible <Double_t,Float_t> (const TVectorD &v1,const TVectorF &v2,Int_t verbose);
2459 template Bool_t AreCompatible <Double_t,Double_t>(const TMatrixD &m, const TVectorD &v, Int_t verbose);
2460 template Bool_t AreCompatible <Double_t,Double_t>(const TVectorD &v, const TMatrixD &m, Int_t verbose);
2461 
2462 template void Compare <Double_t> (const TVectorD &v1,const TVectorD &v2);
2463 template Bool_t VerifyVectorValue <Double_t> (const TVectorD &m, Double_t val,Int_t verbose,Double_t maxDevAllow);
2464 template Bool_t VerifyVectorIdentity<Double_t> (const TVectorD &m1,const TVectorD &m2, Int_t verbose,Double_t maxDevAllow);