Logo ROOT   6.08/07
Reference Guide
TMatrixTUtils.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 TMatrixTUtils
13  \ingroup Matrix
14 
15  Matrix utility classes.
16 
17  Templates of utility classes in the Linear Algebra Package.
18  The following classes are defined here:
19 
20  Different matrix views without copying data elements :
21 ~~~
22  TMatrixTRow_const TMatrixTRow
23  TMatrixTColumn_const TMatrixTColumn
24  TMatrixTDiag_const TMatrixTDiag
25  TMatrixTFlat_const TMatrixTFlat
26  TMatrixTSub_const TMatrixTSub
27  TMatrixTSparseRow_const TMatrixTSparseRow
28  TMatrixTSparseDiag_const TMatrixTSparseDiag
29 
30  TElementActionT
31  TElementPosActionT
32 ~~~
33 */
34 
35 #include "TMatrixTUtils.h"
36 #include "TMatrixT.h"
37 #include "TMatrixTSym.h"
38 #include "TMatrixTSparse.h"
39 #include "TMath.h"
40 #include "TVectorT.h"
41 
42 ////////////////////////////////////////////////////////////////////////////////
43 /// Constructor with row "row" of matrix
44 
45 template<class Element>
47 {
48  R__ASSERT(matrix.IsValid());
49 
50  fRowInd = row-matrix.GetRowLwb();
51  if (fRowInd >= matrix.GetNrows() || fRowInd < 0) {
52  Error("TMatrixTRow_const(const TMatrixT<Element> &,Int_t)","row index out of bounds");
53  fMatrix = 0;
54  fPtr = 0;
55  fInc = 0;
56  return;
57  }
58 
59  fMatrix = &matrix;
60  fPtr = matrix.GetMatrixArray()+fRowInd*matrix.GetNcols();
61  fInc = 1;
62 }
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 /// Constructor with row "row" of symmetric matrix
66 
67 template<class Element>
69 {
70  R__ASSERT(matrix.IsValid());
71 
72  fRowInd = row-matrix.GetRowLwb();
73  if (fRowInd >= matrix.GetNrows() || fRowInd < 0) {
74  Error("TMatrixTRow_const(const TMatrixTSym &,Int_t)","row index out of bounds");
75  fMatrix = 0;
76  fPtr = 0;
77  fInc = 0;
78  return;
79  }
80 
81  fMatrix = &matrix;
82  fPtr = matrix.GetMatrixArray()+fRowInd*matrix.GetNcols();
83  fInc = 1;
84 }
85 
86 ////////////////////////////////////////////////////////////////////////////////
87 /// Constructor with row "row" of symmetric matrix
88 
89 template<class Element>
91  :TMatrixTRow_const<Element>(matrix,row)
92 {
93 }
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Constructor with row "row" of symmetric matrix
97 
98 template<class Element>
100  :TMatrixTRow_const<Element>(matrix,row)
101 {
102 }
103 
104 ////////////////////////////////////////////////////////////////////////////////
105 /// Copy constructor
106 
107 template<class Element>
109 {
110  *this = mr;
111 }
112 
113 ////////////////////////////////////////////////////////////////////////////////
114 /// Assign val to every element of the matrix row.
115 
116 template<class Element>
118 {
119  R__ASSERT(this->fMatrix->IsValid());
120  Element *rp = const_cast<Element *>(this->fPtr);
121  for ( ; rp < this->fPtr+this->fMatrix->GetNcols(); rp += this->fInc)
122  *rp = val;
123 }
124 
125 ////////////////////////////////////////////////////////////////////////////////
126 /// Add val to every element of the matrix row.
127 
128 template<class Element>
130 {
131  R__ASSERT(this->fMatrix->IsValid());
132  Element *rp = const_cast<Element *>(this->fPtr);
133  for ( ; rp < this->fPtr+this->fMatrix->GetNcols(); rp += this->fInc)
134  *rp += val;
135 }
136 
137 ////////////////////////////////////////////////////////////////////////////////
138 /// Multiply every element of the matrix row with val.
139 
140 template<class Element>
142 {
143  R__ASSERT(this->fMatrix->IsValid());
144  Element *rp = const_cast<Element *>(this->fPtr);
145  for ( ; rp < this->fPtr + this->fMatrix->GetNcols(); rp += this->fInc)
146  *rp *= val;
147 }
148 
149 ////////////////////////////////////////////////////////////////////////////////
150 /// Assignment operator
151 
152 template<class Element>
154 {
155  const TMatrixTBase<Element> *mt = mr.GetMatrix();
156  if (this->fMatrix->GetMatrixArray() == mt->GetMatrixArray() && this->fRowInd == mr.GetRowIndex()) return;
157 
158  R__ASSERT(this->fMatrix->IsValid());
159  R__ASSERT(mt->IsValid());
160 
161  if (this->fMatrix->GetNcols() != mt->GetNcols() || this->fMatrix->GetColLwb() != mt->GetColLwb()) {
162  Error("operator=(const TMatrixTRow_const &)", "matrix rows not compatible");
163  return;
164  }
165 
166  Element *rp1 = const_cast<Element *>(this->fPtr);
167  const Element *rp2 = mr.GetPtr();
168  for ( ; rp1 < this->fPtr+this->fMatrix->GetNcols(); rp1 += this->fInc,rp2 += mr.GetInc())
169  *rp1 = *rp2;
170 }
171 
172 ////////////////////////////////////////////////////////////////////////////////
173 /// Assign a vector to a matrix row. The vector is considered row-vector
174 /// to allow the assignment in the strict sense.
175 
176 template<class Element>
178 {
179  R__ASSERT(this->fMatrix->IsValid());
180  R__ASSERT(vec.IsValid());
181 
182  if (this->fMatrix->GetColLwb() != vec.GetLwb() || this->fMatrix->GetNcols() != vec.GetNrows()) {
183  Error("operator=(const TVectorT &)","vector length != matrix-row length");
184  return;
185  }
186 
187  Element *rp = const_cast<Element *>(this->fPtr);
188  const Element *vp = vec.GetMatrixArray();
189  for ( ; rp < this->fPtr+this->fMatrix->GetNcols(); rp += this->fInc)
190  *rp = *vp++;
191 }
192 
193 ////////////////////////////////////////////////////////////////////////////////
194 /// Add to every element of the matrix row the corresponding element of row r.
195 
196 template<class Element>
198 {
199  const TMatrixTBase<Element> *mt = r.GetMatrix();
200 
201  R__ASSERT(this->fMatrix->IsValid());
202  R__ASSERT(mt->IsValid());
203 
204  if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
205  Error("operator+=(const TMatrixTRow_const &)","different row lengths");
206  return;
207  }
208 
209  Element *rp1 = const_cast<Element *>(this->fPtr);
210  const Element *rp2 = r.GetPtr();
211  for ( ; rp1 < this->fPtr+this->fMatrix->GetNcols(); rp1 += this->fInc,rp2 += r.GetInc())
212  *rp1 += *rp2;
213 }
214 
215 ////////////////////////////////////////////////////////////////////////////////
216 /// Multiply every element of the matrix row with the
217 /// corresponding element of row r.
218 
219 template<class Element>
221 {
222  const TMatrixTBase<Element> *mt = r.GetMatrix();
223 
224  R__ASSERT(this->fMatrix->IsValid());
225  R__ASSERT(mt->IsValid());
226 
227  if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
228  Error("operator*=(const TMatrixTRow_const &)","different row lengths");
229  return;
230  }
231 
232  Element *rp1 = const_cast<Element *>(this->fPtr);
233  const Element *rp2 = r.GetPtr();
234  for ( ; rp1 < this->fPtr+this->fMatrix->GetNcols(); rp1 += this->fInc,rp2 += r.GetInc())
235  *rp1 *= *rp2;
236 }
237 
238 ////////////////////////////////////////////////////////////////////////////////
239 /// Constructor with column "col" of matrix
240 
241 template<class Element>
243 {
244  R__ASSERT(matrix.IsValid());
245 
246  this->fColInd = col-matrix.GetColLwb();
247  if (this->fColInd >= matrix.GetNcols() || this->fColInd < 0) {
248  Error("TMatrixTColumn_const(const TMatrixT &,Int_t)","column index out of bounds");
249  fMatrix = 0;
250  fPtr = 0;
251  fInc = 0;
252  return;
253  }
254 
255  fMatrix = &matrix;
256  fPtr = matrix.GetMatrixArray()+fColInd;
257  fInc = matrix.GetNcols();
258 }
259 
260 ////////////////////////////////////////////////////////////////////////////////
261 /// Constructor with column "col" of matrix
262 
263 template<class Element>
265 {
266  R__ASSERT(matrix.IsValid());
267 
268  fColInd = col-matrix.GetColLwb();
269  if (fColInd >= matrix.GetNcols() || fColInd < 0) {
270  Error("TMatrixTColumn_const(const TMatrixTSym &,Int_t)","column index out of bounds");
271  fMatrix = 0;
272  fPtr = 0;
273  fInc = 0;
274  return;
275  }
276 
277  fMatrix = &matrix;
278  fPtr = matrix.GetMatrixArray()+fColInd;
279  fInc = matrix.GetNcols();
280 }
281 
282 ////////////////////////////////////////////////////////////////////////////////
283 /// Constructor with column "col" of matrix
284 
285 template<class Element>
287  :TMatrixTColumn_const<Element>(matrix,col)
288 {
289 }
290 
291 ////////////////////////////////////////////////////////////////////////////////
292 /// Constructor with column "col" of matrix
293 
294 template<class Element>
296  :TMatrixTColumn_const<Element>(matrix,col)
297 {
298 }
299 
300 ////////////////////////////////////////////////////////////////////////////////
301 /// Copy constructor
302 
303 template<class Element>
305 {
306  *this = mc;
307 }
308 
309 ////////////////////////////////////////////////////////////////////////////////
310 /// Assign val to every element of the matrix column.
311 
312 template<class Element>
314 {
315  R__ASSERT(this->fMatrix->IsValid());
316  Element *cp = const_cast<Element *>(this->fPtr);
317  for ( ; cp < this->fPtr+this->fMatrix->GetNoElements(); cp += this->fInc)
318  *cp = val;
319 }
320 
321 ////////////////////////////////////////////////////////////////////////////////
322 /// Add val to every element of the matrix column.
323 
324 template<class Element>
326 {
327  R__ASSERT(this->fMatrix->IsValid());
328  Element *cp = const_cast<Element *>(this->fPtr);
329  for ( ; cp < this->fPtr+this->fMatrix->GetNoElements(); cp += this->fInc)
330  *cp += val;
331 }
332 
333 ////////////////////////////////////////////////////////////////////////////////
334 /// Multiply every element of the matrix column with val.
335 
336 template<class Element>
338 {
339  R__ASSERT(this->fMatrix->IsValid());
340  Element *cp = const_cast<Element *>(this->fPtr);
341  for ( ; cp < this->fPtr+this->fMatrix->GetNoElements(); cp += this->fInc)
342  *cp *= val;
343 }
344 
345 ////////////////////////////////////////////////////////////////////////////////
346 /// Assignment operator
347 
348 template<class Element>
350 {
351  const TMatrixTBase<Element> *mt = mc.GetMatrix();
352  if (this->fMatrix->GetMatrixArray() == mt->GetMatrixArray() && this->fColInd == mc.GetColIndex()) return;
353 
354  R__ASSERT(this->fMatrix->IsValid());
355  R__ASSERT(mt->IsValid());
356 
357  if (this->fMatrix->GetNrows() != mt->GetNrows() || this->fMatrix->GetRowLwb() != mt->GetRowLwb()) {
358  Error("operator=(const TMatrixTColumn_const &)", "matrix columns not compatible");
359  return;
360  }
361 
362  Element *cp1 = const_cast<Element *>(this->fPtr);
363  const Element *cp2 = mc.GetPtr();
364  for ( ; cp1 < this->fPtr+this->fMatrix->GetNoElements(); cp1 += this->fInc,cp2 += mc.GetInc())
365  *cp1 = *cp2;
366 }
367 
368 ////////////////////////////////////////////////////////////////////////////////
369 /// Assign a vector to a matrix column.
370 
371 template<class Element>
373 {
374  R__ASSERT(this->fMatrix->IsValid());
375  R__ASSERT(vec.IsValid());
376 
377  if (this->fMatrix->GetRowLwb() != vec.GetLwb() || this->fMatrix->GetNrows() != vec.GetNrows()) {
378  Error("operator=(const TVectorT &)","vector length != matrix-column length");
379  return;
380  }
381 
382  Element *cp = const_cast<Element *>(this->fPtr);
383  const Element *vp = vec.GetMatrixArray();
384  for ( ; cp < this->fPtr+this->fMatrix->GetNoElements(); cp += this->fInc)
385  *cp = *vp++;
386 
387  R__ASSERT(vp == vec.GetMatrixArray()+vec.GetNrows());
388 }
389 
390 ////////////////////////////////////////////////////////////////////////////////
391 /// Add to every element of the matrix row the corresponding element of row mc.
392 
393 template<class Element>
395 {
396  const TMatrixTBase<Element> *mt = mc.GetMatrix();
397 
398  R__ASSERT(this->fMatrix->IsValid());
399  R__ASSERT(mt->IsValid());
400 
401  if (this->fMatrix->GetRowLwb() != mt->GetRowLwb() || this->fMatrix->GetNrows() != mt->GetNrows()) {
402  Error("operator+=(const TMatrixTColumn_const &)","different row lengths");
403  return;
404  }
405 
406  Element *cp1 = const_cast<Element *>(this->fPtr);
407  const Element *cp2 = mc.GetPtr();
408  for ( ; cp1 < this->fPtr+this->fMatrix->GetNoElements(); cp1 += this->fInc,cp2 += mc.GetInc())
409  *cp1 += *cp2;
410 }
411 
412 ////////////////////////////////////////////////////////////////////////////////
413 /// Multiply every element of the matrix column with the
414 /// corresponding element of column mc.
415 
416 template<class Element>
418 {
419  const TMatrixTBase<Element> *mt = mc.GetMatrix();
420 
421  R__ASSERT(this->fMatrix->IsValid());
422  R__ASSERT(mt->IsValid());
423 
424  if (this->fMatrix->GetRowLwb() != mt->GetRowLwb() || this->fMatrix->GetNrows() != mt->GetNrows()) {
425  Error("operator*=(const TMatrixTColumn_const &)","different row lengths");
426  return;
427  }
428 
429  Element *cp1 = const_cast<Element *>(this->fPtr);
430  const Element *cp2 = mc.GetPtr();
431  for ( ; cp1 < this->fPtr+this->fMatrix->GetNoElements(); cp1 += this->fInc,cp2 += mc.GetInc())
432  *cp1 *= *cp2;
433 }
434 
435 ////////////////////////////////////////////////////////////////////////////////
436 /// Constructor
437 
438 template<class Element>
440 {
441  R__ASSERT(matrix.IsValid());
442 
443  fMatrix = &matrix;
444  fNdiag = TMath::Min(matrix.GetNrows(),matrix.GetNcols());
445  fPtr = matrix.GetMatrixArray();
446  fInc = matrix.GetNcols()+1;
447 }
448 
449 ////////////////////////////////////////////////////////////////////////////////
450 /// Constructor
451 
452 template<class Element>
454 {
455  R__ASSERT(matrix.IsValid());
456 
457  fMatrix = &matrix;
458  fNdiag = TMath::Min(matrix.GetNrows(),matrix.GetNcols());
459  fPtr = matrix.GetMatrixArray();
460  fInc = matrix.GetNcols()+1;
461 }
462 
463 ////////////////////////////////////////////////////////////////////////////////
464 /// Constructor
465 
466 template<class Element>
468  :TMatrixTDiag_const<Element>(matrix)
469 {
470 }
471 
472 ////////////////////////////////////////////////////////////////////////////////
473 /// Constructor
474 
475 template<class Element>
477  :TMatrixTDiag_const<Element>(matrix)
478 {
479 }
480 
481 ////////////////////////////////////////////////////////////////////////////////
482 /// Copy constructor
483 
484 template<class Element>
486 {
487  *this = md;
488 }
489 
490 ////////////////////////////////////////////////////////////////////////////////
491 /// Assign val to every element of the matrix diagonal.
492 
493 template<class Element>
495 {
496  R__ASSERT(this->fMatrix->IsValid());
497  Element *dp = const_cast<Element *>(this->fPtr);
498  for (Int_t i = 0; i < this->fNdiag; i++, dp += this->fInc)
499  *dp = val;
500 }
501 
502 ////////////////////////////////////////////////////////////////////////////////
503 /// Assign val to every element of the matrix diagonal.
504 
505 template<class Element>
507 {
508  R__ASSERT(this->fMatrix->IsValid());
509  Element *dp = const_cast<Element *>(this->fPtr);
510  for (Int_t i = 0; i < this->fNdiag; i++, dp += this->fInc)
511  *dp += val;
512 }
513 
514 ////////////////////////////////////////////////////////////////////////////////
515 /// Assign val to every element of the matrix diagonal.
516 
517 template<class Element>
519 {
520  R__ASSERT(this->fMatrix->IsValid());
521  Element *dp = const_cast<Element *>(this->fPtr);
522  for (Int_t i = 0; i < this->fNdiag; i++, dp += this->fInc)
523  *dp *= val;
524 }
525 
526 ////////////////////////////////////////////////////////////////////////////////
527 /// Assignment operator
528 
529 template<class Element>
531 {
532  const TMatrixTBase<Element> *mt = md.GetMatrix();
533  if (this->fMatrix == mt) return;
534 
535  R__ASSERT(this->fMatrix->IsValid());
536  R__ASSERT(mt->IsValid());
537 
538  if (this->GetNdiags() != md.GetNdiags()) {
539  Error("operator=(const TMatrixTDiag_const &)","diagonals not compatible");
540  return;
541  }
542 
543  Element *dp1 = const_cast<Element *>(this->fPtr);
544  const Element *dp2 = md.GetPtr();
545  for (Int_t i = 0; i < this->fNdiag; i++, dp1 += this->fInc, dp2 += md.GetInc())
546  *dp1 = *dp2;
547 }
548 
549 ////////////////////////////////////////////////////////////////////////////////
550 /// Assign a vector to the matrix diagonal.
551 
552 template<class Element>
554 {
555  R__ASSERT(this->fMatrix->IsValid());
556  R__ASSERT(vec.IsValid());
557 
558  if (this->fNdiag != vec.GetNrows()) {
559  Error("operator=(const TVectorT &)","vector length != matrix-diagonal length");
560  return;
561  }
562 
563  Element *dp = const_cast<Element *>(this->fPtr);
564  const Element *vp = vec.GetMatrixArray();
565  for ( ; vp < vec.GetMatrixArray()+vec.GetNrows(); dp += this->fInc)
566  *dp = *vp++;
567 }
568 
569 ////////////////////////////////////////////////////////////////////////////////
570 /// Add to every element of the matrix diagonal the
571 /// corresponding element of diagonal md.
572 
573 template<class Element>
575 {
576  const TMatrixTBase<Element> *mt = md.GetMatrix();
577 
578  R__ASSERT(this->fMatrix->IsValid());
579  R__ASSERT(mt->IsValid());
580  if (this->fNdiag != md.GetNdiags()) {
581  Error("operator=(const TMatrixTDiag_const &)","matrix-diagonal's different length");
582  return;
583  }
584 
585  Element *dp1 = const_cast<Element *>(this->fPtr);
586  const Element *dp2 = md.GetPtr();
587  for (Int_t i = 0; i < this->fNdiag; i++, dp1 += this->fInc, dp2 += md.GetInc())
588  *dp1 += *dp2;
589 }
590 
591 ////////////////////////////////////////////////////////////////////////////////
592 /// Multiply every element of the matrix diagonal with the
593 /// corresponding element of diagonal md.
594 
595 template<class Element>
597 {
598  const TMatrixTBase<Element> *mt = md.GetMatrix();
599 
600  R__ASSERT(this->fMatrix->IsValid());
601  R__ASSERT(mt->IsValid());
602  if (this->fNdiag != md.GetNdiags()) {
603  Error("operator*=(const TMatrixTDiag_const &)","matrix-diagonal's different length");
604  return;
605  }
606 
607  Element *dp1 = const_cast<Element *>(this->fPtr);
608  const Element *dp2 = md.GetPtr();
609  for (Int_t i = 0; i < this->fNdiag; i++, dp1 += this->fInc, dp2 += md.GetInc())
610  *dp1 *= *dp2;
611 }
612 
613 ////////////////////////////////////////////////////////////////////////////////
614 /// Constructor
615 
616 template<class Element>
618 {
619  R__ASSERT(matrix.IsValid());
620 
621  fMatrix = &matrix;
622  fPtr = matrix.GetMatrixArray();
623  fNelems = matrix.GetNoElements();
624 }
625 
626 ////////////////////////////////////////////////////////////////////////////////
627 /// Constructor
628 
629 template<class Element>
631 {
632  R__ASSERT(matrix.IsValid());
633 
634  fMatrix = &matrix;
635  fPtr = matrix.GetMatrixArray();
636  fNelems = matrix.GetNoElements();
637 }
638 
639 ////////////////////////////////////////////////////////////////////////////////
640 /// Constructor
641 
642 template<class Element>
644  :TMatrixTFlat_const<Element>(matrix)
645 {
646 }
647 
648 ////////////////////////////////////////////////////////////////////////////////
649 /// Constructor
650 
651 template<class Element>
653  :TMatrixTFlat_const<Element>(matrix)
654 {
655 }
656 
657 ////////////////////////////////////////////////////////////////////////////////
658 /// Copy constructor
659 
660 template<class Element>
662 {
663  *this = mf;
664 }
665 
666 ////////////////////////////////////////////////////////////////////////////////
667 /// Assign val to every element of the matrix.
668 
669 template<class Element>
671 {
672  R__ASSERT(this->fMatrix->IsValid());
673  Element *fp = const_cast<Element *>(this->fPtr);
674  while (fp < this->fPtr+this->fMatrix->GetNoElements())
675  *fp++ = val;
676 }
677 
678 ////////////////////////////////////////////////////////////////////////////////
679 /// Add val to every element of the matrix.
680 
681 template<class Element>
683 {
684  R__ASSERT(this->fMatrix->IsValid());
685  Element *fp = const_cast<Element *>(this->fPtr);
686  while (fp < this->fPtr+this->fMatrix->GetNoElements())
687  *fp++ += val;
688 }
689 
690 ////////////////////////////////////////////////////////////////////////////////
691 /// Multiply every element of the matrix with val.
692 
693 template<class Element>
695 {
696  R__ASSERT(this->fMatrix->IsValid());
697  Element *fp = const_cast<Element *>(this->fPtr);
698  while (fp < this->fPtr+this->fMatrix->GetNoElements())
699  *fp++ *= val;
700 }
701 
702 ////////////////////////////////////////////////////////////////////////////////
703 /// Assignment operator
704 
705 template<class Element>
707 {
708  const TMatrixTBase<Element> *mt = mf.GetMatrix();
709  if (this->fMatrix->GetMatrixArray() == mt->GetMatrixArray()) return;
710 
711  R__ASSERT(this->fMatrix->IsValid());
712  R__ASSERT(mt->IsValid());
713  if (this->fMatrix->GetNoElements() != mt->GetNoElements()) {
714  Error("operator=(const TMatrixTFlat_const &)","matrix lengths different");
715  return;
716  }
717 
718  Element *fp1 = const_cast<Element *>(this->fPtr);
719  const Element *fp2 = mf.GetPtr();
720  while (fp1 < this->fPtr+this->fMatrix->GetNoElements())
721  *fp1++ = *fp2++;
722 }
723 
724 ////////////////////////////////////////////////////////////////////////////////
725 /// Assign a vector to the matrix. The matrix is traversed row-wise
726 
727 template<class Element>
729 {
730  R__ASSERT(vec.IsValid());
731 
732  if (this->fMatrix->GetNoElements() != vec.GetNrows()) {
733  Error("operator=(const TVectorT &)","vector length != # matrix-elements");
734  return;
735  }
736 
737  Element *fp = const_cast<Element *>(this->fPtr);
738  const Element *vp = vec.GetMatrixArray();
739  while (fp < this->fPtr+this->fMatrix->GetNoElements())
740  *fp++ = *vp++;
741 }
742 
743 ////////////////////////////////////////////////////////////////////////////////
744 /// Add to every element of the matrix the corresponding element of matrix mf.
745 
746 template<class Element>
748 {
749  const TMatrixTBase<Element> *mt = mf.GetMatrix();
750 
751  R__ASSERT(this->fMatrix->IsValid());
752  R__ASSERT(mt->IsValid());
753  if (this->fMatrix->GetNoElements() != mt->GetNoElements()) {
754  Error("operator+=(const TMatrixTFlat_const &)","matrices lengths different");
755  return;
756  }
757 
758  Element *fp1 = const_cast<Element *>(this->fPtr);
759  const Element *fp2 = mf.GetPtr();
760  while (fp1 < this->fPtr + this->fMatrix->GetNoElements())
761  *fp1++ += *fp2++;
762 }
763 
764 ////////////////////////////////////////////////////////////////////////////////
765 /// Multiply every element of the matrix with the corresponding element of diagonal mf.
766 
767 template<class Element>
769 {
770  const TMatrixTBase<Element> *mt = mf.GetMatrix();
771 
772  R__ASSERT(this->fMatrix->IsValid());
773  R__ASSERT(mt->IsValid());
774  if (this->fMatrix->GetNoElements() != mt->GetNoElements()) {
775  Error("operator*=(const TMatrixTFlat_const &)","matrices lengths different");
776  return;
777  }
778 
779  Element *fp1 = const_cast<Element *>(this->fPtr);
780  const Element *fp2 = mf.GetPtr();
781  while (fp1 < this->fPtr + this->fMatrix->GetNoElements())
782  *fp1++ *= *fp2++;
783 }
784 
785 ////////////////////////////////////////////////////////////////////////////////
786 /// make a reference to submatrix [row_lwbs..row_upbs][col_lwbs..col_upbs];
787 /// The indexing range of the reference is
788 /// [0..row_upbs-row_lwbs+1][0..col_upb-col_lwbs+1] (default)
789 
790 template<class Element>
792  Int_t col_lwbs,Int_t col_upbs)
793 {
794  R__ASSERT(matrix.IsValid());
795 
796  fRowOff = 0;
797  fColOff = 0;
798  fNrowsSub = 0;
799  fNcolsSub = 0;
800  fMatrix = &matrix;
801 
802  if (row_upbs < row_lwbs) {
803  Error("TMatrixTSub_const","Request sub-matrix with row_upbs(%d) < row_lwbs(%d)",row_upbs,row_lwbs);
804  return;
805  }
806  if (col_upbs < col_lwbs) {
807  Error("TMatrixTSub_const","Request sub-matrix with col_upbs(%d) < col_lwbs(%d)",col_upbs,col_lwbs);
808  return;
809  }
810 
811  const Int_t rowLwb = matrix.GetRowLwb();
812  const Int_t rowUpb = matrix.GetRowUpb();
813  const Int_t colLwb = matrix.GetColLwb();
814  const Int_t colUpb = matrix.GetColUpb();
815 
816  if (row_lwbs < rowLwb || row_lwbs > rowUpb) {
817  Error("TMatrixTSub_const","Request row_lwbs sub-matrix(%d) outside matrix range of %d - %d",row_lwbs,rowLwb,rowUpb);
818  return;
819  }
820  if (col_lwbs < colLwb || col_lwbs > colUpb) {
821  Error("TMatrixTSub_const","Request col_lwbs sub-matrix(%d) outside matrix range of %d - %d",col_lwbs,colLwb,colUpb);
822  return;
823  }
824  if (row_upbs < rowLwb || row_upbs > rowUpb) {
825  Error("TMatrixTSub_const","Request row_upbs sub-matrix(%d) outside matrix range of %d - %d",row_upbs,rowLwb,rowUpb);
826  return;
827  }
828  if (col_upbs < colLwb || col_upbs > colUpb) {
829  Error("TMatrixTSub_const","Request col_upbs sub-matrix(%d) outside matrix range of %d - %d",col_upbs,colLwb,colUpb);
830  return;
831  }
832 
833  fRowOff = row_lwbs-rowLwb;
834  fColOff = col_lwbs-colLwb;
835  fNrowsSub = row_upbs-row_lwbs+1;
836  fNcolsSub = col_upbs-col_lwbs+1;
837 }
838 
839 ////////////////////////////////////////////////////////////////////////////////
840 /// make a reference to submatrix [row_lwbs..row_upbs][col_lwbs..col_upbs];
841 /// The indexing range of the reference is
842 /// [0..row_upbs-row_lwbs+1][0..col_upb-col_lwbs+1] (default)
843 
844 template<class Element>
846  Int_t col_lwbs,Int_t col_upbs)
847 {
848  R__ASSERT(matrix.IsValid());
849 
850  fRowOff = 0;
851  fColOff = 0;
852  fNrowsSub = 0;
853  fNcolsSub = 0;
854  fMatrix = &matrix;
855 
856  if (row_upbs < row_lwbs) {
857  Error("TMatrixTSub_const","Request sub-matrix with row_upbs(%d) < row_lwbs(%d)",row_upbs,row_lwbs);
858  return;
859  }
860  if (col_upbs < col_lwbs) {
861  Error("TMatrixTSub_const","Request sub-matrix with col_upbs(%d) < col_lwbs(%d)",col_upbs,col_lwbs);
862  return;
863  }
864 
865  const Int_t rowLwb = matrix.GetRowLwb();
866  const Int_t rowUpb = matrix.GetRowUpb();
867  const Int_t colLwb = matrix.GetColLwb();
868  const Int_t colUpb = matrix.GetColUpb();
869 
870  if (row_lwbs < rowLwb || row_lwbs > rowUpb) {
871  Error("TMatrixTSub_const","Request row_lwbs sub-matrix(%d) outside matrix range of %d - %d",row_lwbs,rowLwb,rowUpb);
872  return;
873  }
874  if (col_lwbs < colLwb || col_lwbs > colUpb) {
875  Error("TMatrixTSub_const","Request col_lwbs sub-matrix(%d) outside matrix range of %d - %d",col_lwbs,colLwb,colUpb);
876  return;
877  }
878  if (row_upbs < rowLwb || row_upbs > rowUpb) {
879  Error("TMatrixTSub_const","Request row_upbs sub-matrix(%d) outside matrix range of %d - %d",row_upbs,rowLwb,rowUpb);
880  return;
881  }
882  if (col_upbs < colLwb || col_upbs > colUpb) {
883  Error("TMatrixTSub_const","Request col_upbs sub-matrix(%d) outside matrix range of %d - %d",col_upbs,colLwb,colUpb);
884  return;
885  }
886 
887  fRowOff = row_lwbs-rowLwb;
888  fColOff = col_lwbs-colLwb;
889  fNrowsSub = row_upbs-row_lwbs+1;
890  fNcolsSub = col_upbs-col_lwbs+1;
891 }
892 
893 ////////////////////////////////////////////////////////////////////////////////
894 /// Constructor
895 
896 template<class Element>
898  Int_t col_lwbs,Int_t col_upbs)
899  :TMatrixTSub_const<Element>(matrix,row_lwbs,row_upbs,col_lwbs,col_upbs)
900 {
901 }
902 
903 ////////////////////////////////////////////////////////////////////////////////
904 /// Constructor
905 
906 template<class Element>
908  Int_t col_lwbs,Int_t col_upbs)
909  :TMatrixTSub_const<Element>(matrix,row_lwbs,row_upbs,col_lwbs,col_upbs)
910 {
911 }
912 
913 ////////////////////////////////////////////////////////////////////////////////
914 /// Copy constructor
915 
916 template<class Element>
918 {
919  *this = ms;
920 }
921 
922 ////////////////////////////////////////////////////////////////////////////////
923 /// Perform a rank 1 operation on the matrix:
924 /// A += alpha * v * v^T
925 
926 template<class Element>
928 {
929  R__ASSERT(this->fMatrix->IsValid());
930  R__ASSERT(v.IsValid());
931 
932  if (v.GetNoElements() < TMath::Max(this->fNrowsSub,this->fNcolsSub)) {
933  Error("Rank1Update","vector too short");
934  return;
935  }
936 
937  const Element * const pv = v.GetMatrixArray();
938  Element *mp = (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->GetMatrixArray();
939 
940  const Int_t ncols = this->fMatrix->GetNcols();
941  for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
942  const Int_t off = (irow+this->fRowOff)*ncols+this->fColOff;
943  const Element tmp = alpha*pv[irow];
944  for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
945  mp[off+icol] += tmp*pv[icol];
946  }
947 }
948 
949 ////////////////////////////////////////////////////////////////////////////////
950 /// Assign val to every element of the sub matrix.
951 
952 template<class Element>
954 {
955  R__ASSERT(this->fMatrix->IsValid());
956 
957  Element *p = (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->GetMatrixArray();
958  const Int_t ncols = this->fMatrix->GetNcols();
959  for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
960  const Int_t off = (irow+this->fRowOff)*ncols+this->fColOff;
961  for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
962  p[off+icol] = val;
963  }
964 }
965 
966 ////////////////////////////////////////////////////////////////////////////////
967 /// Add val to every element of the sub matrix.
968 
969 template<class Element>
971 {
972  R__ASSERT(this->fMatrix->IsValid());
973 
974  Element *p = (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->GetMatrixArray();
975  const Int_t ncols = this->fMatrix->GetNcols();
976  for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
977  const Int_t off = (irow+this->fRowOff)*ncols+this->fColOff;
978  for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
979  p[off+icol] += val;
980  }
981 }
982 
983 ////////////////////////////////////////////////////////////////////////////////
984 /// Multiply every element of the sub matrix by val .
985 
986 template<class Element>
988 {
989  R__ASSERT(this->fMatrix->IsValid());
990 
991  Element *p = (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->GetMatrixArray();
992  const Int_t ncols = this->fMatrix->GetNcols();
993  for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
994  const Int_t off = (irow+this->fRowOff)*ncols+this->fColOff;
995  for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
996  p[off+icol] *= val;
997  }
998 }
999 
1000 ////////////////////////////////////////////////////////////////////////////////
1001 /// Assignment operator
1002 
1003 template<class Element>
1005 {
1006  const TMatrixTBase<Element> *mt = ms.GetMatrix();
1007 
1008  R__ASSERT(this->fMatrix->IsValid());
1009  R__ASSERT(mt->IsValid());
1010 
1011  if (this->fMatrix == mt &&
1012  (this->GetNrows() == ms.GetNrows () && this->GetNcols() == ms.GetNcols () &&
1013  this->GetRowOff() == ms.GetRowOff() && this->GetColOff() == ms.GetColOff()) )
1014  return;
1015 
1016  if (this->GetNrows() != ms.GetNrows() || this->GetNcols() != ms.GetNcols()) {
1017  Error("operator=(const TMatrixTSub_const &)","sub matrices have different size");
1018  return;
1019  }
1020 
1021  const Int_t rowOff2 = ms.GetRowOff();
1022  const Int_t colOff2 = ms.GetColOff();
1023 
1024  Bool_t overlap = (this->fMatrix == mt) &&
1025  ( (rowOff2 >= this->fRowOff && rowOff2 < this->fRowOff+this->fNrowsSub) ||
1026  (colOff2 >= this->fColOff && colOff2 < this->fColOff+this->fNcolsSub) );
1027 
1028  Element *p1 = const_cast<Element *>(this->fMatrix->GetMatrixArray());
1029  if (!overlap) {
1030  const Element *p2 = mt->GetMatrixArray();
1031 
1032  const Int_t ncols1 = this->fMatrix->GetNcols();
1033  const Int_t ncols2 = mt->GetNcols();
1034  for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1035  const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1036  const Int_t off2 = (irow+rowOff2)*ncols2+colOff2;
1037  for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1038  p1[off1+icol] = p2[off2+icol];
1039  }
1040  } else {
1041  const Int_t row_lwbs = rowOff2+mt->GetRowLwb();
1042  const Int_t row_upbs = row_lwbs+this->fNrowsSub-1;
1043  const Int_t col_lwbs = colOff2+mt->GetColLwb();
1044  const Int_t col_upbs = col_lwbs+this->fNcolsSub-1;
1045  TMatrixT<Element> tmp; mt->GetSub(row_lwbs,row_upbs,col_lwbs,col_upbs,tmp);
1046  const Element *p2 = tmp.GetMatrixArray();
1047 
1048  const Int_t ncols1 = this->fMatrix->GetNcols();
1049  const Int_t ncols2 = tmp.GetNcols();
1050  for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1051  const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1052  const Int_t off2 = irow*ncols2;
1053  for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1054  p1[off1+icol] = p2[off2+icol];
1055  }
1056  }
1057 }
1058 
1059 ////////////////////////////////////////////////////////////////////////////////
1060 /// Assignment operator
1061 
1062 template<class Element>
1064 {
1065  R__ASSERT(this->fMatrix->IsValid());
1066  R__ASSERT(m.IsValid());
1067 
1068  if (this->fMatrix->GetMatrixArray() == m.GetMatrixArray()) return;
1069 
1070  if (this->fNrowsSub != m.GetNrows() || this->fNcolsSub != m.GetNcols()) {
1071  Error("operator=(const TMatrixTBase<Element> &)","sub matrices and matrix have different size");
1072  return;
1073  }
1074  const Int_t row_lwbs = this->fRowOff+this->fMatrix->GetRowLwb();
1075  const Int_t col_lwbs = this->fColOff+this->fMatrix->GetColLwb();
1076  (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->SetSub(row_lwbs,col_lwbs,m);
1077 }
1078 
1079 ////////////////////////////////////////////////////////////////////////////////
1080 /// Add to every element of the submatrix the corresponding element of submatrix ms.
1081 
1082 template<class Element>
1084 {
1085  const TMatrixTBase<Element> *mt = ms.GetMatrix();
1086 
1087  R__ASSERT(this->fMatrix->IsValid());
1088  R__ASSERT(mt->IsValid());
1089 
1090  if (this->GetNrows() != ms.GetNrows() || this->GetNcols() != ms.GetNcols()) {
1091  Error("operator+=(const TMatrixTSub_const &)","sub matrices have different size");
1092  return;
1093  }
1094 
1095  const Int_t rowOff2 = ms.GetRowOff();
1096  const Int_t colOff2 = ms.GetColOff();
1097 
1098  Bool_t overlap = (this->fMatrix == mt) &&
1099  ( (rowOff2 >= this->fRowOff && rowOff2 < this->fRowOff+this->fNrowsSub) ||
1100  (colOff2 >= this->fColOff && colOff2 < this->fColOff+this->fNcolsSub) );
1101 
1102  Element *p1 = const_cast<Element *>(this->fMatrix->GetMatrixArray());
1103  if (!overlap) {
1104  const Element *p2 = mt->GetMatrixArray();
1105 
1106  const Int_t ncols1 = this->fMatrix->GetNcols();
1107  const Int_t ncols2 = mt->GetNcols();
1108  for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1109  const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1110  const Int_t off2 = (irow+rowOff2)*ncols2+colOff2;
1111  for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1112  p1[off1+icol] += p2[off2+icol];
1113  }
1114  } else {
1115  const Int_t row_lwbs = rowOff2+mt->GetRowLwb();
1116  const Int_t row_upbs = row_lwbs+this->fNrowsSub-1;
1117  const Int_t col_lwbs = colOff2+mt->GetColLwb();
1118  const Int_t col_upbs = col_lwbs+this->fNcolsSub-1;
1119  TMatrixT<Element> tmp; mt->GetSub(row_lwbs,row_upbs,col_lwbs,col_upbs,tmp);
1120  const Element *p2 = tmp.GetMatrixArray();
1121 
1122  const Int_t ncols1 = this->fMatrix->GetNcols();
1123  const Int_t ncols2 = tmp.GetNcols();
1124  for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1125  const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1126  const Int_t off2 = irow*ncols2;
1127  for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1128  p1[off1+icol] += p2[off2+icol];
1129  }
1130  }
1131 }
1132 
1133 ////////////////////////////////////////////////////////////////////////////////
1134 /// Multiply submatrix with submatrix ms.
1135 
1136 template<class Element>
1138 {
1139  if (this->fNcolsSub != ms.GetNrows() || this->fNcolsSub != ms.GetNcols()) {
1140  Error("operator*=(const TMatrixTSub_const &)","source sub matrix has wrong shape");
1141  return;
1142  }
1143 
1144  const TMatrixTBase<Element> *source = ms.GetMatrix();
1145 
1146  TMatrixT<Element> source_sub;
1147  {
1148  const Int_t row_lwbs = ms.GetRowOff()+source->GetRowLwb();
1149  const Int_t row_upbs = row_lwbs+this->fNrowsSub-1;
1150  const Int_t col_lwbs = ms.GetColOff()+source->GetColLwb();
1151  const Int_t col_upbs = col_lwbs+this->fNcolsSub-1;
1152  source->GetSub(row_lwbs,row_upbs,col_lwbs,col_upbs,source_sub);
1153  }
1154 
1155  const Element *sp = source_sub.GetMatrixArray();
1156  const Int_t ncols = this->fMatrix->GetNcols();
1157 
1158  // One row of the old_target matrix
1159  Element work[kWorkMax];
1160  Bool_t isAllocated = kFALSE;
1161  Element *trp = work;
1162  if (this->fNcolsSub > kWorkMax) {
1163  isAllocated = kTRUE;
1164  trp = new Element[this->fNcolsSub];
1165  }
1166 
1167  Element *cp = const_cast<Element *>(this->fMatrix->GetMatrixArray())+this->fRowOff*ncols+this->fColOff;
1168  const Element *trp0 = cp; // Pointer to target[i,0];
1169  const Element * const trp0_last = trp0+this->fNrowsSub*ncols;
1170  while (trp0 < trp0_last) {
1171  memcpy(trp,trp0,this->fNcolsSub*sizeof(Element)); // copy the i-th row of target, Start at target[i,0]
1172  for (const Element *scp = sp; scp < sp+this->fNcolsSub; ) { // Pointer to the j-th column of source,
1173  // Start scp = source[0,0]
1174  Element cij = 0;
1175  for (Int_t j = 0; j < this->fNcolsSub; j++) {
1176  cij += trp[j] * *scp; // the j-th col of source
1177  scp += this->fNcolsSub;
1178  }
1179  *cp++ = cij;
1180  scp -= source_sub.GetNoElements()-1; // Set bcp to the (j+1)-th col
1181  }
1182  cp += ncols-this->fNcolsSub;
1183  trp0 += ncols; // Set trp0 to the (i+1)-th row
1184  R__ASSERT(trp0 == cp);
1185  }
1186 
1187  R__ASSERT(cp == trp0_last && trp0 == trp0_last);
1188  if (isAllocated)
1189  delete [] trp;
1190 }
1191 
1192 ////////////////////////////////////////////////////////////////////////////////
1193 /// Add to every element of the submatrix the corresponding element of matrix mt.
1194 
1195 template<class Element>
1197 {
1198  R__ASSERT(this->fMatrix->IsValid());
1199  R__ASSERT(mt.IsValid());
1200 
1201  if (this->GetNrows() != mt.GetNrows() || this->GetNcols() != mt.GetNcols()) {
1202  Error("operator+=(const TMatrixTBase<Element> &)","sub matrix and matrix have different size");
1203  return;
1204  }
1205 
1206  Element *p1 = const_cast<Element *>(this->fMatrix->GetMatrixArray());
1207  const Element *p2 = mt.GetMatrixArray();
1208 
1209  const Int_t ncols1 = this->fMatrix->GetNcols();
1210  const Int_t ncols2 = mt.GetNcols();
1211  for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1212  const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1213  const Int_t off2 = irow*ncols2;
1214  for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1215  p1[off1+icol] += p2[off2+icol];
1216  }
1217 }
1218 
1219 ////////////////////////////////////////////////////////////////////////////////
1220 /// Multiply submatrix with matrix source.
1221 
1222 template<class Element>
1224 {
1225  if (this->fNcolsSub != source.GetNrows() || this->fNcolsSub != source.GetNcols()) {
1226  Error("operator*=(const TMatrixT<Element> &)","source matrix has wrong shape");
1227  return;
1228  }
1229 
1230  // Check for A *= A;
1231  const Element *sp;
1232  TMatrixT<Element> tmp;
1233  if (this->fMatrix->GetMatrixArray() == source.GetMatrixArray()) {
1234  tmp.ResizeTo(source);
1235  tmp = source;
1236  sp = tmp.GetMatrixArray();
1237  }
1238  else
1239  sp = source.GetMatrixArray();
1240 
1241  const Int_t ncols = this->fMatrix->GetNcols();
1242 
1243  // One row of the old_target matrix
1244  Element work[kWorkMax];
1245  Bool_t isAllocated = kFALSE;
1246  Element *trp = work;
1247  if (this->fNcolsSub > kWorkMax) {
1248  isAllocated = kTRUE;
1249  trp = new Element[this->fNcolsSub];
1250  }
1251 
1252  Element *cp = const_cast<Element *>(this->fMatrix->GetMatrixArray())+this->fRowOff*ncols+this->fColOff;
1253  const Element *trp0 = cp; // Pointer to target[i,0];
1254  const Element * const trp0_last = trp0+this->fNrowsSub*ncols;
1255  while (trp0 < trp0_last) {
1256  memcpy(trp,trp0,this->fNcolsSub*sizeof(Element)); // copy the i-th row of target, Start at target[i,0]
1257  for (const Element *scp = sp; scp < sp+this->fNcolsSub; ) { // Pointer to the j-th column of source,
1258  // Start scp = source[0,0]
1259  Element cij = 0;
1260  for (Int_t j = 0; j < this->fNcolsSub; j++) {
1261  cij += trp[j] * *scp; // the j-th col of source
1262  scp += this->fNcolsSub;
1263  }
1264  *cp++ = cij;
1265  scp -= source.GetNoElements()-1; // Set bcp to the (j+1)-th col
1266  }
1267  cp += ncols-this->fNcolsSub;
1268  trp0 += ncols; // Set trp0 to the (i+1)-th row
1269  R__ASSERT(trp0 == cp);
1270  }
1271 
1272  R__ASSERT(cp == trp0_last && trp0 == trp0_last);
1273  if (isAllocated)
1274  delete [] trp;
1275 }
1276 
1277 ////////////////////////////////////////////////////////////////////////////////
1278 /// Multiply submatrix with matrix source.
1279 
1280 template<class Element>
1282 {
1283  if (this->fNcolsSub != source.GetNrows() || this->fNcolsSub != source.GetNcols()) {
1284  Error("operator*=(const TMatrixTSym<Element> &)","source matrix has wrong shape");
1285  return;
1286  }
1287 
1288  // Check for A *= A;
1289  const Element *sp;
1291  if (this->fMatrix->GetMatrixArray() == source.GetMatrixArray()) {
1292  tmp.ResizeTo(source);
1293  tmp = source;
1294  sp = tmp.GetMatrixArray();
1295  }
1296  else
1297  sp = source.GetMatrixArray();
1298 
1299  const Int_t ncols = this->fMatrix->GetNcols();
1300 
1301  // One row of the old_target matrix
1302  Element work[kWorkMax];
1303  Bool_t isAllocated = kFALSE;
1304  Element *trp = work;
1305  if (this->fNcolsSub > kWorkMax) {
1306  isAllocated = kTRUE;
1307  trp = new Element[this->fNcolsSub];
1308  }
1309 
1310  Element *cp = const_cast<Element *>(this->fMatrix->GetMatrixArray())+this->fRowOff*ncols+this->fColOff;
1311  const Element *trp0 = cp; // Pointer to target[i,0];
1312  const Element * const trp0_last = trp0+this->fNrowsSub*ncols;
1313  while (trp0 < trp0_last) {
1314  memcpy(trp,trp0,this->fNcolsSub*sizeof(Element)); // copy the i-th row of target, Start at target[i,0]
1315  for (const Element *scp = sp; scp < sp+this->fNcolsSub; ) { // Pointer to the j-th column of source,
1316  // Start scp = source[0,0]
1317  Element cij = 0;
1318  for (Int_t j = 0; j < this->fNcolsSub; j++) {
1319  cij += trp[j] * *scp; // the j-th col of source
1320  scp += this->fNcolsSub;
1321  }
1322  *cp++ = cij;
1323  scp -= source.GetNoElements()-1; // Set bcp to the (j+1)-th col
1324  }
1325  cp += ncols-this->fNcolsSub;
1326  trp0 += ncols; // Set trp0 to the (i+1)-th row
1327  R__ASSERT(trp0 == cp);
1328  }
1329 
1330  R__ASSERT(cp == trp0_last && trp0 == trp0_last);
1331  if (isAllocated)
1332  delete [] trp;
1333 }
1334 
1335 ////////////////////////////////////////////////////////////////////////////////
1336 /// Constructor with row "row" of matrix
1337 
1338 template<class Element>
1340 {
1341  R__ASSERT(matrix.IsValid());
1342 
1343  fRowInd = row-matrix.GetRowLwb();
1344  if (fRowInd >= matrix.GetNrows() || fRowInd < 0) {
1345  Error("TMatrixTSparseRow_const(const TMatrixTSparse &,Int_t)","row index out of bounds");
1346  fMatrix = 0;
1347  fNindex = 0;
1348  fColPtr = 0;
1349  fDataPtr = 0;
1350  return;
1351  }
1352 
1353  const Int_t sIndex = matrix.GetRowIndexArray()[fRowInd];
1354  const Int_t eIndex = matrix.GetRowIndexArray()[fRowInd+1];
1355  fMatrix = &matrix;
1356  fNindex = eIndex-sIndex;
1357  fColPtr = matrix.GetColIndexArray()+sIndex;
1358  fDataPtr = matrix.GetMatrixArray()+sIndex;
1359 }
1360 
1361 ////////////////////////////////////////////////////////////////////////////////
1362 
1363 template<class Element>
1365 {
1367  R__ASSERT(fMatrix->IsValid());
1368  const Int_t acoln = i-fMatrix->GetColLwb();
1369  if (acoln < fMatrix->GetNcols() && acoln >= 0) {
1370  const Int_t index = TMath::BinarySearch(fNindex,fColPtr,acoln);
1371  if (index >= 0 && fColPtr[index] == acoln) return fDataPtr[index];
1372  else return 0.0;
1373  } else {
1374  Error("operator()","Request col(%d) outside matrix range of %d - %d",
1375  i,fMatrix->GetColLwb(),fMatrix->GetColLwb()+fMatrix->GetNcols());
1377  }
1378  }
1379 
1380 ////////////////////////////////////////////////////////////////////////////////
1381 /// Constructor with row "row" of matrix
1382 
1383 template<class Element>
1385  : TMatrixTSparseRow_const<Element>(matrix,row)
1386 {
1387 }
1388 
1389 ////////////////////////////////////////////////////////////////////////////////
1390 /// Copy constructor
1391 
1392 template<class Element>
1394  : TMatrixTSparseRow_const<Element>(mr)
1395 {
1396  *this = mr;
1397 }
1398 
1399 ////////////////////////////////////////////////////////////////////////////////
1400 
1401 template<class Element>
1403 {
1404  if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
1405  R__ASSERT(this->fMatrix->IsValid());
1406  const Int_t acoln = i-this->fMatrix->GetColLwb();
1407  if (acoln < this->fMatrix->GetNcols() && acoln >= 0) {
1408  const Int_t index = TMath::BinarySearch(this->fNindex,this->fColPtr,acoln);
1409  if (index >= 0 && this->fColPtr[index] == acoln) return this->fDataPtr[index];
1410  else return 0.0;
1411  } else {
1412  Error("operator()","Request col(%d) outside matrix range of %d - %d",
1413  i,this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
1415  }
1416 }
1417 
1418 ////////////////////////////////////////////////////////////////////////////////
1419 /// operator() : pick element row(i)
1420 
1421 template<class Element>
1423 {
1424  if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
1425  R__ASSERT(this->fMatrix->IsValid());
1426 
1427  const Int_t acoln = i-this->fMatrix->GetColLwb();
1428  if (acoln >= this->fMatrix->GetNcols() || acoln < 0) {
1429  Error("operator()(Int_t","Requested element %d outside range : %d - %d",i,
1430  this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
1432  }
1433 
1434  Int_t index = TMath::BinarySearch(this->fNindex,this->fColPtr,acoln);
1435  if (index >= 0 && this->fColPtr[index] == acoln)
1436  return (const_cast<Element*>(this->fDataPtr))[index];
1437  else {
1438  TMatrixTSparse<Element> *mt = const_cast<TMatrixTSparse<Element> *>(this->fMatrix);
1439  const Int_t row = this->fRowInd+mt->GetRowLwb();
1440  Element val = 0.;
1441  mt->InsertRow(row,i,&val,1);
1442  const Int_t sIndex = mt->GetRowIndexArray()[this->fRowInd];
1443  const Int_t eIndex = mt->GetRowIndexArray()[this->fRowInd+1];
1444  this->fNindex = eIndex-sIndex;
1445  this->fColPtr = mt->GetColIndexArray()+sIndex;
1446  this->fDataPtr = mt->GetMatrixArray()+sIndex;
1447  index = TMath::BinarySearch(this->fNindex,this->fColPtr,acoln);
1448  if (index >= 0 && this->fColPtr[index] == acoln)
1449  return (const_cast<Element*>(this->fDataPtr))[index];
1450  else {
1451  Error("operator()(Int_t","Insert row failed");
1453  }
1454  }
1455 }
1456 
1457 ////////////////////////////////////////////////////////////////////////////////
1458 /// Assign val to every non-zero (!) element of the matrix row.
1459 
1460 template<class Element>
1462 {
1463  R__ASSERT(this->fMatrix->IsValid());
1464  Element *rp = const_cast<Element *>(this->fDataPtr);
1465  for ( ; rp < this->fDataPtr+this->fNindex; rp++)
1466  *rp = val;
1467 }
1468 
1469 ////////////////////////////////////////////////////////////////////////////////
1470 /// Add val to every non-zero (!) element of the matrix row.
1471 
1472 template<class Element>
1474 {
1475  R__ASSERT(this->fMatrix->IsValid());
1476  Element *rp = const_cast<Element *>(this->fDataPtr);
1477  for ( ; rp < this->fDataPtr+this->fNindex; rp++)
1478  *rp += val;
1479 }
1480 
1481 ////////////////////////////////////////////////////////////////////////////////
1482 /// Multiply every element of the matrix row by val.
1483 
1484 template<class Element>
1486 {
1487  R__ASSERT(this->fMatrix->IsValid());
1488  Element *rp = const_cast<Element *>(this->fDataPtr);
1489  for ( ; rp < this->fDataPtr+this->fNindex; rp++)
1490  *rp *= val;
1491 }
1492 
1493 ////////////////////////////////////////////////////////////////////////////////
1494 /// Assignment operator
1495 
1496 template<class Element>
1498 {
1499  const TMatrixTBase<Element> *mt = mr.GetMatrix();
1500  if (this->fMatrix == mt) return;
1501 
1502  R__ASSERT(this->fMatrix->IsValid());
1503  R__ASSERT(mt->IsValid());
1504  if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
1505  Error("operator=(const TMatrixTSparseRow_const &)","matrix rows not compatible");
1506  return;
1507  }
1508 
1509  const Int_t ncols = this->fMatrix->GetNcols();
1510  const Int_t row1 = this->fRowInd+this->fMatrix->GetRowLwb();
1511  const Int_t row2 = mr.GetRowIndex()+mt->GetRowLwb();
1512  const Int_t col = this->fMatrix->GetColLwb();
1513 
1514  TVectorT<Element> v(ncols);
1515  mt->ExtractRow(row2,col,v.GetMatrixArray());
1516  const_cast<TMatrixTSparse<Element> *>(this->fMatrix)->InsertRow(row1,col,v.GetMatrixArray());
1517 
1518  const Int_t sIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd];
1519  const Int_t eIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd+1];
1520  this->fNindex = eIndex-sIndex;
1521  this->fColPtr = this->fMatrix->GetColIndexArray()+sIndex;
1522  this->fDataPtr = this->fMatrix->GetMatrixArray()+sIndex;
1523 }
1524 
1525 ////////////////////////////////////////////////////////////////////////////////
1526 /// Assign a vector to a matrix row. The vector is considered row-vector
1527 /// to allow the assignment in the strict sense.
1528 
1529 template<class Element>
1531 {
1532  R__ASSERT(this->fMatrix->IsValid());
1533  R__ASSERT(vec.IsValid());
1534 
1535  if (this->fMatrix->GetColLwb() != vec.GetLwb() || this->fMatrix->GetNcols() != vec.GetNrows()) {
1536  Error("operator=(const TVectorT &)","vector length != matrix-row length");
1537  return;
1538  }
1539 
1540  const Element *vp = vec.GetMatrixArray();
1541  const Int_t row = this->fRowInd+this->fMatrix->GetRowLwb();
1542  const Int_t col = this->fMatrix->GetColLwb();
1543  const_cast<TMatrixTSparse<Element> *>(this->fMatrix)->InsertRow(row,col,vp,vec.GetNrows());
1544 
1545  const Int_t sIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd];
1546  const Int_t eIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd+1];
1547  this->fNindex = eIndex-sIndex;
1548  this->fColPtr = this->fMatrix->GetColIndexArray()+sIndex;
1549  this->fDataPtr = this->fMatrix->GetMatrixArray()+sIndex;
1550 }
1551 
1552 ////////////////////////////////////////////////////////////////////////////////
1553 /// Add to every element of the matrix row the corresponding element of row r.
1554 
1555 template<class Element>
1557 {
1558  const TMatrixTBase<Element> *mt = r.GetMatrix();
1559 
1560  R__ASSERT(this->fMatrix->IsValid());
1561  R__ASSERT(mt->IsValid());
1562  if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
1563  Error("operator+=(const TMatrixTRow_const &)","different row lengths");
1564  return;
1565  }
1566 
1567  const Int_t ncols = this->fMatrix->GetNcols();
1568  const Int_t row1 = this->fRowInd+this->fMatrix->GetRowLwb();
1569  const Int_t row2 = r.GetRowIndex()+mt->GetRowLwb();
1570  const Int_t col = this->fMatrix->GetColLwb();
1571 
1572  TVectorT<Element> v1(ncols);
1573  TVectorT<Element> v2(ncols);
1574  this->fMatrix->ExtractRow(row1,col,v1.GetMatrixArray());
1575  mt ->ExtractRow(row2,col,v2.GetMatrixArray());
1576  v1 += v2;
1577  const_cast<TMatrixTSparse<Element> *>(this->fMatrix)->InsertRow(row1,col,v1.GetMatrixArray());
1578 
1579  const Int_t sIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd];
1580  const Int_t eIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd+1];
1581  this->fNindex = eIndex-sIndex;
1582  this->fColPtr = this->fMatrix->GetColIndexArray()+sIndex;
1583  this->fDataPtr = this->fMatrix->GetMatrixArray()+sIndex;
1584 }
1585 
1586 ////////////////////////////////////////////////////////////////////////////////
1587 /// Multiply every element of the matrix row with the
1588 /// corresponding element of row r.
1589 
1590 template<class Element>
1592 {
1593  const TMatrixTBase<Element> *mt = r.GetMatrix();
1594 
1595  R__ASSERT(this->fMatrix->IsValid());
1596  R__ASSERT(mt->IsValid());
1597  if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
1598  Error("operator+=(const TMatrixTRow_const &)","different row lengths");
1599  return;
1600  }
1601 
1602  const Int_t ncols = this->fMatrix->GetNcols();
1603  const Int_t row1 = r.GetRowIndex()+mt->GetRowLwb();
1604  const Int_t row2 = r.GetRowIndex()+mt->GetRowLwb();
1605  const Int_t col = this->fMatrix->GetColLwb();
1606 
1607  TVectorT<Element> v1(ncols);
1608  TVectorT<Element> v2(ncols);
1609  this->fMatrix->ExtractRow(row1,col,v1.GetMatrixArray());
1610  mt ->ExtractRow(row2,col,v2.GetMatrixArray());
1611 
1612  ElementMult(v1,v2);
1613  const_cast<TMatrixTSparse<Element> *>(this->fMatrix)->InsertRow(row1,col,v1.GetMatrixArray());
1614 
1615  const Int_t sIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd];
1616  const Int_t eIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd+1];
1617  this->fNindex = eIndex-sIndex;
1618  this->fColPtr = this->fMatrix->GetColIndexArray()+sIndex;
1619  this->fDataPtr = this->fMatrix->GetMatrixArray()+sIndex;
1620 }
1621 
1622 ////////////////////////////////////////////////////////////////////////////////
1623 /// Constructor
1624 
1625 template<class Element>
1627 {
1628  R__ASSERT(matrix.IsValid());
1629 
1630  fMatrix = &matrix;
1631  fNdiag = TMath::Min(matrix.GetNrows(),matrix.GetNcols());
1632  fDataPtr = matrix.GetMatrixArray();
1633 }
1634 
1635 ////////////////////////////////////////////////////////////////////////////////
1636 
1637 template<class Element>
1639 {
1640  R__ASSERT(fMatrix->IsValid());
1641  if (i < fNdiag && i >= 0) {
1642  const Int_t * const pR = fMatrix->GetRowIndexArray();
1643  const Int_t * const pC = fMatrix->GetColIndexArray();
1644  const Element * const pD = fMatrix->GetMatrixArray();
1645  const Int_t sIndex = pR[i];
1646  const Int_t eIndex = pR[i+1];
1647  const Int_t index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex;
1648  if (index >= sIndex && pC[index] == i) return pD[index];
1649  else return 0.0;
1650  } else {
1651  Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,fNdiag);
1652  return 0.0;
1653  }
1654  return 0.0;
1655 }
1656 
1657 ////////////////////////////////////////////////////////////////////////////////
1658 /// Constructor
1659 
1660 template<class Element>
1662  :TMatrixTSparseDiag_const<Element>(matrix)
1663 {
1664 }
1665 
1666 ////////////////////////////////////////////////////////////////////////////////
1667 /// Constructor
1668 
1669 template<class Element>
1671  : TMatrixTSparseDiag_const<Element>(md)
1672 {
1673  *this = md;
1674 }
1675 
1676 ////////////////////////////////////////////////////////////////////////////////
1677 
1678 template<class Element>
1680 {
1681  R__ASSERT(this->fMatrix->IsValid());
1682  if (i < this->fNdiag && i >= 0) {
1683  const Int_t * const pR = this->fMatrix->GetRowIndexArray();
1684  const Int_t * const pC = this->fMatrix->GetColIndexArray();
1685  const Element * const pD = this->fMatrix->GetMatrixArray();
1686  const Int_t sIndex = pR[i];
1687  const Int_t eIndex = pR[i+1];
1688  const Int_t index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex;
1689  if (index >= sIndex && pC[index] == i) return pD[index];
1690  else return 0.0;
1691  } else {
1692  Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,this->fNdiag);
1693  return 0.0;
1694  }
1695  return 0.0;
1696  }
1697 
1698 ////////////////////////////////////////////////////////////////////////////////
1699 /// operator() : pick element diag(i)
1700 
1701 template<class Element>
1703 {
1704  R__ASSERT(this->fMatrix->IsValid());
1705 
1706  if (i < 0 || i >= this->fNdiag) {
1707  Error("operator()(Int_t","Requested element %d outside range : 0 - %d",i,this->fNdiag);
1708  return (const_cast<Element*>(this->fDataPtr))[0];
1709  }
1710 
1711  TMatrixTSparse<Element> *mt = const_cast<TMatrixTSparse<Element> *>(this->fMatrix);
1712  const Int_t *pR = mt->GetRowIndexArray();
1713  const Int_t *pC = mt->GetColIndexArray();
1714  Int_t sIndex = pR[i];
1715  Int_t eIndex = pR[i+1];
1716  Int_t index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex;
1717  if (index >= sIndex && pC[index] == i)
1718  return (const_cast<Element*>(this->fDataPtr))[index];
1719  else {
1720  const Int_t row = i+mt->GetRowLwb();
1721  const Int_t col = i+mt->GetColLwb();
1722  Element val = 0.;
1723  mt->InsertRow(row,col,&val,1);
1724  this->fDataPtr = mt->GetMatrixArray();
1725  pR = mt->GetRowIndexArray();
1726  pC = mt->GetColIndexArray();
1727  sIndex = pR[i];
1728  eIndex = pR[i+1];
1729  index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex;
1730  if (index >= sIndex && pC[index] == i)
1731  return (const_cast<Element*>(this->fDataPtr))[index];
1732  else {
1733  Error("operator()(Int_t","Insert row failed");
1734  return (const_cast<Element*>(this->fDataPtr))[0];
1735  }
1736  }
1737 }
1738 
1739 ////////////////////////////////////////////////////////////////////////////////
1740 /// Assign val to every element of the matrix diagonal.
1741 
1742 template<class Element>
1744 {
1745  R__ASSERT(this->fMatrix->IsValid());
1746  for (Int_t i = 0; i < this->fNdiag; i++)
1747  (*this)(i) = val;
1748 }
1749 
1750 ////////////////////////////////////////////////////////////////////////////////
1751 /// Add val to every element of the matrix diagonal.
1752 
1753 template<class Element>
1755 {
1756  R__ASSERT(this->fMatrix->IsValid());
1757  for (Int_t i = 0; i < this->fNdiag; i++)
1758  (*this)(i) += val;
1759 }
1760 
1761 ////////////////////////////////////////////////////////////////////////////////
1762 /// Multiply every element of the matrix diagonal by val.
1763 
1764 template<class Element>
1766 {
1767  R__ASSERT(this->fMatrix->IsValid());
1768  for (Int_t i = 0; i < this->fNdiag; i++)
1769  (*this)(i) *= val;
1770 }
1771 
1772 ////////////////////////////////////////////////////////////////////////////////
1773 /// Assignment operator
1774 
1775 template<class Element>
1777 {
1778  const TMatrixTBase<Element> *mt = md.GetMatrix();
1779  if (this->fMatrix == mt) return;
1780 
1781  R__ASSERT(this->fMatrix->IsValid());
1782  R__ASSERT(mt->IsValid());
1783  if (this->fNdiag != md.GetNdiags()) {
1784  Error("operator=(const TMatrixTSparseDiag_const &)","matrix-diagonal's different length");
1785  return;
1786  }
1787 
1788  for (Int_t i = 0; i < this->fNdiag; i++)
1789  (*this)(i) = md(i);
1790 }
1791 
1792 ////////////////////////////////////////////////////////////////////////////////
1793 /// Assign a vector to the matrix diagonal.
1794 
1795 template<class Element>
1797 {
1798  R__ASSERT(this->fMatrix->IsValid());
1799  R__ASSERT(vec.IsValid());
1800 
1801  if (this->fNdiag != vec.GetNrows()) {
1802  Error("operator=(const TVectorT &)","vector length != matrix-diagonal length");
1803  return;
1804  }
1805 
1806  const Element *vp = vec.GetMatrixArray();
1807  for (Int_t i = 0; i < this->fNdiag; i++)
1808  (*this)(i) = vp[i];
1809 }
1810 
1811 ////////////////////////////////////////////////////////////////////////////////
1812 /// Add to every element of the matrix diagonal the
1813 /// corresponding element of diagonal md.
1814 
1815 template<class Element>
1817 {
1818  const TMatrixTBase<Element> *mt = md.GetMatrix();
1819 
1820  R__ASSERT(this->fMatrix->IsValid());
1821  R__ASSERT(mt->IsValid());
1822  if (this->fNdiag != md.GetNdiags()) {
1823  Error("operator+=(const TMatrixTSparseDiag_const &)","matrix-diagonal's different length");
1824  return;
1825  }
1826 
1827  for (Int_t i = 0; i < this->fNdiag; i++)
1828  (*this)(i) += md(i);
1829 }
1830 
1831 ////////////////////////////////////////////////////////////////////////////////
1832 /// Multiply every element of the matrix diagonal with the
1833 /// corresponding element of diagonal md.
1834 
1835 template<class Element>
1837 {
1838  const TMatrixTBase<Element> *mt = md.GetMatrix();
1839 
1840  R__ASSERT(this->fMatrix->IsValid());
1841  R__ASSERT(mt->IsValid());
1842  if (this->fNdiag != md.GetNdiags()) {
1843  Error("operator*=(const TMatrixTSparseDiag_const &)","matrix-diagonal's different length");
1844  return;
1845  }
1846 
1847  for (Int_t i = 0; i < this->fNdiag; i++)
1848  (*this)(i) *= md(i);
1849 }
1850 
1851 ////////////////////////////////////////////////////////////////////////////////
1852 /// Random number generator [0....1] with seed ix
1853 
1855 {
1856  const Double_t a = 16807.0;
1857  const Double_t b15 = 32768.0;
1858  const Double_t b16 = 65536.0;
1859  const Double_t p = 2147483647.0;
1860  Double_t xhi = ix/b16;
1861  Int_t xhiint = (Int_t) xhi;
1862  xhi = xhiint;
1863  Double_t xalo = (ix-xhi*b16)*a;
1864 
1865  Double_t leftlo = xalo/b16;
1866  Int_t leftloint = (int) leftlo;
1867  leftlo = leftloint;
1868  Double_t fhi = xhi*a+leftlo;
1869  Double_t k = fhi/b15;
1870  Int_t kint = (Int_t) k;
1871  k = kint;
1872  ix = (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k;
1873  if (ix < 0.0) ix = ix+p;
1874 
1875  return (ix*4.656612875e-10);
1876 }
1877 
1878 template class TMatrixTRow_const <Float_t>;
1879 template class TMatrixTColumn_const <Float_t>;
1880 template class TMatrixTDiag_const <Float_t>;
1881 template class TMatrixTFlat_const <Float_t>;
1882 template class TMatrixTSub_const <Float_t>;
1883 template class TMatrixTSparseRow_const <Float_t>;
1884 template class TMatrixTSparseDiag_const<Float_t>;
1885 template class TMatrixTRow <Float_t>;
1886 template class TMatrixTColumn <Float_t>;
1887 template class TMatrixTDiag <Float_t>;
1888 template class TMatrixTFlat <Float_t>;
1889 template class TMatrixTSub <Float_t>;
1890 template class TMatrixTSparseRow <Float_t>;
1891 template class TMatrixTSparseDiag <Float_t>;
1892 template class TElementActionT <Float_t>;
1893 template class TElementPosActionT <Float_t>;
1894 
1895 template class TMatrixTRow_const <Double_t>;
1896 template class TMatrixTColumn_const <Double_t>;
1897 template class TMatrixTDiag_const <Double_t>;
1898 template class TMatrixTFlat_const <Double_t>;
1899 template class TMatrixTSub_const <Double_t>;
1900 template class TMatrixTSparseRow_const <Double_t>;
1901 template class TMatrixTSparseDiag_const<Double_t>;
1902 template class TMatrixTRow <Double_t>;
1903 template class TMatrixTColumn <Double_t>;
1904 template class TMatrixTDiag <Double_t>;
1905 template class TMatrixTFlat <Double_t>;
1906 template class TMatrixTSub <Double_t>;
1907 template class TMatrixTSparseRow <Double_t>;
1908 template class TMatrixTSparseDiag <Double_t>;
1909 template class TElementActionT <Double_t>;
1910 template class TElementPosActionT <Double_t>;
void Rank1Update(const TVectorT< Element > &vec, Element alpha=1.0)
Perform a rank 1 operation on the matrix: A += alpha * v * v^T.
Element operator()(Int_t i) const
const TMatrixTSparse< Element > * fMatrix
virtual const Element * GetMatrixArray() const =0
const TMatrixTBase< Element > * GetMatrix() const
virtual const Element * GetMatrixArray() const
virtual TMatrixTBase< Element > & InsertRow(Int_t row, Int_t col, const Element *v, Int_t n=-1)
Insert in row rown, n elements of array v at column coln.
const Element * GetPtr() const
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...
void operator+=(Element val)
Add val to every element of the matrix column.
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:132
Int_t GetInc() const
Int_t GetRowIndex() const
virtual const Int_t * GetRowIndexArray() const
const Element * fPtr
void operator*=(Element val)
Multiply every element of the matrix with val.
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:223
const TMatrixTBase< Element > * fMatrix
Int_t GetNcols() const
Definition: TMatrixTBase.h:137
Int_t GetLwb() const
Definition: TVectorT.h:79
void operator+=(Element val)
Assign val to every element of the matrix diagonal.
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
TVectorT.
Definition: TMatrixTBase.h:89
void operator*=(Element val)
Multiply every element of the sub matrix by val .
#define R__ASSERT(e)
Definition: TError.h:98
Int_t GetNrows() const
Definition: TVectorT.h:81
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:170
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
const Bool_t kFALSE
Definition: Rtypes.h:92
void operator=(Element val)
Assign val to every element of the matrix diagonal.
const Element * fDataPtr
Element operator()(Int_t i) const
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSym.h:193
void operator+=(Element val)
Add val to every non-zero (!) element of the matrix row.
void operator*=(Element val)
Multiply every element of the matrix row with val.
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
TMatrixT.
Definition: TMatrixDfwd.h:24
Int_t GetNoElements() const
Definition: TMatrixTBase.h:138
Int_t GetNrows() const
const Element * GetPtr() const
const Element * fPtr
Int_t GetColLwb() const
Definition: TMatrixTBase.h:135
const TMatrixTBase< Element > * fMatrix
void operator*=(Element val)
Multiply every element of the matrix diagonal by val.
const Element * GetPtr() const
static double p2(double t, double a, double b, double c)
const Element * fDataPtr
Int_t GetRowIndex() const
TMatrixTSym.
void Error(const char *location, const char *msgfmt,...)
Int_t GetColOff() const
void operator=(Element val)
Assign val to every element of the matrix column.
Element * GetMatrixArray()
Definition: TVectorT.h:84
void operator=(Element val)
Assign val to every element of the matrix.
TRandom2 r(17)
SVector< double, 2 > v
Definition: Dict.h:5
TMatrixTSparse.
Int_t GetNcols() const
virtual void ExtractRow(Int_t row, Int_t col, Element *v, Int_t n=-1) const
Store in array v, n matrix elements of row rown starting at column coln.
TMatrixT< Element > & ElementMult(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Multiply target by the source, element-by-element.
Definition: TMatrixT.cxx:2981
Element operator()(Int_t i) const
Bool_t IsValid() const
Definition: TVectorT.h:89
TMarker * m
Definition: textangle.C:8
const Element * fPtr
Int_t GetNoElements() const
Definition: TVectorT.h:82
static double p1(double t, double a, double b)
Int_t GetColUpb() const
Definition: TMatrixTBase.h:136
const TMatrixTBase< Element > * fMatrix
void operator+=(Element val)
Add val to every element of the sub matrix.
Linear Algebra Package.
const TMatrixTBase< Element > * GetMatrix() const
virtual const Int_t * GetColIndexArray() const
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
const TMatrixTBase< Element > * GetMatrix() const
void operator+=(Element val)
Add val to every element of the matrix.
void operator=(Element val)
Assign val to every element of the sub matrix.
void operator+=(Element val)
Add val to every element of the matrix diagonal.
double Double_t
Definition: RtypesCore.h:55
Element operator()(Int_t i) const
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
const TMatrixTBase< Element > * GetMatrix() const
const TMatrixTSparse< Element > * fMatrix
void operator=(Element val)
Assign val to every element of the matrix row.
const TMatrixTBase< Element > * fMatrix
Int_t GetRowOff() const
const TMatrixTBase< Element > * GetMatrix() const
void operator*=(Element val)
Multiply every element of the matrix row by val.
const Element * fPtr
Int_t GetColIndex() const
void operator*=(Element val)
Multiply every element of the matrix column with val.
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
static Element & NaNValue()
void operator=(Element val)
Assign val to every non-zero (!) element of the matrix row.
Int_t GetRowUpb() const
Definition: TMatrixTBase.h:133
const Bool_t kTRUE
Definition: Rtypes.h:91
const TMatrixTBase< Element > * fMatrix
void operator*=(Element val)
Assign val to every element of the matrix diagonal.
virtual TMatrixTBase< Element > & GetSub(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, TMatrixTBase< Element > &target, Option_t *option="S") const =0
void operator+=(Element val)
Add val to every element of the matrix row.
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMath.h:931
const TMatrixTBase< Element > * GetMatrix() const
const TMatrixTBase< Element > * GetMatrix() const
Bool_t IsValid() const
Definition: TMatrixTBase.h:157
Int_t GetNdiags() const
void operator=(Element val)
Assign val to every element of the matrix diagonal.
Int_t GetInc() const