Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TMatrixTUtils.h
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 #ifndef ROOT_TMatrixTUtils
13 #define ROOT_TMatrixTUtils
14 
15 //////////////////////////////////////////////////////////////////////////
16 // //
17 // Matrix utility classes. //
18 // //
19 // Templates of utility classes in the Linear Algebra Package. //
20 // The following classes are defined here: //
21 // //
22 // Different matrix views without copying data elements : //
23 // TMatrixTRow_const TMatrixTRow //
24 // TMatrixTColumn_const TMatrixTColumn //
25 // TMatrixTDiag_const TMatrixTDiag //
26 // TMatrixTFlat_const TMatrixTFlat //
27 // TMatrixTSub_const TMatrixTSub //
28 // TMatrixTSparseRow_const TMatrixTSparseRow //
29 // TMatrixTSparseDiag_const TMatrixTSparseDiag //
30 // //
31 // TElementActionT //
32 // TElementPosActionT //
33 // //
34 //////////////////////////////////////////////////////////////////////////
35 
36 #include "TMatrixTBase.h"
37 
38 #include <initializer_list>
39 
40 template<class Element> class TVectorT;
41 template<class Element> class TMatrixT;
42 template<class Element> class TMatrixTSym;
43 template<class Element> class TMatrixTSparse;
44 
45 //////////////////////////////////////////////////////////////////////////
46 // //
47 // TElementActionT //
48 // //
49 // A class to do a specific operation on every vector or matrix element //
50 // (regardless of it position) as the object is being traversed. //
51 // This is an abstract class. Derived classes need to implement the //
52 // action function Operation(). //
53 // //
54 //////////////////////////////////////////////////////////////////////////
55 
56 template<class Element> class TElementActionT {
57 
58 #ifndef __CINT__
59 friend class TMatrixTBase <Element>;
60 friend class TMatrixT <Element>;
61 friend class TMatrixTSym <Element>;
62 friend class TMatrixTSparse<Element>;
63 friend class TVectorT <Element>;
64 #endif
65 
66 protected:
67  virtual ~TElementActionT() { }
68  virtual void Operation(Element &element) const = 0;
69 
70 private:
71  TElementActionT& operator=(const TElementActionT<Element> &) {return *this;}
72 };
73 
74 //////////////////////////////////////////////////////////////////////////
75 // //
76 // TElementPosActionT //
77 // //
78 // A class to do a specific operation on every vector or matrix element //
79 // as the object is being traversed. This is an abstract class. //
80 // Derived classes need to implement the action function Operation(). //
81 // In the action function the location of the current element is //
82 // known (fI=row, fJ=columns). //
83 // //
84 //////////////////////////////////////////////////////////////////////////
85 
86 template<class Element> class TElementPosActionT {
87 
88 #ifndef __CINT__
89 friend class TMatrixTBase <Element>;
90 friend class TMatrixT <Element>;
91 friend class TMatrixTSym <Element>;
92 friend class TMatrixTSparse<Element>;
93 friend class TVectorT <Element>;
94 #endif
95 
96 protected:
97  mutable Int_t fI; // i position of element being passed to Operation()
98  mutable Int_t fJ; // j position of element being passed to Operation()
99  virtual ~TElementPosActionT() { }
100  virtual void Operation(Element &element) const = 0;
101 
102 private:
103  TElementPosActionT<Element>& operator=(const TElementPosActionT<Element> &) {return *this;}
104 };
105 
106 //////////////////////////////////////////////////////////////////////////
107 // //
108 // TMatrixTRow_const //
109 // //
110 // Template class represents a row of a TMatrixT/TMatrixTSym //
111 // //
112 //////////////////////////////////////////////////////////////////////////
113 
114 template<class Element> class TMatrixTRow_const {
115 
116 protected:
117  const TMatrixTBase<Element> *fMatrix; // the matrix I am a row of
118  Int_t fRowInd; // effective row index
119  Int_t fInc; // if ptr = @a[row,i], then ptr+inc = @a[row,i+1]
120  const Element *fPtr; // pointer to the a[row,0]
121 
122 public:
123  TMatrixTRow_const() { fMatrix = 0; fRowInd = 0; fInc = 0; fPtr = 0; }
124  TMatrixTRow_const(const TMatrixT <Element> &matrix,Int_t row);
125  TMatrixTRow_const(const TMatrixTSym<Element> &matrix,Int_t row);
126  TMatrixTRow_const(const TMatrixTRow_const<Element>& trc):
127  fMatrix(trc.fMatrix), fRowInd(trc.fRowInd), fInc(trc.fInc), fPtr(trc.fPtr) { }
128  TMatrixTRow_const<Element>& operator=(const TMatrixTRow_const<Element>& trc) {
129  if(this != &trc) { fMatrix=trc.fMatrix; fRowInd=trc.fRowInd; fInc=trc.fInc; fPtr=trc.fPtr; } return *this;}
130  virtual ~TMatrixTRow_const() { }
131 
132  inline const TMatrixTBase<Element> *GetMatrix () const { return fMatrix; }
133  inline Int_t GetRowIndex() const { return fRowInd; }
134  inline Int_t GetInc () const { return fInc; }
135  inline const Element *GetPtr () const { return fPtr; }
136  inline const Element &operator ()(Int_t i) const {
137  if (!fMatrix) return TMatrixTBase<Element>::NaNValue();
138  R__ASSERT(fMatrix->IsValid());
139  const Int_t acoln = i-fMatrix->GetColLwb();
140  if (acoln < fMatrix->GetNcols() && acoln >= 0)
141  return fPtr[acoln];
142  else {
143  Error("operator()","Request col(%d) outside matrix range of %d - %d",
144  i,fMatrix->GetColLwb(),fMatrix->GetColLwb()+fMatrix->GetNcols());
145  return TMatrixTBase<Element>::NaNValue();
146  }
147  }
148  inline const Element &operator [](Int_t i) const { return (*(const TMatrixTRow_const<Element> *)this)(i); }
149 
150  ClassDef(TMatrixTRow_const,0) // Template of General Matrix Row Access class
151 };
152 
153 template<class Element> class TMatrixTRow : public TMatrixTRow_const<Element> {
154 
155 public:
156  TMatrixTRow() {}
157  TMatrixTRow(TMatrixT <Element> &matrix,Int_t row);
158  TMatrixTRow(TMatrixTSym<Element> &matrix,Int_t row);
159  TMatrixTRow(const TMatrixTRow<Element> &mr);
160 
161  inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
162 
163  inline const Element &operator()(Int_t i) const {
164  if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
165  R__ASSERT(this->fMatrix->IsValid());
166  const Int_t acoln = i-this->fMatrix->GetColLwb();
167  if (acoln < this->fMatrix->GetNcols() || acoln >= 0)
168  return (this->fPtr)[acoln];
169  else {
170  Error("operator()","Request col(%d) outside matrix range of %d - %d",
171  i,this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
172  return TMatrixTBase<Element>::NaNValue();
173  }
174  }
175  inline Element &operator()(Int_t i) {
176  if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
177  R__ASSERT(this->fMatrix->IsValid());
178  const Int_t acoln = i-this->fMatrix->GetColLwb();
179  if (acoln < this->fMatrix->GetNcols() && acoln >= 0)
180  return (const_cast<Element *>(this->fPtr))[acoln];
181  else {
182  Error("operator()","Request col(%d) outside matrix range of %d - %d",
183  i,this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
184  //return (const_cast<Element *>(this->fPtr))[0];
185  return TMatrixTBase<Element>::NaNValue();
186  }
187  }
188  inline const Element &operator[](Int_t i) const { return (*(const TMatrixTRow<Element> *)this)(i); }
189  inline Element &operator[](Int_t i) { return (*( TMatrixTRow<Element> *)this)(i); }
190 
191  void Assign (Element val);
192  void operator= (std::initializer_list<Element> l);
193  void operator+=(Element val);
194  void operator*=(Element val);
195 
196  void operator=(const TMatrixTRow_const<Element> &r);
197  TMatrixTRow<Element>& operator=(const TMatrixTRow <Element> &r) { operator=((TMatrixTRow_const<Element> &)r); return *this;}
198  void operator=(const TVectorT <Element> &vec);
199 
200  void operator+=(const TMatrixTRow_const<Element> &r);
201  void operator*=(const TMatrixTRow_const<Element> &r);
202 
203  ClassDef(TMatrixTRow,0) // Template of General Matrix Row Access class
204 };
205 
206 //////////////////////////////////////////////////////////////////////////
207 // //
208 // TMatrixTColumn_const //
209 // //
210 // Template class represents a column of a TMatrixT/TMatrixTSym //
211 // //
212 //////////////////////////////////////////////////////////////////////////
213 
214 template<class Element> class TMatrixTColumn_const {
215 
216 protected:
217  const TMatrixTBase<Element> *fMatrix; // the matrix I am a column of
218  Int_t fColInd; // effective column index
219  Int_t fInc; // if ptr = @a[i,col], then ptr+inc = @a[i+1,col]
220  const Element *fPtr; // pointer to the a[0,col] column
221 
222 public:
223  TMatrixTColumn_const() { fMatrix = 0; fColInd = 0; fInc = 0; fPtr = 0; }
224  TMatrixTColumn_const(const TMatrixT <Element> &matrix,Int_t col);
225  TMatrixTColumn_const(const TMatrixTSym<Element> &matrix,Int_t col);
226  TMatrixTColumn_const(const TMatrixTColumn_const<Element>& trc):
227  fMatrix(trc.fMatrix), fColInd(trc.fColInd), fInc(trc.fInc), fPtr(trc.fPtr) { }
228  TMatrixTColumn_const<Element>& operator=(const TMatrixTColumn_const<Element>& trc) {
229  if(this != &trc) { fMatrix=trc.fMatrix; fColInd=trc.fColInd; fInc=trc.fInc; fPtr=trc.fPtr; } return *this;}
230  virtual ~TMatrixTColumn_const() { }
231 
232  inline const TMatrixTBase <Element> *GetMatrix () const { return fMatrix; }
233  inline Int_t GetColIndex() const { return fColInd; }
234  inline Int_t GetInc () const { return fInc; }
235  inline const Element *GetPtr () const { return fPtr; }
236  inline const Element &operator ()(Int_t i) const {
237  if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
238  R__ASSERT(fMatrix->IsValid());
239  const Int_t arown = i-fMatrix->GetRowLwb();
240  if (arown < fMatrix->GetNrows() && arown >= 0)
241  return fPtr[arown*fInc];
242  else {
243  Error("operator()","Request row(%d) outside matrix range of %d - %d",
244  i,fMatrix->GetRowLwb(),fMatrix->GetRowLwb()+fMatrix->GetNrows());
245  return TMatrixTBase<Element>::NaNValue();
246  }
247  }
248  inline const Element &operator [](Int_t i) const { return (*(const TMatrixTColumn_const<Element> *)this)(i); }
249 
250  ClassDef(TMatrixTColumn_const,0) // Template of General Matrix Column Access class
251 };
252 
253 template<class Element> class TMatrixTColumn : public TMatrixTColumn_const<Element> {
254 
255 public:
256  TMatrixTColumn() {}
257  TMatrixTColumn(TMatrixT <Element>&matrix,Int_t col);
258  TMatrixTColumn(TMatrixTSym<Element>&matrix,Int_t col);
259  TMatrixTColumn(const TMatrixTColumn <Element>&mc);
260 
261  inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
262 
263  inline const Element &operator()(Int_t i) const {
264  if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
265  R__ASSERT(this->fMatrix->IsValid());
266  const Int_t arown = i-this->fMatrix->GetRowLwb();
267  if (arown < this->fMatrix->GetNrows() && arown >= 0)
268  return (this->fPtr)[arown*this->fInc];
269  else {
270  Error("operator()","Request row(%d) outside matrix range of %d - %d",
271  i,this->fMatrix->GetRowLwb(),this->fMatrix->GetRowLwb()+this->fMatrix->GetNrows());
272  return TMatrixTBase<Element>::NaNValue();
273  }
274  }
275  inline Element &operator()(Int_t i) {
276  if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
277  R__ASSERT(this->fMatrix->IsValid());
278  const Int_t arown = i-this->fMatrix->GetRowLwb();
279 
280  if (arown < this->fMatrix->GetNrows() && arown >= 0)
281  return (const_cast<Element *>(this->fPtr))[arown*this->fInc];
282  else {
283  Error("operator()","Request row(%d) outside matrix range of %d - %d",
284  i,this->fMatrix->GetRowLwb(),this->fMatrix->GetRowLwb()+this->fMatrix->GetNrows());
285  return TMatrixTBase<Element>::NaNValue();
286  }
287  }
288  inline const Element &operator[](Int_t i) const { return (*(const TMatrixTColumn<Element> *)this)(i); }
289  inline Element &operator[](Int_t i) { return (*( TMatrixTColumn<Element> *)this)(i); }
290 
291  void Assign (Element val);
292  // keep it for backward compatibility (but it has been removed for TMatrixTRow)
293  void operator= (Element val) { return Assign(val); }
294  void operator= (std::initializer_list<Element> l);
295  void operator+=(Element val);
296  void operator*=(Element val);
297 
298  void operator=(const TMatrixTColumn_const<Element> &c);
299  TMatrixTColumn<Element>& operator=(const TMatrixTColumn <Element> &c) { operator=((TMatrixTColumn_const<Element> &)c); return *this;}
300  void operator=(const TVectorT <Element> &vec);
301 
302  void operator+=(const TMatrixTColumn_const<Element> &c);
303  void operator*=(const TMatrixTColumn_const<Element> &c);
304 
305  ClassDef(TMatrixTColumn,0) // Template of General Matrix Column Access class
306 };
307 
308 //////////////////////////////////////////////////////////////////////////
309 // //
310 // TMatrixTDiag_const //
311 // //
312 // Template class represents the diagonal of a TMatrixT/TMatrixTSym //
313 // //
314 //////////////////////////////////////////////////////////////////////////
315 
316 template<class Element> class TMatrixTDiag_const {
317 
318 protected:
319  const TMatrixTBase<Element> *fMatrix; // the matrix I am the diagonal of
320  Int_t fInc; // if ptr=@a[i,i], then ptr+inc = @a[i+1,i+1]
321  Int_t fNdiag; // number of diag elems, min(nrows,ncols)
322  const Element *fPtr; // pointer to the a[0,0]
323 
324 public:
325  TMatrixTDiag_const() { fMatrix = 0; fInc = 0; fNdiag = 0; fPtr = 0; }
326  TMatrixTDiag_const(const TMatrixT <Element> &matrix);
327  TMatrixTDiag_const(const TMatrixTSym<Element> &matrix);
328  TMatrixTDiag_const(const TMatrixTDiag_const<Element>& trc):
329  fMatrix(trc.fMatrix), fInc(trc.fInc), fNdiag(trc.fNdiag), fPtr(trc.fPtr) { }
330  TMatrixTDiag_const<Element>& operator=(const TMatrixTDiag_const<Element>& trc) {
331  if(this != &trc) { fMatrix=trc.fMatrix; fInc=trc.fInc; fNdiag=trc.fNdiag; fPtr=trc.fPtr; } return *this;}
332  virtual ~TMatrixTDiag_const() { }
333 
334  inline const TMatrixTBase<Element> *GetMatrix() const { return fMatrix; }
335  inline const Element *GetPtr () const { return fPtr; }
336  inline Int_t GetInc () const { return fInc; }
337  inline const Element &operator ()(Int_t i) const {
338  R__ASSERT(fMatrix->IsValid());
339  if (i < fNdiag && i >= 0)
340  return fPtr[i*fInc];
341  else {
342  Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,fNdiag);
343  return TMatrixTBase<Element>::NaNValue();
344  }
345  }
346  inline const Element &operator [](Int_t i) const { return (*(const TMatrixTDiag_const<Element> *)this)(i); }
347 
348  Int_t GetNdiags() const { return fNdiag; }
349 
350  ClassDef(TMatrixTDiag_const,0) // Template of General Matrix Diagonal Access class
351 };
352 
353 template<class Element> class TMatrixTDiag : public TMatrixTDiag_const<Element> {
354 
355 public:
356  TMatrixTDiag() {}
357  TMatrixTDiag(TMatrixT <Element>&matrix);
358  TMatrixTDiag(TMatrixTSym<Element>&matrix);
359  TMatrixTDiag(const TMatrixTDiag<Element> &md);
360 
361  inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
362 
363  inline const Element &operator()(Int_t i) const {
364  R__ASSERT(this->fMatrix->IsValid());
365  if (i < this->fNdiag && i >= 0)
366  return (this->fPtr)[i*this->fInc];
367  else {
368  Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,this->fNdiag);
369  return TMatrixTBase<Element>::NaNValue();
370  }
371  }
372  inline Element &operator()(Int_t i) {
373  R__ASSERT(this->fMatrix->IsValid());
374  if (i < this->fNdiag && i >= 0)
375  return (const_cast<Element *>(this->fPtr))[i*this->fInc];
376  else {
377  Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,this->fNdiag);
378  return (const_cast<Element *>(this->fPtr))[0];
379  }
380  }
381  inline const Element &operator[](Int_t i) const { return (*(const TMatrixTDiag<Element> *)this)(i); }
382  inline Element &operator[](Int_t i) { return (*( TMatrixTDiag *)this)(i); }
383 
384  void operator= (Element val);
385  void operator+=(Element val);
386  void operator*=(Element val);
387 
388  void operator=(const TMatrixTDiag_const<Element> &d);
389  TMatrixTDiag<Element>& operator=(const TMatrixTDiag <Element> &d) { operator=((TMatrixTDiag_const<Element> &)d); return *this;}
390  void operator=(const TVectorT <Element> &vec);
391 
392  void operator+=(const TMatrixTDiag_const<Element> &d);
393  void operator*=(const TMatrixTDiag_const<Element> &d);
394 
395  ClassDef(TMatrixTDiag,0) // Template of General Matrix Diagonal Access class
396 };
397 
398 //////////////////////////////////////////////////////////////////////////
399 // //
400 // TMatrixTFlat_const //
401 // //
402 // Template class represents a flat TMatrixT/TMatrixTSym //
403 // //
404 //////////////////////////////////////////////////////////////////////////
405 
406 template<class Element> class TMatrixTFlat_const {
407 
408 protected:
409  const TMatrixTBase<Element> *fMatrix; // the matrix I am the diagonal of
410  Int_t fNelems; //
411  const Element *fPtr; // pointer to the a[0,0]
412 
413 public:
414  TMatrixTFlat_const() { fMatrix = 0; fNelems = 0; fPtr = 0; }
415  TMatrixTFlat_const(const TMatrixT <Element> &matrix);
416  TMatrixTFlat_const(const TMatrixTSym<Element> &matrix);
417  TMatrixTFlat_const(const TMatrixTFlat_const<Element>& trc):
418  fMatrix(trc.fMatrix), fNelems(trc.fNelems), fPtr(trc.fPtr) { }
419  TMatrixTFlat_const<Element>& operator=(const TMatrixTFlat_const<Element>& trc) {
420  if(this != &trc) { fMatrix=trc.fMatrix; fNelems=trc.fNelems; fPtr=trc.fPtr; } return *this;}
421  virtual ~TMatrixTFlat_const() { }
422 
423  inline const TMatrixTBase<Element> *GetMatrix() const { return fMatrix; }
424  inline const Element *GetPtr () const { return fPtr; }
425  inline const Element &operator ()(Int_t i) const {
426  R__ASSERT(fMatrix->IsValid());
427  if (i < fNelems && i >= 0)
428  return fPtr[i];
429  else {
430  Error("operator()","Request element(%d) outside matrix range of 0 - %d",i,fNelems);
431  return TMatrixTBase<Element>::NaNValue();
432  }
433  }
434  inline const Element &operator [](Int_t i) const { return (*(const TMatrixTFlat_const<Element> *)this)(i); }
435 
436  ClassDef(TMatrixTFlat_const,0) // Template of General Matrix Flat Representation class
437 };
438 
439 template<class Element> class TMatrixTFlat : public TMatrixTFlat_const<Element> {
440 
441 public:
442  TMatrixTFlat() {}
443  TMatrixTFlat(TMatrixT <Element> &matrix);
444  TMatrixTFlat(TMatrixTSym<Element> &matrix);
445  TMatrixTFlat(const TMatrixTFlat<Element> &mf);
446 
447  inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
448 
449  inline const Element &operator()(Int_t i) const {
450  R__ASSERT(this->fMatrix->IsValid());
451  if (i < this->fNelems && i >= 0)
452  return (this->fPtr)[i];
453  else {
454  Error("operator()","Request element(%d) outside matrix range of 0 - %d",i,this->fNelems);
455  return TMatrixTBase<Element>::NaNValue();
456  }
457  }
458  inline Element &operator()(Int_t i) {
459  R__ASSERT(this->fMatrix->IsValid());
460  if (i < this->fNelems && i >= 0)
461  return (const_cast<Element *>(this->fPtr))[i];
462  else {
463  Error("operator()","Request element(%d) outside matrix range of 0 - %d",i,this->fNelems);
464  return TMatrixTBase<Element>::NaNValue();
465  }
466  }
467  inline const Element &operator[](Int_t i) const { return (*(const TMatrixTFlat<Element> *)this)(i); }
468  inline Element &operator[](Int_t i) { return (*( TMatrixTFlat<Element> *)this)(i); }
469 
470  void operator= (Element val);
471  void operator+=(Element val);
472  void operator*=(Element val);
473 
474  void operator=(const TMatrixTFlat_const<Element> &f);
475  TMatrixTFlat<Element>& operator=(const TMatrixTFlat <Element> &f) { operator=((TMatrixTFlat_const<Element> &)f); return *this;}
476  void operator=(const TVectorT <Element> &vec);
477 
478  void operator+=(const TMatrixTFlat_const<Element> &f);
479  void operator*=(const TMatrixTFlat_const<Element> &f);
480 
481  ClassDef(TMatrixTFlat,0) // Template of General Matrix Flat Representation class
482 };
483 
484 //////////////////////////////////////////////////////////////////////////
485 // //
486 // TMatrixTSub_const //
487 // //
488 // Template class represents a sub matrix of TMatrixT/TMatrixTSym //
489 // //
490 //////////////////////////////////////////////////////////////////////////
491 
492 template<class Element> class TMatrixTSub_const {
493 
494 protected:
495  const TMatrixTBase<Element> *fMatrix; // the matrix I am a submatrix of
496  Int_t fRowOff; //
497  Int_t fColOff; //
498  Int_t fNrowsSub; //
499  Int_t fNcolsSub; //
500 
501 public:
502  TMatrixTSub_const() { fRowOff = fColOff = fNrowsSub = fNcolsSub = 0; fMatrix = 0; }
503  TMatrixTSub_const(const TMatrixT <Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
504  TMatrixTSub_const(const TMatrixTSym<Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
505  virtual ~TMatrixTSub_const() { }
506 
507  inline const TMatrixTBase<Element> *GetMatrix() const { return fMatrix; }
508  inline Int_t GetRowOff() const { return fRowOff; }
509  inline Int_t GetColOff() const { return fColOff; }
510  inline Int_t GetNrows () const { return fNrowsSub; }
511  inline Int_t GetNcols () const { return fNcolsSub; }
512  inline const Element &operator ()(Int_t rown,Int_t coln) const {
513  R__ASSERT(fMatrix->IsValid());
514 
515  const Element *ptr = fMatrix->GetMatrixArray();
516  if (rown >= fNrowsSub || rown < 0) {
517  Error("operator()","Request row(%d) outside matrix range of 0 - %d",rown,fNrowsSub);
518  return TMatrixTBase<Element>::NaNValue();
519  }
520  if (coln >= fNcolsSub || coln < 0) {
521  Error("operator()","Request column(%d) outside matrix range of 0 - %d",coln,fNcolsSub);
522  return TMatrixTBase<Element>::NaNValue();
523  }
524  const Int_t index = (rown+fRowOff)*fMatrix->GetNcols()+coln+fColOff;
525  return ptr[index];
526  }
527 
528  ClassDef(TMatrixTSub_const,0) // Template of Sub Matrix Access class
529 };
530 
531 template<class Element> class TMatrixTSub : public TMatrixTSub_const<Element> {
532 
533 public:
534 
535  enum {kWorkMax = 100};
536 
537  TMatrixTSub() {}
538  TMatrixTSub(TMatrixT <Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
539  TMatrixTSub(TMatrixTSym<Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
540  TMatrixTSub(const TMatrixTSub<Element> &ms);
541 
542  inline Element &operator()(Int_t rown,Int_t coln) {
543  R__ASSERT(this->fMatrix->IsValid());
544 
545  const Element *ptr = this->fMatrix->GetMatrixArray();
546  if (rown >= this->fNrowsSub || rown < 0) {
547  Error("operator()","Request row(%d) outside matrix range of 0 - %d",rown,this->fNrowsSub);
548  return TMatrixTBase<Element>::NaNValue();
549  }
550  if (coln >= this->fNcolsSub || coln < 0) {
551  Error("operator()","Request column(%d) outside matrix range of 0 - %d",coln,this->fNcolsSub);
552  return TMatrixTBase<Element>::NaNValue();
553  }
554  const Int_t index = (rown+this->fRowOff)*this->fMatrix->GetNcols()+coln+this->fColOff;
555  return (const_cast<Element *>(ptr))[index];
556  }
557 
558  void Rank1Update(const TVectorT<Element> &vec,Element alpha=1.0);
559 
560  void operator= (Element val);
561  void operator+=(Element val);
562  void operator*=(Element val);
563 
564  void operator=(const TMatrixTSub_const<Element> &s);
565  TMatrixTSub<Element>& operator=(const TMatrixTSub <Element> &s) { operator=((TMatrixTSub_const<Element> &)s); return *this;}
566  void operator=(const TMatrixTBase <Element> &m);
567 
568  void operator+=(const TMatrixTSub_const<Element> &s);
569  void operator*=(const TMatrixTSub_const<Element> &s);
570  void operator+=(const TMatrixTBase <Element> &m);
571  void operator*=(const TMatrixT <Element> &m);
572  void operator*=(const TMatrixTSym <Element> &m);
573 
574  ClassDef(TMatrixTSub,0) // Template of Sub Matrix Access class
575 };
576 
577 //////////////////////////////////////////////////////////////////////////
578 // //
579 // TMatrixTSparseRow_const //
580 // //
581 // Template class represents a row of TMatrixTSparse //
582 // //
583 //////////////////////////////////////////////////////////////////////////
584 
585 template<class Element> class TMatrixTSparseRow_const {
586 
587 protected:
588  const TMatrixTSparse<Element> *fMatrix; // the matrix I am a row of
589  Int_t fRowInd; // effective row index
590  Int_t fNindex; // index range
591  const Int_t *fColPtr; // column index pointer
592  const Element *fDataPtr; // data pointer
593 
594 public:
595  TMatrixTSparseRow_const() { fMatrix = 0; fRowInd = 0; fNindex = 0; fColPtr = 0; fDataPtr = 0; }
596  TMatrixTSparseRow_const(const TMatrixTSparse<Element> &matrix,Int_t row);
597  TMatrixTSparseRow_const(const TMatrixTSparseRow_const<Element>& trc):
598  fMatrix(trc.fMatrix), fRowInd(trc.fRowInd), fNindex(trc.fNindex), fColPtr(trc.fColPtr), fDataPtr(trc.fDataPtr) { }
599  TMatrixTSparseRow_const<Element>& operator=(const TMatrixTSparseRow_const<Element>& trc) {
600  if(this != &trc) { fMatrix=trc.fMatrix; fRowInd=trc.fRowInd; fNindex=trc.fNindex; fColPtr=trc.fColPtr; fDataPtr=trc.fDataPtr; } return *this;}
601  virtual ~TMatrixTSparseRow_const() { }
602 
603  inline const TMatrixTBase<Element> *GetMatrix () const { return fMatrix; }
604  inline const Element *GetDataPtr () const { return fDataPtr; }
605  inline const Int_t *GetColPtr () const { return fColPtr; }
606  inline Int_t GetRowIndex() const { return fRowInd; }
607  inline Int_t GetNindex () const { return fNindex; }
608 
609  Element operator()(Int_t i) const;
610  inline Element operator[](Int_t i) const { return (*(const TMatrixTSparseRow_const<Element> *)this)(i); }
611 
612  ClassDef(TMatrixTSparseRow_const,0) // Template of Sparse Matrix Row Access class
613 };
614 
615 template<class Element> class TMatrixTSparseRow : public TMatrixTSparseRow_const<Element> {
616 
617 public:
618  TMatrixTSparseRow() {}
619  TMatrixTSparseRow(TMatrixTSparse<Element> &matrix,Int_t row);
620  TMatrixTSparseRow(const TMatrixTSparseRow<Element> &mr);
621 
622  inline Element *GetDataPtr() const { return const_cast<Element *>(this->fDataPtr); }
623 
624  Element operator()(Int_t i) const;
625  Element &operator()(Int_t i);
626  inline Element operator[](Int_t i) const { return (*(const TMatrixTSparseRow<Element> *)this)(i); }
627  inline Element &operator[](Int_t i) { return (*(TMatrixTSparseRow<Element> *)this)(i); }
628 
629  void operator= (Element val);
630  void operator+=(Element val);
631  void operator*=(Element val);
632 
633  void operator=(const TMatrixTSparseRow_const<Element> &r);
634  TMatrixTSparseRow<Element>& operator=(const TMatrixTSparseRow <Element> &r) { operator=((TMatrixTSparseRow_const<Element> &)r); return *this;}
635  void operator=(const TVectorT <Element> &vec);
636 
637  void operator+=(const TMatrixTSparseRow_const<Element> &r);
638  void operator*=(const TMatrixTSparseRow_const<Element> &r);
639 
640  ClassDef(TMatrixTSparseRow,0) // Template of Sparse Matrix Row Access class
641 };
642 
643 //////////////////////////////////////////////////////////////////////////
644 // //
645 // TMatrixTSparseDiag_const //
646 // //
647 // Template class represents the diagonal of TMatrixTSparse //
648 // //
649 //////////////////////////////////////////////////////////////////////////
650 
651 template<class Element> class TMatrixTSparseDiag_const {
652 
653 protected:
654  const TMatrixTSparse<Element> *fMatrix; // the matrix I am the diagonal of
655  Int_t fNdiag; // number of diag elems, min(nrows,ncols)
656  const Element *fDataPtr; // data pointer
657 
658 public:
659  TMatrixTSparseDiag_const() { fMatrix = 0; fNdiag = 0; fDataPtr = 0; }
660  TMatrixTSparseDiag_const(const TMatrixTSparse<Element> &matrix);
661  TMatrixTSparseDiag_const(const TMatrixTSparseDiag_const<Element>& trc):
662  fMatrix(trc.fMatrix), fNdiag(trc.fNdiag), fDataPtr(trc.fDataPtr) { }
663  TMatrixTSparseDiag_const<Element>& operator=(const TMatrixTSparseDiag_const<Element>& trc) {
664  if(this != &trc) { fMatrix=trc.fMatrix; fNdiag=trc.fNdiag; fDataPtr=trc.fDataPtr; } return *this;}
665  virtual ~TMatrixTSparseDiag_const() { }
666 
667  inline const TMatrixTBase<Element> *GetMatrix () const { return fMatrix; }
668  inline const Element *GetDataPtr() const { return fDataPtr; }
669  inline Int_t GetNdiags () const { return fNdiag; }
670 
671  Element operator ()(Int_t i) const;
672  inline Element operator [](Int_t i) const { return (*(const TMatrixTSparseRow_const<Element> *)this)(i); }
673 
674  ClassDef(TMatrixTSparseDiag_const,0) // Template of Sparse Matrix Diagonal Access class
675 };
676 
677 template<class Element> class TMatrixTSparseDiag : public TMatrixTSparseDiag_const<Element> {
678 
679 public:
680  TMatrixTSparseDiag() {}
681  TMatrixTSparseDiag(TMatrixTSparse<Element> &matrix);
682  TMatrixTSparseDiag(const TMatrixTSparseDiag<Element> &md);
683 
684  inline Element *GetDataPtr() const { return const_cast<Element *>(this->fDataPtr); }
685 
686  Element operator()(Int_t i) const;
687  Element &operator()(Int_t i);
688  inline Element operator[](Int_t i) const { return (*(const TMatrixTSparseDiag<Element> *)this)(i); }
689  inline Element &operator[](Int_t i) { return (*(TMatrixTSparseDiag<Element> *)this)(i); }
690 
691  void operator= (Element val);
692  void operator+=(Element val);
693  void operator*=(Element val);
694 
695  void operator=(const TMatrixTSparseDiag_const<Element> &d);
696  TMatrixTSparseDiag<Element>& operator=(const TMatrixTSparseDiag <Element> &d) { operator=((TMatrixTSparseDiag_const<Element> &)d); return *this;}
697  void operator=(const TVectorT <Element> &vec);
698 
699  void operator+=(const TMatrixTSparseDiag_const<Element> &d);
700  void operator*=(const TMatrixTSparseDiag_const<Element> &d);
701 
702  ClassDef(TMatrixTSparseDiag,0) // Template of Sparse Matrix Diagonal Access class
703 };
704 
705 Double_t Drand(Double_t &ix);
706 #endif