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