ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
TMatrixTSparse.cxx
Go to the documentation of this file.
1 // @(#)root/matrix:$Id$
2 // Authors: Fons Rademakers, Eddy Offermann Feb 2004
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 //////////////////////////////////////////////////////////////////////////
13 // //
14 // TMatrixTSparse //
15 // //
16 // Template class of a general sparse matrix in the Harwell-Boeing //
17 // format //
18 // //
19 // Besides the usual shape/size decsriptors of a matrix like fNrows, //
20 // fRowLwb,fNcols and fColLwb, we also store a row index, fRowIndex and //
21 // column index, fColIndex only for those elements unequal zero: //
22 // //
23 // fRowIndex[0,..,fNrows]: Stores for each row the index range of //
24 // the elements in the data and column array //
25 // fColIndex[0,..,fNelems-1]: Stores the column number for each data //
26 // element != 0 //
27 // //
28 // As an example how to access all sparse data elements: //
29 // //
30 // for (Int_t irow = 0; irow < this->fNrows; irow++) { //
31 // const Int_t sIndex = fRowIndex[irow]; //
32 // const Int_t eIndex = fRowIndex[irow+1]; //
33 // for (Int_t index = sIndex; index < eIndex; index++) { //
34 // const Int_t icol = fColIndex[index]; //
35 // const Element data = fElements[index]; //
36 // printf("data(%d,%d) = %.4e\n",irow+this->fRowLwb,icol+ //
37 // this->fColLwb,data); //
38 // } //
39 // } //
40 // //
41 // When checking whether sparse matrices are compatible (like in an //
42 // assigment !), not only the shape parameters are compared but also //
43 // the sparse structure through fRowIndex and fColIndex . //
44 // //
45 // Several methods exist to fill a sparse matrix with data entries. //
46 // Most are the same like for dense matrices but some care has to be //
47 // taken with regard to performance. In the constructor, always the //
48 // shape of the matrix has to be specified in some form . Data can be //
49 // entered through the following methods : //
50 // 1. constructor //
51 // TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t dol_lwb, //
52 // Int_t col_upb,Int_t nr_nonzeros, //
53 // Int_t *row, Int_t *col,Element *data); //
54 // It uses SetMatrixArray(..), see below //
55 // 2. copy constructors //
56 // 3. SetMatrixArray(Int_t nr,Int_t *irow,Int_t *icol,Element *data) //
57 // where it is expected that the irow,icol and data array contain //
58 // nr entries . Only the entries with non-zero data[i] value are //
59 // inserted. Be aware that the input data array will be modified //
60 // inside the routine for doing the necessary sorting of indices ! //
61 // 4. TMatrixTSparse a(n,m); for(....) { a(i,j) = .... //
62 // This is a very flexible method but expensive : //
63 // - if no entry for slot (i,j) is found in the sparse index table //
64 // it will be entered, which involves some memory management ! //
65 // - before invoking this method in a loop it is smart to first //
66 // set the index table through a call to SetSparseIndex(..) //
67 // 5. SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixTBase &source) //
68 // the matrix to be inserted at position (row_lwb,col_lwb) can be //
69 // both dense or sparse . //
70 // //
71 //////////////////////////////////////////////////////////////////////////
72 
73 #include "TMatrixTSparse.h"
74 #include "TBuffer.h"
75 #include "TMatrixT.h"
76 #include "TMath.h"
77 
78 templateClassImp(TMatrixTSparse)
79 
80 
81 ////////////////////////////////////////////////////////////////////////////////
82 /// Space is allocated for row/column indices and data, but the sparse structure
83 /// information has still to be set !
84 
85 template<class Element>
86 TMatrixTSparse<Element>::TMatrixTSparse(Int_t no_rows,Int_t no_cols)
87 {
88  Allocate(no_rows,no_cols,0,0,1);
89 }
90 
91 ////////////////////////////////////////////////////////////////////////////////
92 /// Space is allocated for row/column indices and data, but the sparse structure
93 /// information has still to be set !
94 
95 template<class Element>
96 TMatrixTSparse<Element>::TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb)
97 {
98  Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb,1);
99 }
100 
101 ////////////////////////////////////////////////////////////////////////////////
102 /// Space is allocated for row/column indices and data. Sparse row/column index
103 /// structure together with data is coming from the arrays, row, col and data, resp .
104 
105 template<class Element>
106 TMatrixTSparse<Element>::TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
107  Int_t nr,Int_t *row, Int_t *col,Element *data)
108 {
109  const Int_t irowmin = TMath::LocMin(nr,row);
110  const Int_t irowmax = TMath::LocMax(nr,row);
111  const Int_t icolmin = TMath::LocMin(nr,col);
112  const Int_t icolmax = TMath::LocMax(nr,col);
113 
114  if (row[irowmin] < row_lwb || row[irowmax] > row_upb) {
115  Error("TMatrixTSparse","Inconsistency between row index and its range");
116  if (row[irowmin] < row_lwb) {
117  Info("TMatrixTSparse","row index lower bound adjusted to %d",row[irowmin]);
118  row_lwb = row[irowmin];
119  }
120  if (row[irowmax] > row_upb) {
121  Info("TMatrixTSparse","row index upper bound adjusted to %d",row[irowmax]);
122  col_lwb = col[irowmax];
123  }
124  }
125  if (col[icolmin] < col_lwb || col[icolmax] > col_upb) {
126  Error("TMatrixTSparse","Inconsistency between column index and its range");
127  if (col[icolmin] < col_lwb) {
128  Info("TMatrixTSparse","column index lower bound adjusted to %d",col[icolmin]);
129  col_lwb = col[icolmin];
130  }
131  if (col[icolmax] > col_upb) {
132  Info("TMatrixTSparse","column index upper bound adjusted to %d",col[icolmax]);
133  col_upb = col[icolmax];
134  }
135  }
136 
137  Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb,1,nr);
138 
139  SetMatrixArray(nr,row,col,data);
140 }
141 
142 ////////////////////////////////////////////////////////////////////////////////
143 
144 template<class Element>
145 TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixTSparse<Element> &another) : TMatrixTBase<Element>(another)
146 {
147  Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb(),1,
148  another.GetNoElements());
149  memcpy(fRowIndex,another.GetRowIndexArray(),this->fNrowIndex*sizeof(Int_t));
150  memcpy(fColIndex,another.GetColIndexArray(),this->fNelems*sizeof(Int_t));
151 
152  *this = another;
153 }
154 
155 ////////////////////////////////////////////////////////////////////////////////
156 
157 template<class Element>
158 TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixT<Element> &another) : TMatrixTBase<Element>(another)
159 {
160  const Int_t nr_nonzeros = another.NonZeros();
161  Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb(),1,nr_nonzeros);
162  SetSparseIndex(another);
163  *this = another;
164 }
165 
166 ////////////////////////////////////////////////////////////////////////////////
167 /// Create a matrix applying a specific operation to the prototype.
168 /// Supported operations are: kZero, kUnit, kTransposed and kAtA
169 
170 template<class Element>
171 TMatrixTSparse<Element>::TMatrixTSparse(EMatrixCreatorsOp1 op,const TMatrixTSparse<Element> &prototype)
172 {
173  R__ASSERT(prototype.IsValid());
174 
175  Int_t nr_nonzeros = 0;
176 
177  switch(op) {
178  case kZero:
179  {
180  Allocate(prototype.GetNrows(),prototype.GetNcols(),
181  prototype.GetRowLwb(),prototype.GetColLwb(),1,nr_nonzeros);
182  break;
183  }
184  case kUnit:
185  {
186  const Int_t nrows = prototype.GetNrows();
187  const Int_t ncols = prototype.GetNcols();
188  const Int_t rowLwb = prototype.GetRowLwb();
189  const Int_t colLwb = prototype.GetColLwb();
190  for (Int_t i = rowLwb; i <= rowLwb+nrows-1; i++)
191  for (Int_t j = colLwb; j <= colLwb+ncols-1; j++)
192  if (i==j) nr_nonzeros++;
193  Allocate(nrows,ncols,rowLwb,colLwb,1,nr_nonzeros);
194  UnitMatrix();
195  break;
196  }
197  case kTransposed:
198  {
199  Allocate(prototype.GetNcols(), prototype.GetNrows(),
200  prototype.GetColLwb(),prototype.GetRowLwb(),1,prototype.GetNoElements());
201  Transpose(prototype);
202  break;
203  }
204  case kAtA:
205  {
206  const TMatrixTSparse<Element> at(TMatrixTSparse<Element>::kTransposed,prototype);
207  AMultBt(at,at,1);
208  break;
209  }
210 
211  default:
212  Error("TMatrixTSparse(EMatrixCreatorOp1)","operation %d not yet implemented", op);
213  }
214 }
215 
216 ////////////////////////////////////////////////////////////////////////////////
217 /// Create a matrix applying a specific operation to two prototypes.
218 /// Supported operations are: kMult (a*b), kMultTranspose (a*b'), kPlus (a+b), kMinus (a-b)
219 
220 template<class Element>
221 TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b)
222 {
223  R__ASSERT(a.IsValid());
224  R__ASSERT(b.IsValid());
225 
226  switch(op) {
227  case kMult:
228  AMultB(a,b,1);
229  break;
230 
231  case kMultTranspose:
232  AMultBt(a,b,1);
233  break;
234 
235  case kPlus:
236  APlusB(a,b,1);
237  break;
238 
239  case kMinus:
240  AMinusB(a,b,1);
241  break;
242 
243  default:
244  Error("TMatrixTSparse(EMatrixCreatorOp2)", "operation %d not yet implemented",op);
245  }
246 }
247 
248 ////////////////////////////////////////////////////////////////////////////////
249 /// Create a matrix applying a specific operation to two prototypes.
250 /// Supported operations are: kMult (a*b), kMultTranspose (a*b'), kPlus (a+b), kMinus (a-b)
251 
252 template<class Element>
253 TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixT<Element> &b)
254 {
255  R__ASSERT(a.IsValid());
256  R__ASSERT(b.IsValid());
257 
258  switch(op) {
259  case kMult:
260  AMultB(a,b,1);
261  break;
262 
263  case kMultTranspose:
264  AMultBt(a,b,1);
265  break;
266 
267  case kPlus:
268  APlusB(a,b,1);
269  break;
270 
271  case kMinus:
272  AMinusB(a,b,1);
273  break;
274 
275  default:
276  Error("TMatrixTSparse(EMatrixCreatorOp2)", "operation %d not yet implemented",op);
277  }
278 }
279 
280 ////////////////////////////////////////////////////////////////////////////////
281 /// Create a matrix applying a specific operation to two prototypes.
282 /// Supported operations are: kMult (a*b), kMultTranspose (a*b'), kPlus (a+b), kMinus (a-b)
283 
284 template<class Element>
285 TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixT<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b)
286 {
287  R__ASSERT(a.IsValid());
288  R__ASSERT(b.IsValid());
289 
290  switch(op) {
291  case kMult:
292  AMultB(a,b,1);
293  break;
294 
295  case kMultTranspose:
296  AMultBt(a,b,1);
297  break;
298 
299  case kPlus:
300  APlusB(a,b,1);
301  break;
302 
303  case kMinus:
304  AMinusB(a,b,1);
305  break;
306 
307  default:
308  Error("TMatrixTSparse(EMatrixCreatorOp2)", "operation %d not yet implemented",op);
309  }
310 }
311 
312 ////////////////////////////////////////////////////////////////////////////////
313 /// Allocate new matrix. Arguments are number of rows, columns, row lowerbound (0 default)
314 /// and column lowerbound (0 default), 0 initialization flag and number of non-zero
315 /// elements (only relevant for sparse format).
316 
317 template<class Element>
318 void TMatrixTSparse<Element>::Allocate(Int_t no_rows,Int_t no_cols,Int_t row_lwb,Int_t col_lwb,
319  Int_t init,Int_t nr_nonzeros)
320 {
321  if ( (nr_nonzeros > 0 && (no_rows == 0 || no_cols == 0)) ||
322  (no_rows < 0 || no_cols < 0 || nr_nonzeros < 0) )
323  {
324  Error("Allocate","no_rows=%d no_cols=%d non_zeros=%d",no_rows,no_cols,nr_nonzeros);
325  this->Invalidate();
326  return;
327  }
328 
329  this->MakeValid();
330  this->fNrows = no_rows;
331  this->fNcols = no_cols;
332  this->fRowLwb = row_lwb;
333  this->fColLwb = col_lwb;
334  this->fNrowIndex = this->fNrows+1;
335  this->fNelems = nr_nonzeros;
336  this->fIsOwner = kTRUE;
337  this->fTol = std::numeric_limits<Element>::epsilon();
338 
339  fRowIndex = new Int_t[this->fNrowIndex];
340  if (init)
341  memset(fRowIndex,0,this->fNrowIndex*sizeof(Int_t));
342 
343  if (this->fNelems > 0) {
344  fElements = new Element[this->fNelems];
345  fColIndex = new Int_t [this->fNelems];
346  if (init) {
347  memset(fElements,0,this->fNelems*sizeof(Element));
348  memset(fColIndex,0,this->fNelems*sizeof(Int_t));
349  }
350  } else {
351  fElements = 0;
352  fColIndex = 0;
353  }
354 }
355 
356 ////////////////////////////////////////////////////////////////////////////////
357 /// Insert in row rown, n elements of array v at column coln
358 
359 template<class Element>
360 TMatrixTBase<Element> &TMatrixTSparse<Element>::InsertRow(Int_t rown,Int_t coln,const Element *v,Int_t n)
361 {
362  const Int_t arown = rown-this->fRowLwb;
363  const Int_t acoln = coln-this->fColLwb;
364  const Int_t nr = (n > 0) ? n : this->fNcols;
365 
366  if (gMatrixCheck) {
367  if (arown >= this->fNrows || arown < 0) {
368  Error("InsertRow","row %d out of matrix range",rown);
369  return *this;
370  }
371 
372  if (acoln >= this->fNcols || acoln < 0) {
373  Error("InsertRow","column %d out of matrix range",coln);
374  return *this;
375  }
376 
377  if (acoln+nr > this->fNcols || nr < 0) {
378  Error("InsertRow","row length %d out of range",nr);
379  return *this;
380  }
381  }
382 
383  const Int_t sIndex = fRowIndex[arown];
384  const Int_t eIndex = fRowIndex[arown+1];
385 
386  // check first how many slots are available from [acoln,..,acoln+nr-1]
387  // also note lIndex and rIndex so that [sIndex..lIndex] and [rIndex..eIndex-1]
388  // contain the row entries except for the region to be inserted
389 
390  Int_t nslots = 0;
391  Int_t lIndex = sIndex-1;
392  Int_t rIndex = sIndex-1;
393  Int_t index;
394  for (index = sIndex; index < eIndex; index++) {
395  const Int_t icol = fColIndex[index];
396  rIndex++;
397  if (icol >= acoln+nr) break;
398  if (icol >= acoln) nslots++;
399  else lIndex++;
400  }
401 
402  const Int_t nelems_old = this->fNelems;
403  const Int_t *colIndex_old = fColIndex;
404  const Element *elements_old = fElements;
405 
406  const Int_t ndiff = nr-nslots;
407  this->fNelems += ndiff;
408  fColIndex = new Int_t[this->fNelems];
409  fElements = new Element[this->fNelems];
410 
411  for (Int_t irow = arown+1; irow < this->fNrows+1; irow++)
412  fRowIndex[irow] += ndiff;
413 
414  if (lIndex+1 > 0) {
415  memmove(fColIndex,colIndex_old,(lIndex+1)*sizeof(Int_t));
416  memmove(fElements,elements_old,(lIndex+1)*sizeof(Element));
417  }
418 
419  if (nelems_old > 0 && nelems_old-rIndex > 0) {
420  memmove(fColIndex+rIndex+ndiff,colIndex_old+rIndex,(nelems_old-rIndex)*sizeof(Int_t));
421  memmove(fElements+rIndex+ndiff,elements_old+rIndex,(nelems_old-rIndex)*sizeof(Element));
422  }
423 
424  index = lIndex+1;
425  for (Int_t i = 0; i < nr; i++) {
426  fColIndex[index] = acoln+i;
427  fElements[index] = v[i];
428  index++;
429  }
430 
431  if (colIndex_old) delete [] (Int_t*) colIndex_old;
432  if (elements_old) delete [] (Element*) elements_old;
433 
434  R__ASSERT(this->fNelems == fRowIndex[this->fNrowIndex-1]);
435 
436  return *this;
437 }
438 
439 ////////////////////////////////////////////////////////////////////////////////
440 /// Store in array v, n matrix elements of row rown starting at column coln
441 
442 template<class Element>
443 void TMatrixTSparse<Element>::ExtractRow(Int_t rown, Int_t coln, Element *v,Int_t n) const
444 {
445  const Int_t arown = rown-this->fRowLwb;
446  const Int_t acoln = coln-this->fColLwb;
447  const Int_t nr = (n > 0) ? n : this->fNcols;
448 
449  if (gMatrixCheck) {
450  if (arown >= this->fNrows || arown < 0) {
451  Error("ExtractRow","row %d out of matrix range",rown);
452  return;
453  }
454 
455  if (acoln >= this->fNcols || acoln < 0) {
456  Error("ExtractRow","column %d out of matrix range",coln);
457  return;
458  }
459 
460  if (acoln+nr > this->fNcols || nr < 0) {
461  Error("ExtractRow","row length %d out of range",nr);
462  return;
463  }
464  }
465 
466  const Int_t sIndex = fRowIndex[arown];
467  const Int_t eIndex = fRowIndex[arown+1];
468 
469  memset(v,0,nr*sizeof(Element));
470  const Int_t * const pColIndex = GetColIndexArray();
471  const Element * const pData = GetMatrixArray();
472  for (Int_t index = sIndex; index < eIndex; index++) {
473  const Int_t icol = pColIndex[index];
474  if (icol < acoln || icol >= acoln+nr) continue;
475  v[icol-acoln] = pData[index];
476  }
477 }
478 
479 ////////////////////////////////////////////////////////////////////////////////
480 /// General matrix multiplication. Create a matrix C such that C = A * B'.
481 /// Note, matrix C is allocated for constr=1.
482 
483 template<class Element>
484 void TMatrixTSparse<Element>::AMultBt(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
485 {
486  if (gMatrixCheck) {
487  R__ASSERT(a.IsValid());
488  R__ASSERT(b.IsValid());
489 
490  if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
491  Error("AMultBt","A and B columns incompatible");
492  return;
493  }
494 
495  if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
496  Error("AMultB","this = &a");
497  return;
498  }
499 
500  if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
501  Error("AMultB","this = &b");
502  return;
503  }
504  }
505 
506  const Int_t * const pRowIndexa = a.GetRowIndexArray();
507  const Int_t * const pColIndexa = a.GetColIndexArray();
508  const Int_t * const pRowIndexb = b.GetRowIndexArray();
509  const Int_t * const pColIndexb = b.GetColIndexArray();
510 
511  Int_t *pRowIndexc;
512  Int_t *pColIndexc;
513  if (constr) {
514  // make a best guess of the sparse structure; it will guarantee
515  // enough allocated space !
516 
517  Int_t nr_nonzero_rowa = 0;
518  {
519  for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++)
520  if (pRowIndexa[irowa] < pRowIndexa[irowa+1])
521  nr_nonzero_rowa++;
522  }
523  Int_t nr_nonzero_rowb = 0;
524  {
525  for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++)
526  if (pRowIndexb[irowb] < pRowIndexb[irowb+1])
527  nr_nonzero_rowb++;
528  }
529 
530  const Int_t nc = nr_nonzero_rowa*nr_nonzero_rowb; // best guess
531  Allocate(a.GetNrows(),b.GetNrows(),a.GetRowLwb(),b.GetRowLwb(),1,nc);
532 
533  pRowIndexc = this->GetRowIndexArray();
534  pColIndexc = this->GetColIndexArray();
535 
536  pRowIndexc[0] = 0;
537  Int_t ielem = 0;
538  for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++) {
539  pRowIndexc[irowa+1] = pRowIndexc[irowa];
540  if (pRowIndexa[irowa] >= pRowIndexa[irowa+1]) continue;
541  for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++) {
542  if (pRowIndexb[irowb] >= pRowIndexb[irowb+1]) continue;
543  pRowIndexc[irowa+1]++;
544  pColIndexc[ielem++] = irowb;
545  }
546  }
547  } else {
548  pRowIndexc = this->GetRowIndexArray();
549  pColIndexc = this->GetColIndexArray();
550  }
551 
552  const Element * const pDataa = a.GetMatrixArray();
553  const Element * const pDatab = b.GetMatrixArray();
554  Element * const pDatac = this->GetMatrixArray();
555  Int_t indexc_r = 0;
556  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
557  const Int_t sIndexa = pRowIndexa[irowc];
558  const Int_t eIndexa = pRowIndexa[irowc+1];
559  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
560  const Int_t sIndexb = pRowIndexb[icolc];
561  const Int_t eIndexb = pRowIndexb[icolc+1];
562  Element sum = 0.0;
563  Int_t indexb = sIndexb;
564  for (Int_t indexa = sIndexa; indexa < eIndexa && indexb < eIndexb; indexa++) {
565  const Int_t icola = pColIndexa[indexa];
566  while (indexb < eIndexb && pColIndexb[indexb] <= icola) {
567  if (icola == pColIndexb[indexb]) {
568  sum += pDataa[indexa]*pDatab[indexb];
569  break;
570  }
571  indexb++;
572  }
573  }
574  if (sum != 0.0) {
575  pColIndexc[indexc_r] = icolc;
576  pDatac[indexc_r] = sum;
577  indexc_r++;
578  }
579  }
580  pRowIndexc[irowc+1] = indexc_r;
581  }
582 
583  if (constr)
584  SetSparseIndex(indexc_r);
585 }
586 
587 ////////////////////////////////////////////////////////////////////////////////
588 /// General matrix multiplication. Create a matrix C such that C = A * B'.
589 /// Note, matrix C is allocated for constr=1.
590 
591 template<class Element>
592 void TMatrixTSparse<Element>::AMultBt(const TMatrixTSparse<Element> &a,const TMatrixT<Element> &b,Int_t constr)
593 {
594  if (gMatrixCheck) {
595  R__ASSERT(a.IsValid());
596  R__ASSERT(b.IsValid());
597 
598  if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
599  Error("AMultBt","A and B columns incompatible");
600  return;
601  }
602 
603  if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
604  Error("AMultB","this = &a");
605  return;
606  }
607 
608  if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
609  Error("AMultB","this = &b");
610  return;
611  }
612  }
613 
614  const Int_t * const pRowIndexa = a.GetRowIndexArray();
615  const Int_t * const pColIndexa = a.GetColIndexArray();
616 
617  Int_t *pRowIndexc;
618  Int_t *pColIndexc;
619  if (constr) {
620  // make a best guess of the sparse structure; it will guarantee
621  // enough allocated space !
622 
623  Int_t nr_nonzero_rowa = 0;
624  {
625  for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++)
626  if (pRowIndexa[irowa] < pRowIndexa[irowa+1])
627  nr_nonzero_rowa++;
628  }
629  Int_t nr_nonzero_rowb = b.GetNrows();
630 
631  const Int_t nc = nr_nonzero_rowa*nr_nonzero_rowb; // best guess
632  Allocate(a.GetNrows(),b.GetNrows(),a.GetRowLwb(),b.GetRowLwb(),1,nc);
633 
634  pRowIndexc = this->GetRowIndexArray();
635  pColIndexc = this->GetColIndexArray();
636 
637  pRowIndexc[0] = 0;
638  Int_t ielem = 0;
639  for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++) {
640  pRowIndexc[irowa+1] = pRowIndexc[irowa];
641  for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++) {
642  pRowIndexc[irowa+1]++;
643  pColIndexc[ielem++] = irowb;
644  }
645  }
646  } else {
647  pRowIndexc = this->GetRowIndexArray();
648  pColIndexc = this->GetColIndexArray();
649  }
650 
651  const Element * const pDataa = a.GetMatrixArray();
652  const Element * const pDatab = b.GetMatrixArray();
653  Element * const pDatac = this->GetMatrixArray();
654  Int_t indexc_r = 0;
655  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
656  const Int_t sIndexa = pRowIndexa[irowc];
657  const Int_t eIndexa = pRowIndexa[irowc+1];
658  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
659  const Int_t off = icolc*b.GetNcols();
660  Element sum = 0.0;
661  for (Int_t indexa = sIndexa; indexa < eIndexa; indexa++) {
662  const Int_t icola = pColIndexa[indexa];
663  sum += pDataa[indexa]*pDatab[off+icola];
664  }
665  if (sum != 0.0) {
666  pColIndexc[indexc_r] = icolc;
667  pDatac[indexc_r] = sum;
668  indexc_r++;
669  }
670  }
671  pRowIndexc[irowc+1]= indexc_r;
672  }
673 
674  if (constr)
675  SetSparseIndex(indexc_r);
676 }
677 
678 ////////////////////////////////////////////////////////////////////////////////
679 /// General matrix multiplication. Create a matrix C such that C = A * B'.
680 /// Note, matrix C is allocated for constr=1.
681 
682 template<class Element>
683 void TMatrixTSparse<Element>::AMultBt(const TMatrixT<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
684 {
685  if (gMatrixCheck) {
686  R__ASSERT(a.IsValid());
687  R__ASSERT(b.IsValid());
688 
689  if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
690  Error("AMultBt","A and B columns incompatible");
691  return;
692  }
693 
694  if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
695  Error("AMultB","this = &a");
696  return;
697  }
698 
699  if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
700  Error("AMultB","this = &b");
701  return;
702  }
703  }
704 
705  const Int_t * const pRowIndexb = b.GetRowIndexArray();
706  const Int_t * const pColIndexb = b.GetColIndexArray();
707 
708  Int_t *pRowIndexc;
709  Int_t *pColIndexc;
710  if (constr) {
711  // make a best guess of the sparse structure; it will guarantee
712  // enough allocated space !
713 
714  Int_t nr_nonzero_rowa = a.GetNrows();
715  Int_t nr_nonzero_rowb = 0;
716  {
717  for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++)
718  if (pRowIndexb[irowb] < pRowIndexb[irowb+1])
719  nr_nonzero_rowb++;
720  }
721 
722  const Int_t nc = nr_nonzero_rowa*nr_nonzero_rowb; // best guess
723  Allocate(a.GetNrows(),b.GetNrows(),a.GetRowLwb(),b.GetRowLwb(),1,nc);
724 
725  pRowIndexc = this->GetRowIndexArray();
726  pColIndexc = this->GetColIndexArray();
727 
728  pRowIndexc[0] = 0;
729  Int_t ielem = 0;
730  for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++) {
731  pRowIndexc[irowa+1] = pRowIndexc[irowa];
732  for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++) {
733  if (pRowIndexb[irowb] >= pRowIndexb[irowb+1]) continue;
734  pRowIndexc[irowa+1]++;
735  pColIndexc[ielem++] = irowb;
736  }
737  }
738  } else {
739  pRowIndexc = this->GetRowIndexArray();
740  pColIndexc = this->GetColIndexArray();
741  }
742 
743  const Element * const pDataa = a.GetMatrixArray();
744  const Element * const pDatab = b.GetMatrixArray();
745  Element * const pDatac = this->GetMatrixArray();
746  Int_t indexc_r = 0;
747  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
748  const Int_t off = irowc*a.GetNcols();
749  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
750  const Int_t sIndexb = pRowIndexb[icolc];
751  const Int_t eIndexb = pRowIndexb[icolc+1];
752  Element sum = 0.0;
753  for (Int_t indexb = sIndexb; indexb < eIndexb; indexb++) {
754  const Int_t icolb = pColIndexb[indexb];
755  sum += pDataa[off+icolb]*pDatab[indexb];
756  }
757  if (sum != 0.0) {
758  pColIndexc[indexc_r] = icolc;
759  pDatac[indexc_r] = sum;
760  indexc_r++;
761  }
762  }
763  pRowIndexc[irowc+1] = indexc_r;
764  }
765 
766  if (constr)
767  SetSparseIndex(indexc_r);
768 }
769 
770 ////////////////////////////////////////////////////////////////////////////////
771 /// General matrix addition. Create a matrix C such that C = A + B.
772 /// Note, matrix C is allocated for constr=1.
773 
774 template<class Element>
775 void TMatrixTSparse<Element>::APlusB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
776 {
777  if (gMatrixCheck) {
778  R__ASSERT(a.IsValid());
779  R__ASSERT(b.IsValid());
780 
781  if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
782  a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
783  Error("APlusB(const TMatrixTSparse &,const TMatrixTSparse &","matrices not compatible");
784  return;
785  }
786 
787  if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
788  Error("APlusB","this = &a");
789  return;
790  }
791 
792  if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
793  Error("APlusB","this = &b");
794  return;
795  }
796  }
797 
798  const Int_t * const pRowIndexa = a.GetRowIndexArray();
799  const Int_t * const pRowIndexb = b.GetRowIndexArray();
800  const Int_t * const pColIndexa = a.GetColIndexArray();
801  const Int_t * const pColIndexb = b.GetColIndexArray();
802 
803  if (constr) {
804  Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
805  SetSparseIndexAB(a,b);
806  }
807 
808  Int_t * const pRowIndexc = this->GetRowIndexArray();
809  Int_t * const pColIndexc = this->GetColIndexArray();
810 
811  const Element * const pDataa = a.GetMatrixArray();
812  const Element * const pDatab = b.GetMatrixArray();
813  Element * const pDatac = this->GetMatrixArray();
814  Int_t indexc_r = 0;
815  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
816  const Int_t sIndexa = pRowIndexa[irowc];
817  const Int_t eIndexa = pRowIndexa[irowc+1];
818  const Int_t sIndexb = pRowIndexb[irowc];
819  const Int_t eIndexb = pRowIndexb[irowc+1];
820  Int_t indexa = sIndexa;
821  Int_t indexb = sIndexb;
822  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
823  Element sum = 0.0;
824  while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
825  if (icolc == pColIndexa[indexa]) {
826  sum += pDataa[indexa];
827  break;
828  }
829  indexa++;
830  }
831  while (indexb < eIndexb && pColIndexb[indexb] <= icolc) {
832  if (icolc == pColIndexb[indexb]) {
833  sum += pDatab[indexb];
834  break;
835  }
836  indexb++;
837  }
838 
839  if (sum != 0.0) {
840  pColIndexc[indexc_r] = icolc;
841  pDatac[indexc_r] = sum;
842  indexc_r++;
843  }
844  }
845  pRowIndexc[irowc+1] = indexc_r;
846  }
847 
848  if (constr)
849  SetSparseIndex(indexc_r);
850 }
851 
852 ////////////////////////////////////////////////////////////////////////////////
853 /// General matrix addition. Create a matrix C such that C = A + B.
854 /// Note, matrix C is allocated for constr=1.
855 
856 template<class Element>
857 void TMatrixTSparse<Element>::APlusB(const TMatrixTSparse<Element> &a,const TMatrixT<Element> &b,Int_t constr)
858 {
859  if (gMatrixCheck) {
860  R__ASSERT(a.IsValid());
861  R__ASSERT(b.IsValid());
862 
863  if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
864  a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
865  Error("APlusB(const TMatrixTSparse &,const TMatrixT &","matrices not compatible");
866  return;
867  }
868 
869  if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
870  Error("APlusB","this = &a");
871  return;
872  }
873 
874  if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
875  Error("APlusB","this = &b");
876  return;
877  }
878  }
879 
880  if (constr) {
881  Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
882  SetSparseIndexAB(a,b);
883  }
884 
885  Int_t * const pRowIndexc = this->GetRowIndexArray();
886  Int_t * const pColIndexc = this->GetColIndexArray();
887 
888  const Int_t * const pRowIndexa = a.GetRowIndexArray();
889  const Int_t * const pColIndexa = a.GetColIndexArray();
890 
891  const Element * const pDataa = a.GetMatrixArray();
892  const Element * const pDatab = b.GetMatrixArray();
893  Element * const pDatac = this->GetMatrixArray();
894  Int_t indexc_r = 0;
895  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
896  const Int_t sIndexa = pRowIndexa[irowc];
897  const Int_t eIndexa = pRowIndexa[irowc+1];
898  const Int_t off = irowc*this->GetNcols();
899  Int_t indexa = sIndexa;
900  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
901  Element sum = pDatab[off+icolc];
902  while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
903  if (icolc == pColIndexa[indexa]) {
904  sum += pDataa[indexa];
905  break;
906  }
907  indexa++;
908  }
909 
910  if (sum != 0.0) {
911  pColIndexc[indexc_r] = icolc;
912  pDatac[indexc_r] = sum;
913  indexc_r++;
914  }
915  }
916  pRowIndexc[irowc+1] = indexc_r;
917  }
918 
919  if (constr)
920  SetSparseIndex(indexc_r);
921 }
922 
923 ////////////////////////////////////////////////////////////////////////////////
924 /// General matrix subtraction. Create a matrix C such that C = A - B.
925 /// Note, matrix C is allocated for constr=1.
926 
927 template<class Element>
928 void TMatrixTSparse<Element>::AMinusB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
929 {
930  if (gMatrixCheck) {
931  R__ASSERT(a.IsValid());
932  R__ASSERT(b.IsValid());
933 
934  if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
935  a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
936  Error("AMinusB(const TMatrixTSparse &,const TMatrixTSparse &","matrices not compatible");
937  return;
938  }
939 
940  if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
941  Error("AMinusB","this = &a");
942  return;
943  }
944 
945  if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
946  Error("AMinusB","this = &b");
947  return;
948  }
949  }
950 
951  const Int_t * const pRowIndexa = a.GetRowIndexArray();
952  const Int_t * const pRowIndexb = b.GetRowIndexArray();
953  const Int_t * const pColIndexa = a.GetColIndexArray();
954  const Int_t * const pColIndexb = b.GetColIndexArray();
955 
956  if (constr) {
957  Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
958  SetSparseIndexAB(a,b);
959  }
960 
961  Int_t * const pRowIndexc = this->GetRowIndexArray();
962  Int_t * const pColIndexc = this->GetColIndexArray();
963 
964  const Element * const pDataa = a.GetMatrixArray();
965  const Element * const pDatab = b.GetMatrixArray();
966  Element * const pDatac = this->GetMatrixArray();
967  Int_t indexc_r = 0;
968  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
969  const Int_t sIndexa = pRowIndexa[irowc];
970  const Int_t eIndexa = pRowIndexa[irowc+1];
971  const Int_t sIndexb = pRowIndexb[irowc];
972  const Int_t eIndexb = pRowIndexb[irowc+1];
973  Int_t indexa = sIndexa;
974  Int_t indexb = sIndexb;
975  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
976  Element sum = 0.0;
977  while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
978  if (icolc == pColIndexa[indexa]) {
979  sum += pDataa[indexa];
980  break;
981  }
982  indexa++;
983  }
984  while (indexb < eIndexb && pColIndexb[indexb] <= icolc) {
985  if (icolc == pColIndexb[indexb]) {
986  sum -= pDatab[indexb];
987  break;
988  }
989  indexb++;
990  }
991 
992  if (sum != 0.0) {
993  pColIndexc[indexc_r] = icolc;
994  pDatac[indexc_r] = sum;
995  indexc_r++;
996  }
997  }
998  pRowIndexc[irowc+1] = indexc_r;
999  }
1000 
1001  if (constr)
1002  SetSparseIndex(indexc_r);
1003 }
1004 
1005 ////////////////////////////////////////////////////////////////////////////////
1006 /// General matrix subtraction. Create a matrix C such that C = A - B.
1007 /// Note, matrix C is allocated for constr=1.
1008 
1009 template<class Element>
1010 void TMatrixTSparse<Element>::AMinusB(const TMatrixTSparse<Element> &a,const TMatrixT<Element> &b,Int_t constr)
1011 {
1012  if (gMatrixCheck) {
1013  R__ASSERT(a.IsValid());
1014  R__ASSERT(b.IsValid());
1015 
1016  if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
1017  a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
1018  Error("AMinusB(const TMatrixTSparse &,const TMatrixT &","matrices not compatible");
1019  return;
1020  }
1021 
1022  if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
1023  Error("AMinusB","this = &a");
1024  return;
1025  }
1026 
1027  if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
1028  Error("AMinusB","this = &b");
1029  return;
1030  }
1031  }
1032 
1033  if (constr) {
1034  Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
1035  SetSparseIndexAB(a,b);
1036  }
1037 
1038  Int_t * const pRowIndexc = this->GetRowIndexArray();
1039  Int_t * const pColIndexc = this->GetColIndexArray();
1040 
1041  const Int_t * const pRowIndexa = a.GetRowIndexArray();
1042  const Int_t * const pColIndexa = a.GetColIndexArray();
1043 
1044  const Element * const pDataa = a.GetMatrixArray();
1045  const Element * const pDatab = b.GetMatrixArray();
1046  Element * const pDatac = this->GetMatrixArray();
1047  Int_t indexc_r = 0;
1048  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
1049  const Int_t sIndexa = pRowIndexa[irowc];
1050  const Int_t eIndexa = pRowIndexa[irowc+1];
1051  const Int_t off = irowc*this->GetNcols();
1052  Int_t indexa = sIndexa;
1053  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
1054  Element sum = -pDatab[off+icolc];
1055  while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
1056  if (icolc == pColIndexa[indexa]) {
1057  sum += pDataa[indexa];
1058  break;
1059  }
1060  indexa++;
1061  }
1062 
1063  if (sum != 0.0) {
1064  pColIndexc[indexc_r] = icolc;
1065  pDatac[indexc_r] = sum;
1066  indexc_r++;
1067  }
1068  }
1069  pRowIndexc[irowc+1] = indexc_r;
1070  }
1071 
1072  if (constr)
1073  SetSparseIndex(indexc_r);
1074 }
1075 
1076 ////////////////////////////////////////////////////////////////////////////////
1077 /// General matrix subtraction. Create a matrix C such that C = A - B.
1078 /// Note, matrix C is allocated for constr=1.
1079 
1080 template<class Element>
1081 void TMatrixTSparse<Element>::AMinusB(const TMatrixT<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
1082 {
1083  if (gMatrixCheck) {
1084  R__ASSERT(a.IsValid());
1085  R__ASSERT(b.IsValid());
1086 
1087  if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
1088  a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
1089  Error("AMinusB(const TMatrixT &,const TMatrixTSparse &","matrices not compatible");
1090  return;
1091  }
1092 
1093  if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
1094  Error("AMinusB","this = &a");
1095  return;
1096  }
1097 
1098  if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
1099  Error("AMinusB","this = &b");
1100  return;
1101  }
1102  }
1103 
1104  if (constr) {
1105  Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
1106  SetSparseIndexAB(a,b);
1107  }
1108 
1109  Int_t * const pRowIndexc = this->GetRowIndexArray();
1110  Int_t * const pColIndexc = this->GetColIndexArray();
1111 
1112  const Int_t * const pRowIndexb = b.GetRowIndexArray();
1113  const Int_t * const pColIndexb = b.GetColIndexArray();
1114 
1115  const Element * const pDataa = a.GetMatrixArray();
1116  const Element * const pDatab = b.GetMatrixArray();
1117  Element * const pDatac = this->GetMatrixArray();
1118  Int_t indexc_r = 0;
1119  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
1120  const Int_t sIndexb = pRowIndexb[irowc];
1121  const Int_t eIndexb = pRowIndexb[irowc+1];
1122  const Int_t off = irowc*this->GetNcols();
1123  Int_t indexb = sIndexb;
1124  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
1125  Element sum = pDataa[off+icolc];
1126  while (indexb < eIndexb && pColIndexb[indexb] <= icolc) {
1127  if (icolc == pColIndexb[indexb]) {
1128  sum -= pDatab[indexb];
1129  break;
1130  }
1131  indexb++;
1132  }
1133 
1134  if (sum != 0.0) {
1135  pColIndexc[indexc_r] = icolc;
1136  pDatac[indexc_r] = sum;
1137  indexc_r++;
1138  }
1139  }
1140  pRowIndexc[irowc+1] = indexc_r;
1141  }
1142 
1143  if (constr)
1144  SetSparseIndex(indexc_r);
1145 }
1146 
1147 ////////////////////////////////////////////////////////////////////////////////
1148 /// Copy matrix data to array . It is assumed that array is of size >= fNelems
1149 
1150 template<class Element>
1151 void TMatrixTSparse<Element>::GetMatrix2Array(Element *data,Option_t * /*option*/) const
1152 {
1153  R__ASSERT(this->IsValid());
1154 
1155  const Element * const elem = GetMatrixArray();
1156  memcpy(data,elem,this->fNelems*sizeof(Element));
1157 }
1158 
1159 ////////////////////////////////////////////////////////////////////////////////
1160 /// Copy nr elements from row/col index and data array to matrix . It is assumed
1161 /// that arrays are of size >= nr
1162 /// Note that the input arrays are not passed as const since they will be modified !
1163 
1164 template<class Element>
1165 TMatrixTBase<Element> &TMatrixTSparse<Element>::SetMatrixArray(Int_t nr,Int_t *row,Int_t *col,Element *data)
1166 {
1167  R__ASSERT(this->IsValid());
1168  if (nr <= 0) {
1169  Error("SetMatrixArray(Int_t,Int_t*,Int_t*,Element*","nr <= 0");
1170  return *this;
1171  }
1172 
1173  const Int_t irowmin = TMath::LocMin(nr,row);
1174  const Int_t irowmax = TMath::LocMax(nr,row);
1175  const Int_t icolmin = TMath::LocMin(nr,col);
1176  const Int_t icolmax = TMath::LocMax(nr,col);
1177 
1178  R__ASSERT(row[irowmin] >= this->fRowLwb && row[irowmax] <= this->fRowLwb+this->fNrows-1);
1179  R__ASSERT(col[icolmin] >= this->fColLwb && col[icolmax] <= this->fColLwb+this->fNcols-1);
1180 
1181  if (row[irowmin] < this->fRowLwb || row[irowmax] > this->fRowLwb+this->fNrows-1) {
1182  Error("SetMatrixArray","Inconsistency between row index and its range");
1183  if (row[irowmin] < this->fRowLwb) {
1184  Info("SetMatrixArray","row index lower bound adjusted to %d",row[irowmin]);
1185  this->fRowLwb = row[irowmin];
1186  }
1187  if (row[irowmax] > this->fRowLwb+this->fNrows-1) {
1188  Info("SetMatrixArray","row index upper bound adjusted to %d",row[irowmax]);
1189  this->fNrows = row[irowmax]-this->fRowLwb+1;
1190  }
1191  }
1192  if (col[icolmin] < this->fColLwb || col[icolmax] > this->fColLwb+this->fNcols-1) {
1193  Error("SetMatrixArray","Inconsistency between column index and its range");
1194  if (col[icolmin] < this->fColLwb) {
1195  Info("SetMatrixArray","column index lower bound adjusted to %d",col[icolmin]);
1196  this->fColLwb = col[icolmin];
1197  }
1198  if (col[icolmax] > this->fColLwb+this->fNcols-1) {
1199  Info("SetMatrixArray","column index upper bound adjusted to %d",col[icolmax]);
1200  this->fNcols = col[icolmax]-this->fColLwb+1;
1201  }
1202  }
1203 
1204  TMatrixTBase<Element>::DoubleLexSort(nr,row,col,data);
1205 
1206  Int_t nr_nonzeros = 0;
1207  const Element *ep = data;
1208  const Element * const fp = data+nr;
1209 
1210  while (ep < fp)
1211  if (*ep++ != 0.0) nr_nonzeros++;
1212 
1213  if (nr_nonzeros != this->fNelems) {
1214  if (fColIndex) { delete [] fColIndex; fColIndex = 0; }
1215  if (fElements) { delete [] fElements; fElements = 0; }
1216  this->fNelems = nr_nonzeros;
1217  if (this->fNelems > 0) {
1218  fColIndex = new Int_t[nr_nonzeros];
1219  fElements = new Element[nr_nonzeros];
1220  } else {
1221  fColIndex = 0;
1222  fElements = 0;
1223  }
1224  }
1225 
1226  if (this->fNelems <= 0)
1227  return *this;
1228 
1229  fRowIndex[0] = 0;
1230  Int_t ielem = 0;
1231  nr_nonzeros = 0;
1232  for (Int_t irow = 1; irow < this->fNrows+1; irow++) {
1233  if (ielem < nr && row[ielem] < irow) {
1234  while (ielem < nr) {
1235  if (data[ielem] != 0.0) {
1236  fColIndex[nr_nonzeros] = col[ielem]-this->fColLwb;
1237  fElements[nr_nonzeros] = data[ielem];
1238  nr_nonzeros++;
1239  }
1240  ielem++;
1241  if (ielem >= nr || row[ielem] != row[ielem-1])
1242  break;
1243  }
1244  }
1245  fRowIndex[irow] = nr_nonzeros;
1246  }
1247 
1248  return *this;
1249 }
1250 
1251 ////////////////////////////////////////////////////////////////////////////////
1252 /// Increase/decrease the number of non-zero elements to nelems_new
1253 
1254 template<class Element>
1255 TMatrixTSparse<Element> &TMatrixTSparse<Element>::SetSparseIndex(Int_t nelems_new)
1256 {
1257  if (nelems_new != this->fNelems) {
1258  Int_t nr = TMath::Min(nelems_new,this->fNelems);
1259  Int_t *oIp = fColIndex;
1260  fColIndex = new Int_t[nelems_new];
1261  if (oIp) {
1262  memmove(fColIndex,oIp,nr*sizeof(Int_t));
1263  delete [] oIp;
1264  }
1265  Element *oDp = fElements;
1266  fElements = new Element[nelems_new];
1267  if (oDp) {
1268  memmove(fElements,oDp,nr*sizeof(Element));
1269  delete [] oDp;
1270  }
1271  this->fNelems = nelems_new;
1272  if (nelems_new > nr) {
1273  memset(fElements+nr,0,(nelems_new-nr)*sizeof(Element));
1274  memset(fColIndex+nr,0,(nelems_new-nr)*sizeof(Int_t));
1275  } else {
1276  for (Int_t irow = 0; irow < this->fNrowIndex; irow++)
1277  if (fRowIndex[irow] > nelems_new)
1278  fRowIndex[irow] = nelems_new;
1279  }
1280  }
1281 
1282  return *this;
1283 }
1284 
1285 ////////////////////////////////////////////////////////////////////////////////
1286 /// Use non-zero data of matrix source to set the sparse structure
1287 
1288 template<class Element>
1289 TMatrixTSparse<Element> &TMatrixTSparse<Element>::SetSparseIndex(const TMatrixTBase<Element> &source)
1290 {
1291  if (gMatrixCheck) {
1292  R__ASSERT(source.IsValid());
1293  if (this->GetNrows() != source.GetNrows() || this->GetNcols() != source.GetNcols() ||
1294  this->GetRowLwb() != source.GetRowLwb() || this->GetColLwb() != source.GetColLwb()) {
1295  Error("SetSparseIndex","matrices not compatible");
1296  return *this;
1297  }
1298  }
1299 
1300  const Int_t nr_nonzeros = source.NonZeros();
1301 
1302  if (nr_nonzeros != this->fNelems)
1303  SetSparseIndex(nr_nonzeros);
1304 
1305  if (source.GetRowIndexArray() && source.GetColIndexArray()) {
1306  memmove(fRowIndex,source.GetRowIndexArray(),this->fNrowIndex*sizeof(Int_t));
1307  memmove(fColIndex,source.GetColIndexArray(),this->fNelems*sizeof(Int_t));
1308  } else {
1309  const Element *ep = source.GetMatrixArray();
1310  Int_t nr = 0;
1311  for (Int_t irow = 0; irow < this->fNrows; irow++) {
1312  fRowIndex[irow] = nr;
1313  for (Int_t icol = 0; icol < this->fNcols; icol++) {
1314  if (*ep != 0.0) {
1315  fColIndex[nr] = icol;
1316  nr++;
1317  }
1318  ep++;
1319  }
1320  }
1321  fRowIndex[this->fNrows] = nr;
1322  }
1323 
1324  return *this;
1325 }
1326 
1327 ////////////////////////////////////////////////////////////////////////////////
1328 /// Set the row/column indices to the "sum" of matrices a and b
1329 /// It is checked that enough space has been allocated
1330 
1331 template<class Element>
1332 TMatrixTSparse<Element> &TMatrixTSparse<Element>::SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b)
1333 {
1334  if (gMatrixCheck) {
1335  R__ASSERT(a.IsValid());
1336  R__ASSERT(b.IsValid());
1337 
1338  if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
1339  a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
1340  Error("SetSparseIndexAB","source matrices not compatible");
1341  return *this;
1342  }
1343 
1344  if (this->GetNrows() != a.GetNrows() || this->GetNcols() != a.GetNcols() ||
1345  this->GetRowLwb() != a.GetRowLwb() || this->GetColLwb() != a.GetColLwb()) {
1346  Error("SetSparseIndexAB","matrix not compatible with source matrices");
1347  return *this;
1348  }
1349  }
1350 
1351  const Int_t * const pRowIndexa = a.GetRowIndexArray();
1352  const Int_t * const pRowIndexb = b.GetRowIndexArray();
1353  const Int_t * const pColIndexa = a.GetColIndexArray();
1354  const Int_t * const pColIndexb = b.GetColIndexArray();
1355 
1356  // Count first the number of non-zero slots that are needed
1357  Int_t nc = 0;
1358  for (Int_t irowc = 0; irowc < a.GetNrows(); irowc++) {
1359  const Int_t sIndexa = pRowIndexa[irowc];
1360  const Int_t eIndexa = pRowIndexa[irowc+1];
1361  const Int_t sIndexb = pRowIndexb[irowc];
1362  const Int_t eIndexb = pRowIndexb[irowc+1];
1363  nc += eIndexa-sIndexa;
1364  Int_t indexb = sIndexb;
1365  for (Int_t indexa = sIndexa; indexa < eIndexa; indexa++) {
1366  const Int_t icola = pColIndexa[indexa];
1367  for (; indexb < eIndexb; indexb++) {
1368  if (pColIndexb[indexb] >= icola) {
1369  if (pColIndexb[indexb] == icola)
1370  indexb++;
1371  break;
1372  }
1373  nc++;
1374  }
1375  }
1376  while (indexb < eIndexb) {
1377  const Int_t icola = (eIndexa > sIndexa && eIndexa > 0) ? pColIndexa[eIndexa-1] : -1;
1378  if (pColIndexb[indexb++] > icola)
1379  nc++;
1380  }
1381  }
1382 
1383  // Allocate the necessary space in fRowIndex and fColIndex
1384  if (this->NonZeros() != nc)
1385  SetSparseIndex(nc);
1386 
1387  Int_t * const pRowIndexc = this->GetRowIndexArray();
1388  Int_t * const pColIndexc = this->GetColIndexArray();
1389 
1390  nc = 0;
1391  pRowIndexc[0] = 0;
1392  for (Int_t irowc = 0; irowc < a.GetNrows(); irowc++) {
1393  const Int_t sIndexa = pRowIndexa[irowc];
1394  const Int_t eIndexa = pRowIndexa[irowc+1];
1395  const Int_t sIndexb = pRowIndexb[irowc];
1396  const Int_t eIndexb = pRowIndexb[irowc+1];
1397  Int_t indexb = sIndexb;
1398  for (Int_t indexa = sIndexa; indexa < eIndexa; indexa++) {
1399  const Int_t icola = pColIndexa[indexa];
1400  for (; indexb < eIndexb; indexb++) {
1401  if (pColIndexb[indexb] >= icola) {
1402  if (pColIndexb[indexb] == icola)
1403  indexb++;
1404  break;
1405  }
1406  pColIndexc[nc++] = pColIndexb[indexb];
1407  }
1408  pColIndexc[nc++] = pColIndexa[indexa];
1409  }
1410  while (indexb < eIndexb) {
1411  const Int_t icola = (eIndexa > 0) ? pColIndexa[eIndexa-1] : -1;
1412  if (pColIndexb[indexb++] > icola)
1413  pColIndexc[nc++] = pColIndexb[indexb-1];
1414  }
1415  pRowIndexc[irowc+1] = nc;
1416  }
1417 
1418  return *this;
1419 }
1420 
1421 ////////////////////////////////////////////////////////////////////////////////
1422 /// Set the row/column indices to the "sum" of matrices a and b
1423 /// It is checked that enough space has been allocated
1424 
1425 template<class Element>
1426 TMatrixTSparse<Element> &TMatrixTSparse<Element>::SetSparseIndexAB(const TMatrixT<Element> &a,const TMatrixTSparse<Element> &b)
1427 {
1428  if (gMatrixCheck) {
1429  R__ASSERT(a.IsValid());
1430  R__ASSERT(b.IsValid());
1431 
1432  if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
1433  a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
1434  Error("SetSparseIndexAB","source matrices not compatible");
1435  return *this;
1436  }
1437 
1438  if (this->GetNrows() != a.GetNrows() || this->GetNcols() != a.GetNcols() ||
1439  this->GetRowLwb() != a.GetRowLwb() || this->GetColLwb() != a.GetColLwb()) {
1440  Error("SetSparseIndexAB","matrix not compatible with source matrices");
1441  return *this;
1442  }
1443  }
1444 
1445  const Element * const pDataa = a.GetMatrixArray();
1446 
1447  const Int_t * const pRowIndexb = b.GetRowIndexArray();
1448  const Int_t * const pColIndexb = b.GetColIndexArray();
1449 
1450  // Count first the number of non-zero slots that are needed
1451  Int_t nc = a.NonZeros();
1452  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
1453  const Int_t sIndexb = pRowIndexb[irowc];
1454  const Int_t eIndexb = pRowIndexb[irowc+1];
1455  const Int_t off = irowc*this->GetNcols();
1456  Int_t indexb = sIndexb;
1457  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
1458  if (pDataa[off+icolc] != 0.0 || pColIndexb[indexb] > icolc) continue;
1459  for (; indexb < eIndexb; indexb++) {
1460  if (pColIndexb[indexb] >= icolc) {
1461  if (pColIndexb[indexb] == icolc) {
1462  nc++;
1463  indexb++;
1464  }
1465  break;
1466  }
1467  }
1468  }
1469  }
1470 
1471  // Allocate thre necessary space in fRowIndex and fColIndex
1472  if (this->NonZeros() != nc)
1473  SetSparseIndex(nc);
1474 
1475  Int_t * const pRowIndexc = this->GetRowIndexArray();
1476  Int_t * const pColIndexc = this->GetColIndexArray();
1477 
1478  nc = 0;
1479  pRowIndexc[0] = 0;
1480  for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
1481  const Int_t sIndexb = pRowIndexb[irowc];
1482  const Int_t eIndexb = pRowIndexb[irowc+1];
1483  const Int_t off = irowc*this->GetNcols();
1484  Int_t indexb = sIndexb;
1485  for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
1486  if (pDataa[off+icolc] != 0.0)
1487  pColIndexc[nc++] = icolc;
1488  else if (pColIndexb[indexb] <= icolc) {
1489  for (; indexb < eIndexb; indexb++) {
1490  if (pColIndexb[indexb] >= icolc) {
1491  if (pColIndexb[indexb] == icolc)
1492  pColIndexc[nc++] = pColIndexb[indexb++];
1493  break;
1494  }
1495  }
1496  }
1497  }
1498  pRowIndexc[irowc+1] = nc;
1499  }
1500 
1501  return *this;
1502 }
1503 
1504 ////////////////////////////////////////////////////////////////////////////////
1505 /// Set size of the matrix to nrows x ncols with nr_nonzeros non-zero entries
1506 /// if nr_nonzeros > 0 .
1507 /// New dynamic elements are created, the overlapping part of the old ones are
1508 /// copied to the new structures, then the old elements are deleted.
1509 
1510 template<class Element>
1511 TMatrixTBase<Element> &TMatrixTSparse<Element>::ResizeTo(Int_t nrows,Int_t ncols,Int_t nr_nonzeros)
1512 {
1513  R__ASSERT(this->IsValid());
1514  if (!this->fIsOwner) {
1515  Error("ResizeTo(Int_t,Int_t,Int_t)","Not owner of data array,cannot resize");
1516  return *this;
1517  }
1518 
1519  if (this->fNelems > 0) {
1520  if (this->fNrows == nrows && this->fNcols == ncols &&
1521  (this->fNelems == nr_nonzeros || nr_nonzeros < 0))
1522  return *this;
1523  else if (nrows == 0 || ncols == 0 || nr_nonzeros == 0) {
1524  this->fNrows = nrows; this->fNcols = ncols;
1525  Clear();
1526  return *this;
1527  }
1528 
1529  const Element *elements_old = GetMatrixArray();
1530  const Int_t *rowIndex_old = GetRowIndexArray();
1531  const Int_t *colIndex_old = GetColIndexArray();
1532 
1533  Int_t nrows_old = this->fNrows;
1534  Int_t nrowIndex_old = this->fNrowIndex;
1535 
1536  Int_t nelems_new;
1537  if (nr_nonzeros > 0)
1538  nelems_new = nr_nonzeros;
1539  else {
1540  nelems_new = 0;
1541  for (Int_t irow = 0; irow < nrows_old; irow++) {
1542  if (irow >= nrows) continue;
1543  const Int_t sIndex = rowIndex_old[irow];
1544  const Int_t eIndex = rowIndex_old[irow+1];
1545  for (Int_t index = sIndex; index < eIndex; index++) {
1546  const Int_t icol = colIndex_old[index];
1547  if (icol < ncols)
1548  nelems_new++;
1549  }
1550  }
1551  }
1552 
1553  Allocate(nrows,ncols,0,0,1,nelems_new);
1554  R__ASSERT(this->IsValid());
1555 
1556  Element *elements_new = GetMatrixArray();
1557  Int_t *rowIndex_new = GetRowIndexArray();
1558  Int_t *colIndex_new = GetColIndexArray();
1559 
1560  Int_t nelems_copy = 0;
1561  rowIndex_new[0] = 0;
1562  Bool_t cont = kTRUE;
1563  for (Int_t irow = 0; irow < nrows_old && cont; irow++) {
1564  if (irow >= nrows) continue;
1565  const Int_t sIndex = rowIndex_old[irow];
1566  const Int_t eIndex = rowIndex_old[irow+1];
1567  for (Int_t index = sIndex; index < eIndex; index++) {
1568  const Int_t icol = colIndex_old[index];
1569  if (icol < ncols) {
1570  rowIndex_new[irow+1] = nelems_copy+1;
1571  colIndex_new[nelems_copy] = icol;
1572  elements_new[nelems_copy] = elements_old[index];
1573  nelems_copy++;
1574  }
1575  if (nelems_copy >= nelems_new) {
1576  cont = kFALSE;
1577  break;
1578  }
1579  }
1580  }
1581 
1582  if (rowIndex_old) delete [] (Int_t*) rowIndex_old;
1583  if (colIndex_old) delete [] (Int_t*) colIndex_old;
1584  if (elements_old) delete [] (Element*) elements_old;
1585 
1586  if (nrowIndex_old < this->fNrowIndex) {
1587  for (Int_t irow = nrowIndex_old; irow < this->fNrowIndex; irow++)
1588  rowIndex_new[irow] = rowIndex_new[nrowIndex_old-1];
1589  }
1590  } else {
1591  const Int_t nelems_new = (nr_nonzeros >= 0) ? nr_nonzeros : 0;
1592  Allocate(nrows,ncols,0,0,1,nelems_new);
1593  }
1594 
1595  return *this;
1596 }
1597 
1598 ////////////////////////////////////////////////////////////////////////////////
1599 /// Set size of the matrix to [row_lwb:row_upb] x [col_lwb:col_upb] with nr_nonzeros
1600 /// non-zero entries if nr_nonzeros > 0 .
1601 /// New dynamic elemenst are created, the overlapping part of the old ones are
1602 /// copied to the new structures, then the old elements are deleted.
1603 
1604 template<class Element>
1605 TMatrixTBase<Element> &TMatrixTSparse<Element>::ResizeTo(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
1606  Int_t nr_nonzeros)
1607 {
1608  R__ASSERT(this->IsValid());
1609  if (!this->fIsOwner) {
1610  Error("ResizeTo(Int_t,Int_t,Int_t,Int_t,Int_t)","Not owner of data array,cannot resize");
1611  return *this;
1612  }
1613 
1614  const Int_t new_nrows = row_upb-row_lwb+1;
1615  const Int_t new_ncols = col_upb-col_lwb+1;
1616 
1617  if (this->fNelems > 0) {
1618  if (this->fNrows == new_nrows && this->fNcols == new_ncols &&
1619  this->fRowLwb == row_lwb && this->fColLwb == col_lwb &&
1620  (this->fNelems == nr_nonzeros || nr_nonzeros < 0))
1621  return *this;
1622  else if (new_nrows == 0 || new_ncols == 0 || nr_nonzeros == 0) {
1623  this->fNrows = new_nrows; this->fNcols = new_ncols;
1624  this->fRowLwb = row_lwb; this->fColLwb = col_lwb;
1625  Clear();
1626  return *this;
1627  }
1628 
1629  const Int_t *rowIndex_old = GetRowIndexArray();
1630  const Int_t *colIndex_old = GetColIndexArray();
1631  const Element *elements_old = GetMatrixArray();
1632 
1633  const Int_t nrowIndex_old = this->fNrowIndex;
1634 
1635  const Int_t nrows_old = this->fNrows;
1636  const Int_t rowLwb_old = this->fRowLwb;
1637  const Int_t colLwb_old = this->fColLwb;
1638 
1639  Int_t nelems_new;
1640  if (nr_nonzeros > 0)
1641  nelems_new = nr_nonzeros;
1642  else {
1643  nelems_new = 0;
1644  for (Int_t irow = 0; irow < nrows_old; irow++) {
1645  if (irow+rowLwb_old > row_upb || irow+rowLwb_old < row_lwb) continue;
1646  const Int_t sIndex = rowIndex_old[irow];
1647  const Int_t eIndex = rowIndex_old[irow+1];
1648  for (Int_t index = sIndex; index < eIndex; index++) {
1649  const Int_t icol = colIndex_old[index];
1650  if (icol+colLwb_old <= col_upb && icol+colLwb_old >= col_lwb)
1651  nelems_new++;
1652  }
1653  }
1654  }
1655 
1656  Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1,nelems_new);
1657  R__ASSERT(this->IsValid());
1658 
1659  Int_t *rowIndex_new = GetRowIndexArray();
1660  Int_t *colIndex_new = GetColIndexArray();
1661  Element *elements_new = GetMatrixArray();
1662 
1663  Int_t nelems_copy = 0;
1664  rowIndex_new[0] = 0;
1665  const Int_t row_off = rowLwb_old-row_lwb;
1666  const Int_t col_off = colLwb_old-col_lwb;
1667  for (Int_t irow = 0; irow < nrows_old; irow++) {
1668  if (irow+rowLwb_old > row_upb || irow+rowLwb_old < row_lwb) continue;
1669  const Int_t sIndex = rowIndex_old[irow];
1670  const Int_t eIndex = rowIndex_old[irow+1];
1671  for (Int_t index = sIndex; index < eIndex; index++) {
1672  const Int_t icol = colIndex_old[index];
1673  if (icol+colLwb_old <= col_upb && icol+colLwb_old >= col_lwb) {
1674  rowIndex_new[irow+row_off+1] = nelems_copy+1;
1675  colIndex_new[nelems_copy] = icol+col_off;
1676  elements_new[nelems_copy] = elements_old[index];
1677  nelems_copy++;
1678  }
1679  if (nelems_copy >= nelems_new) {
1680  break;
1681  }
1682  }
1683  }
1684 
1685  if (rowIndex_old) delete [] (Int_t*) rowIndex_old;
1686  if (colIndex_old) delete [] (Int_t*) colIndex_old;
1687  if (elements_old) delete [] (Element*) elements_old;
1688 
1689  if (nrowIndex_old < this->fNrowIndex) {
1690  for (Int_t irow = nrowIndex_old; irow < this->fNrowIndex; irow++)
1691  rowIndex_new[irow] = rowIndex_new[nrowIndex_old-1];
1692  }
1693  } else {
1694  const Int_t nelems_new = (nr_nonzeros >= 0) ? nr_nonzeros : 0;
1695  Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1,nelems_new);
1696  }
1697 
1698  return *this;
1699 }
1700 
1701 ////////////////////////////////////////////////////////////////////////////////
1702 
1703 template<class Element>
1704 TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
1705  Int_t nr_nonzeros,Int_t *pRowIndex,Int_t *pColIndex,Element *pData)
1706 {
1707  if (gMatrixCheck) {
1708  if (row_upb < row_lwb) {
1709  Error("Use","row_upb=%d < row_lwb=%d",row_upb,row_lwb);
1710  return *this;
1711  }
1712  if (col_upb < col_lwb) {
1713  Error("Use","col_upb=%d < col_lwb=%d",col_upb,col_lwb);
1714  return *this;
1715  }
1716  }
1717 
1718  Clear();
1719 
1720  this->fNrows = row_upb-row_lwb+1;
1721  this->fNcols = col_upb-col_lwb+1;
1722  this->fRowLwb = row_lwb;
1723  this->fColLwb = col_lwb;
1724  this->fNrowIndex = this->fNrows+1;
1725  this->fNelems = nr_nonzeros;
1726  this->fIsOwner = kFALSE;
1727  this->fTol = std::numeric_limits<Element>::epsilon();
1728 
1729  fElements = pData;
1730  fRowIndex = pRowIndex;
1731  fColIndex = pColIndex;
1732 
1733  return *this;
1734 }
1735 
1736 ////////////////////////////////////////////////////////////////////////////////
1737 /// Get submatrix [row_lwb..row_upb][col_lwb..col_upb]; The indexing range of the
1738 /// returned matrix depends on the argument option:
1739 ///
1740 /// option == "S" : return [0..row_upb-row_lwb+1][0..col_upb-col_lwb+1] (default)
1741 /// else : return [row_lwb..row_upb][col_lwb..col_upb]
1742 
1743 template<class Element>
1744 TMatrixTBase<Element> &TMatrixTSparse<Element>::GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
1745  TMatrixTBase<Element> &target,Option_t *option) const
1746 {
1747  if (gMatrixCheck) {
1748  R__ASSERT(this->IsValid());
1749  if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
1750  Error("GetSub","row_lwb out-of-bounds");
1751  return target;
1752  }
1753  if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
1754  Error("GetSub","col_lwb out-of-bounds");
1755  return target;
1756  }
1757  if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
1758  Error("GetSub","row_upb out-of-bounds");
1759  return target;
1760  }
1761  if (col_upb < this->fColLwb || col_upb > this->fColLwb+this->fNcols-1) {
1762  Error("GetSub","col_upb out-of-bounds");
1763  return target;
1764  }
1765  if (row_upb < row_lwb || col_upb < col_lwb) {
1766  Error("GetSub","row_upb < row_lwb || col_upb < col_lwb");
1767  return target;
1768  }
1769  }
1770 
1771  TString opt(option);
1772  opt.ToUpper();
1773  const Int_t shift = (opt.Contains("S")) ? 1 : 0;
1774 
1775  const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
1776  const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
1777  const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
1778  const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
1779 
1780  Int_t nr_nonzeros = 0;
1781  for (Int_t irow = 0; irow < this->fNrows; irow++) {
1782  if (irow+this->fRowLwb > row_upb || irow+this->fRowLwb < row_lwb) continue;
1783  const Int_t sIndex = fRowIndex[irow];
1784  const Int_t eIndex = fRowIndex[irow+1];
1785  for (Int_t index = sIndex; index < eIndex; index++) {
1786  const Int_t icol = fColIndex[index];
1787  if (icol+this->fColLwb <= col_upb && icol+this->fColLwb >= col_lwb)
1788  nr_nonzeros++;
1789  }
1790  }
1791 
1792  target.ResizeTo(row_lwb_sub,row_upb_sub,col_lwb_sub,col_upb_sub,nr_nonzeros);
1793 
1794  const Element *ep = this->GetMatrixArray();
1795 
1796  Int_t *rowIndex_sub = target.GetRowIndexArray();
1797  Int_t *colIndex_sub = target.GetColIndexArray();
1798  Element *ep_sub = target.GetMatrixArray();
1799 
1800  if (target.GetRowIndexArray() && target.GetColIndexArray()) {
1801  Int_t nelems_copy = 0;
1802  rowIndex_sub[0] = 0;
1803  const Int_t row_off = this->fRowLwb-row_lwb;
1804  const Int_t col_off = this->fColLwb-col_lwb;
1805  for (Int_t irow = 0; irow < this->fNrows; irow++) {
1806  if (irow+this->fRowLwb > row_upb || irow+this->fRowLwb < row_lwb) continue;
1807  const Int_t sIndex = fRowIndex[irow];
1808  const Int_t eIndex = fRowIndex[irow+1];
1809  for (Int_t index = sIndex; index < eIndex; index++) {
1810  const Int_t icol = fColIndex[index];
1811  if (icol+this->fColLwb <= col_upb && icol+this->fColLwb >= col_lwb) {
1812  rowIndex_sub[irow+row_off+1] = nelems_copy+1;
1813  colIndex_sub[nelems_copy] = icol+col_off;
1814  ep_sub[nelems_copy] = ep[index];
1815  nelems_copy++;
1816  }
1817  }
1818  }
1819  } else {
1820  const Int_t row_off = this->fRowLwb-row_lwb;
1821  const Int_t col_off = this->fColLwb-col_lwb;
1822  const Int_t ncols_sub = col_upb_sub-col_lwb_sub+1;
1823  for (Int_t irow = 0; irow < this->fNrows; irow++) {
1824  if (irow+this->fRowLwb > row_upb || irow+this->fRowLwb < row_lwb) continue;
1825  const Int_t sIndex = fRowIndex[irow];
1826  const Int_t eIndex = fRowIndex[irow+1];
1827  const Int_t off = (irow+row_off)*ncols_sub;
1828  for (Int_t index = sIndex; index < eIndex; index++) {
1829  const Int_t icol = fColIndex[index];
1830  if (icol+this->fColLwb <= col_upb && icol+this->fColLwb >= col_lwb)
1831  ep_sub[off+icol+col_off] = ep[index];
1832  }
1833  }
1834  }
1835 
1836  return target;
1837 }
1838 
1839 ////////////////////////////////////////////////////////////////////////////////
1840 /// Insert matrix source starting at [row_lwb][col_lwb], thereby overwriting the part
1841 /// [row_lwb..row_lwb+nrows_source-1][col_lwb..col_lwb+ncols_source-1];
1842 
1843 template<class Element>
1844 TMatrixTBase<Element> &TMatrixTSparse<Element>::SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source)
1845 {
1846  if (gMatrixCheck) {
1847  R__ASSERT(this->IsValid());
1848  R__ASSERT(source.IsValid());
1849 
1850  if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
1851  Error("SetSub","row_lwb out-of-bounds");
1852  return *this;
1853  }
1854  if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
1855  Error("SetSub","col_lwb out-of-bounds");
1856  return *this;
1857  }
1858  if (row_lwb+source.GetNrows() > this->fRowLwb+this->fNrows || col_lwb+source.GetNcols() > this->fColLwb+this->fNcols) {
1859  Error("SetSub","source matrix too large");
1860  return *this;
1861  }
1862  }
1863 
1864  const Int_t nRows_source = source.GetNrows();
1865  const Int_t nCols_source = source.GetNcols();
1866 
1867  // Determine how many non-zero's are already available in
1868  // [row_lwb..row_lwb+nrows_source-1][col_lwb..col_lwb+ncols_source-1]
1869  Int_t nr_nonzeros = 0;
1870  for (Int_t irow = 0; irow < this->fNrows; irow++) {
1871  if (irow+this->fRowLwb >= row_lwb+nRows_source || irow+this->fRowLwb < row_lwb) continue;
1872  const Int_t sIndex = fRowIndex[irow];
1873  const Int_t eIndex = fRowIndex[irow+1];
1874  for (Int_t index = sIndex; index < eIndex; index++) {
1875  const Int_t icol = fColIndex[index];
1876  if (icol+this->fColLwb < col_lwb+nCols_source && icol+this->fColLwb >= col_lwb)
1877  nr_nonzeros++;
1878  }
1879  }
1880 
1881  const Int_t *rowIndex_s = source.GetRowIndexArray();
1882  const Int_t *colIndex_s = source.GetColIndexArray();
1883  const Element *elements_s = source.GetMatrixArray();
1884 
1885  const Int_t nelems_old = this->fNelems;
1886  const Int_t *rowIndex_old = GetRowIndexArray();
1887  const Int_t *colIndex_old = GetColIndexArray();
1888  const Element *elements_old = GetMatrixArray();
1889 
1890  const Int_t nelems_new = nelems_old+source.NonZeros()-nr_nonzeros;
1891  fRowIndex = new Int_t[this->fNrowIndex];
1892  fColIndex = new Int_t[nelems_new];
1893  fElements = new Element[nelems_new];
1894  this->fNelems = nelems_new;
1895 
1896  Int_t *rowIndex_new = GetRowIndexArray();
1897  Int_t *colIndex_new = GetColIndexArray();
1898  Element *elements_new = GetMatrixArray();
1899 
1900  const Int_t row_off = row_lwb-this->fRowLwb;
1901  const Int_t col_off = col_lwb-this->fColLwb;
1902 
1903  Int_t nr = 0;
1904  rowIndex_new[0] = 0;
1905  for (Int_t irow = 0; irow < this->fNrows; irow++) {
1906  rowIndex_new[irow+1] = rowIndex_new[irow];
1907  Bool_t flagRow = kFALSE;
1908  if (irow+this->fRowLwb < row_lwb+nRows_source && irow+this->fRowLwb >= row_lwb)
1909  flagRow = kTRUE;
1910 
1911  const Int_t sIndex_o = rowIndex_old[irow];
1912  const Int_t eIndex_o = rowIndex_old[irow+1];
1913 
1914  if (flagRow) {
1915  const Int_t icol_left = col_off-1;
1916  const Int_t left = TMath::BinarySearch(eIndex_o-sIndex_o,colIndex_old+sIndex_o,icol_left)+sIndex_o;
1917  for (Int_t index = sIndex_o; index <= left; index++) {
1918  rowIndex_new[irow+1]++;
1919  colIndex_new[nr] = colIndex_old[index];
1920  elements_new[nr] = elements_old[index];
1921  nr++;
1922  }
1923 
1924  if (rowIndex_s && colIndex_s) {
1925  const Int_t sIndex_s = rowIndex_s[irow-row_off];
1926  const Int_t eIndex_s = rowIndex_s[irow-row_off+1];
1927  for (Int_t index = sIndex_s; index < eIndex_s; index++) {
1928  rowIndex_new[irow+1]++;
1929  colIndex_new[nr] = colIndex_s[index]+col_off;
1930  elements_new[nr] = elements_s[index];
1931  nr++;
1932  }
1933  } else {
1934  const Int_t off = (irow-row_off)*nCols_source;
1935  for (Int_t icol = 0; icol < nCols_source; icol++) {
1936  rowIndex_new[irow+1]++;
1937  colIndex_new[nr] = icol+col_off;
1938  elements_new[nr] = elements_s[off+icol];
1939  nr++;
1940  }
1941  }
1942 
1943  const Int_t icol_right = col_off+nCols_source-1;
1944  if (colIndex_old) {
1945  Int_t right = TMath::BinarySearch(eIndex_o-sIndex_o,colIndex_old+sIndex_o,icol_right)+sIndex_o;
1946  while (right < eIndex_o-1 && colIndex_old[right+1] <= icol_right)
1947  right++;
1948  right++;
1949 
1950  for (Int_t index = right; index < eIndex_o; index++) {
1951  rowIndex_new[irow+1]++;
1952  colIndex_new[nr] = colIndex_old[index];
1953  elements_new[nr] = elements_old[index];
1954  nr++;
1955  }
1956  }
1957  } else {
1958  for (Int_t index = sIndex_o; index < eIndex_o; index++) {
1959  rowIndex_new[irow+1]++;
1960  colIndex_new[nr] = colIndex_old[index];
1961  elements_new[nr] = elements_old[index];
1962  nr++;
1963  }
1964  }
1965  }
1966 
1967  R__ASSERT(this->fNelems == fRowIndex[this->fNrowIndex-1]);
1968 
1969  if (rowIndex_old) delete [] (Int_t*) rowIndex_old;
1970  if (colIndex_old) delete [] (Int_t*) colIndex_old;
1971  if (elements_old) delete [] (Element*) elements_old;
1972 
1973  return *this;
1974 }
1975 
1976 ////////////////////////////////////////////////////////////////////////////////
1977 /// Transpose a matrix.
1978 
1979 template<class Element>
1980 TMatrixTSparse<Element> &TMatrixTSparse<Element>::Transpose(const TMatrixTSparse<Element> &source)
1981 {
1982  if (gMatrixCheck) {
1983  R__ASSERT(this->IsValid());
1984  R__ASSERT(source.IsValid());
1985 
1986  if (this->fNrows != source.GetNcols() || this->fNcols != source.GetNrows() ||
1987  this->fRowLwb != source.GetColLwb() || this->fColLwb != source.GetRowLwb()) {
1988  Error("Transpose","matrix has wrong shape");
1989  return *this;
1990  }
1991 
1992  if (source.NonZeros() <= 0)
1993  return *this;
1994  }
1995 
1996  const Int_t nr_nonzeros = source.NonZeros();
1997 
1998  const Int_t * const pRowIndex_s = source.GetRowIndexArray();
1999  const Int_t * const pColIndex_s = source.GetColIndexArray();
2000  const Element * const pData_s = source.GetMatrixArray();
2001 
2002  Int_t *rownr = new Int_t [nr_nonzeros];
2003  Int_t *colnr = new Int_t [nr_nonzeros];
2004  Element *pData_t = new Element[nr_nonzeros];
2005 
2006  Int_t ielem = 0;
2007  for (Int_t irow_s = 0; irow_s < source.GetNrows(); irow_s++) {
2008  const Int_t sIndex = pRowIndex_s[irow_s];
2009  const Int_t eIndex = pRowIndex_s[irow_s+1];
2010  for (Int_t index = sIndex; index < eIndex; index++) {
2011  rownr[ielem] = pColIndex_s[index]+this->fRowLwb;
2012  colnr[ielem] = irow_s+this->fColLwb;
2013  pData_t[ielem] = pData_s[index];
2014  ielem++;
2015  }
2016  }
2017 
2018  R__ASSERT(nr_nonzeros >= ielem);
2019 
2020  TMatrixTBase<Element>::DoubleLexSort(nr_nonzeros,rownr,colnr,pData_t);
2021  SetMatrixArray(nr_nonzeros,rownr,colnr,pData_t);
2022 
2023  R__ASSERT(this->fNelems == fRowIndex[this->fNrowIndex-1]);
2024 
2025  if (pData_t) delete [] pData_t;
2026  if (rownr) delete [] rownr;
2027  if (colnr) delete [] colnr;
2028 
2029  return *this;
2030 }
2031 
2032 ////////////////////////////////////////////////////////////////////////////////
2033 
2034 template<class Element>
2035 TMatrixTBase<Element> &TMatrixTSparse<Element>::Zero()
2036 {
2037  R__ASSERT(this->IsValid());
2038 
2039  if (fElements) { delete [] fElements; fElements = 0; }
2040  if (fColIndex) { delete [] fColIndex; fColIndex = 0; }
2041  this->fNelems = 0;
2042  memset(this->GetRowIndexArray(),0,this->fNrowIndex*sizeof(Int_t));
2043 
2044  return *this;
2045 }
2046 
2047 ////////////////////////////////////////////////////////////////////////////////
2048 /// Make a unit matrix (matrix need not be a square one).
2049 
2050 template<class Element>
2051 TMatrixTBase<Element> &TMatrixTSparse<Element>::UnitMatrix()
2052 {
2053  R__ASSERT(this->IsValid());
2054 
2055  Int_t i;
2056 
2057  Int_t nr_nonzeros = 0;
2058  for (i = this->fRowLwb; i <= this->fRowLwb+this->fNrows-1; i++)
2059  for (Int_t j = this->fColLwb; j <= this->fColLwb+this->fNcols-1; j++)
2060  if (i == j) nr_nonzeros++;
2061 
2062  if (nr_nonzeros != this->fNelems) {
2063  this->fNelems = nr_nonzeros;
2064  Int_t *oIp = fColIndex;
2065  fColIndex = new Int_t[nr_nonzeros];
2066  if (oIp) delete [] oIp;
2067  Element *oDp = fElements;
2068  fElements = new Element[nr_nonzeros];
2069  if (oDp) delete [] oDp;
2070  }
2071 
2072  Int_t ielem = 0;
2073  fRowIndex[0] = 0;
2074  for (i = this->fRowLwb; i <= this->fRowLwb+this->fNrows-1; i++) {
2075  for (Int_t j = this->fColLwb; j <= this->fColLwb+this->fNcols-1; j++) {
2076  if (i == j) {
2077  const Int_t irow = i-this->fRowLwb;
2078  fRowIndex[irow+1] = ielem+1;
2079  fElements[ielem] = 1.0;
2080  fColIndex[ielem++] = j-this->fColLwb;
2081  }
2082  }
2083  }
2084 
2085  return *this;
2086 }
2087 
2088 ////////////////////////////////////////////////////////////////////////////////
2089 /// Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
2090 /// The norm is induced by the infinity vector norm.
2091 
2092 template<class Element>
2093 Element TMatrixTSparse<Element>::RowNorm() const
2094 {
2095  R__ASSERT(this->IsValid());
2096 
2097  const Element * ep = GetMatrixArray();
2098  const Element * const fp = ep+this->fNelems;
2099  const Int_t * const pR = GetRowIndexArray();
2100  Element norm = 0;
2101 
2102  // Scan the matrix row-after-row
2103  for (Int_t irow = 0; irow < this->fNrows; irow++) {
2104  const Int_t sIndex = pR[irow];
2105  const Int_t eIndex = pR[irow+1];
2106  Element sum = 0;
2107  for (Int_t index = sIndex; index < eIndex; index++)
2108  sum += TMath::Abs(*ep++);
2109  norm = TMath::Max(norm,sum);
2110  }
2111 
2112  R__ASSERT(ep == fp);
2113 
2114  return norm;
2115 }
2116 
2117 ////////////////////////////////////////////////////////////////////////////////
2118 /// Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
2119 /// The norm is induced by the 1 vector norm.
2120 
2121 template<class Element>
2122 Element TMatrixTSparse<Element>::ColNorm() const
2123 {
2124  R__ASSERT(this->IsValid());
2125 
2126  const TMatrixTSparse<Element> mt(kTransposed,*this);
2127 
2128  const Element * ep = mt.GetMatrixArray();
2129  const Element * const fp = ep+this->fNcols;
2130  Element norm = 0;
2131 
2132  // Scan the matrix col-after-col
2133  while (ep < fp) {
2134  Element sum = 0;
2135  // Scan a col to compute the sum
2136  for (Int_t i = 0; i < this->fNrows; i++,ep += this->fNcols)
2137  sum += TMath::Abs(*ep);
2138  ep -= this->fNelems-1; // Point ep to the beginning of the next col
2139  norm = TMath::Max(norm,sum);
2140  }
2141 
2142  R__ASSERT(ep == fp);
2143 
2144  return norm;
2145 }
2146 
2147 ////////////////////////////////////////////////////////////////////////////////
2148 
2149 template<class Element>
2150 Element &TMatrixTSparse<Element>::operator()(Int_t rown,Int_t coln)
2151 {
2152  R__ASSERT(this->IsValid());
2153 
2154  const Int_t arown = rown-this->fRowLwb;
2155  const Int_t acoln = coln-this->fColLwb;
2156  if (arown >= this->fNrows || arown < 0) {
2157  Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
2158  return fElements[0];
2159  }
2160  if (acoln >= this->fNcols || acoln < 0) {
2161  Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
2162  return fElements[0];
2163  }
2164 
2165  Int_t index = -1;
2166  Int_t sIndex = 0;
2167  Int_t eIndex = 0;
2168  if (this->fNrowIndex > 0 && fRowIndex[this->fNrowIndex-1] != 0) {
2169  sIndex = fRowIndex[arown];
2170  eIndex = fRowIndex[arown+1];
2171  index = TMath::BinarySearch(eIndex-sIndex,fColIndex+sIndex,acoln)+sIndex;
2172  }
2173 
2174  if (index >= sIndex && fColIndex[index] == acoln)
2175  return fElements[index];
2176  else {
2177  Element val = 0.;
2178  InsertRow(rown,coln,&val,1);
2179  sIndex = fRowIndex[arown];
2180  eIndex = fRowIndex[arown+1];
2181  index = TMath::BinarySearch(eIndex-sIndex,fColIndex+sIndex,acoln)+sIndex;
2182  if (index >= sIndex && fColIndex[index] == acoln)
2183  return fElements[index];
2184  else {
2185  Error("operator()(Int_t,Int_t","Insert row failed");
2186  return fElements[0];
2187  }
2188  }
2189 }
2190 
2191 ////////////////////////////////////////////////////////////////////////////////
2192 
2193 template <class Element>
2194 Element TMatrixTSparse<Element>::operator()(Int_t rown,Int_t coln) const
2195 {
2196  R__ASSERT(this->IsValid());
2197  if (this->fNrowIndex > 0 && this->fRowIndex[this->fNrowIndex-1] == 0) {
2198  Error("operator()(Int_t,Int_t) const","row/col indices are not set");
2199  Info("operator()","fNrowIndex = %d fRowIndex[fNrowIndex-1] = %d\n",this->fNrowIndex,this->fRowIndex[this->fNrowIndex-1]);
2200  return 0.0;
2201  }
2202  const Int_t arown = rown-this->fRowLwb;
2203  const Int_t acoln = coln-this->fColLwb;
2204 
2205  if (arown >= this->fNrows || arown < 0) {
2206  Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
2207  return 0.0;
2208  }
2209  if (acoln >= this->fNcols || acoln < 0) {
2210  Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
2211  return 0.0;
2212  }
2213 
2214  const Int_t sIndex = fRowIndex[arown];
2215  const Int_t eIndex = fRowIndex[arown+1];
2216  const Int_t index = (Int_t)TMath::BinarySearch(eIndex-sIndex,fColIndex+sIndex,acoln)+sIndex;
2217  if (index < sIndex || fColIndex[index] != acoln) return 0.0;
2218  else return fElements[index];
2219 }
2220 
2221 ////////////////////////////////////////////////////////////////////////////////
2222 /// Notice that the sparsity of the matrix is NOT changed : its fRowIndex/fColIndex
2223 /// are used !
2224 
2225 template<class Element>
2226 TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator=(const TMatrixTSparse<Element> &source)
2227 {
2228  if (gMatrixCheck && !AreCompatible(*this,source)) {
2229  Error("operator=(const TMatrixTSparse &)","matrices not compatible");
2230  return *this;
2231  }
2232 
2233  if (this->GetMatrixArray() != source.GetMatrixArray()) {
2234  TObject::operator=(source);
2235 
2236  const Element * const sp = source.GetMatrixArray();
2237  Element * const tp = this->GetMatrixArray();
2238  memcpy(tp,sp,this->fNelems*sizeof(Element));
2239  this->fTol = source.GetTol();
2240  }
2241  return *this;
2242 }
2243 
2244 ////////////////////////////////////////////////////////////////////////////////
2245 /// Notice that the sparsity of the matrix is NOT changed : its fRowIndex/fColIndex
2246 /// are used !
2247 
2248 template<class Element>
2249 TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator=(const TMatrixT<Element> &source)
2250 {
2251  if (gMatrixCheck && !AreCompatible(*this,(TMatrixTBase<Element> &)source)) {
2252  Error("operator=(const TMatrixT &)","matrices not compatible");
2253  return *this;
2254  }
2255 
2256  if (this->GetMatrixArray() != source.GetMatrixArray()) {
2257  TObject::operator=(source);
2258 
2259  const Element * const sp = source.GetMatrixArray();
2260  Element * const tp = this->GetMatrixArray();
2261  for (Int_t irow = 0; irow < this->fNrows; irow++) {
2262  const Int_t sIndex = fRowIndex[irow];
2263  const Int_t eIndex = fRowIndex[irow+1];
2264  const Int_t off = irow*this->fNcols;
2265  for (Int_t index = sIndex; index < eIndex; index++) {
2266  const Int_t icol = fColIndex[index];
2267  tp[index] = sp[off+icol];
2268  }
2269  }
2270  this->fTol = source.GetTol();
2271  }
2272  return *this;
2273 }
2274 
2275 ////////////////////////////////////////////////////////////////////////////////
2276 /// Assign val to every element of the matrix. Check that the row/col
2277 /// indices are set !
2278 
2279 template<class Element>
2280 TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator=(Element val)
2281 {
2282  R__ASSERT(this->IsValid());
2283 
2284  if (fRowIndex[this->fNrowIndex-1] == 0) {
2285  Error("operator=(Element","row/col indices are not set");
2286  return *this;
2287  }
2288 
2289  Element *ep = this->GetMatrixArray();
2290  const Element * const ep_last = ep+this->fNelems;
2291  while (ep < ep_last)
2292  *ep++ = val;
2293 
2294  return *this;
2295 }
2296 
2297 ////////////////////////////////////////////////////////////////////////////////
2298 /// Add val to every element of the matrix.
2299 
2300 template<class Element>
2301 TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator+=(Element val)
2302 {
2303  R__ASSERT(this->IsValid());
2304 
2305  Element *ep = this->GetMatrixArray();
2306  const Element * const ep_last = ep+this->fNelems;
2307  while (ep < ep_last)
2308  *ep++ += val;
2309 
2310  return *this;
2311 }
2312 
2313 ////////////////////////////////////////////////////////////////////////////////
2314 /// Subtract val from every element of the matrix.
2315 
2316 template<class Element>
2317 TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator-=(Element val)
2318 {
2319  R__ASSERT(this->IsValid());
2320 
2321  Element *ep = this->GetMatrixArray();
2322  const Element * const ep_last = ep+this->fNelems;
2323  while (ep < ep_last)
2324  *ep++ -= val;
2325 
2326  return *this;
2327 }
2328 
2329 ////////////////////////////////////////////////////////////////////////////////
2330 /// Multiply every element of the matrix with val.
2331 
2332 template<class Element>
2333 TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator*=(Element val)
2334 {
2335  R__ASSERT(this->IsValid());
2336 
2337  Element *ep = this->GetMatrixArray();
2338  const Element * const ep_last = ep+this->fNelems;
2339  while (ep < ep_last)
2340  *ep++ *= val;
2341 
2342  return *this;
2343 }
2344 
2345 ////////////////////////////////////////////////////////////////////////////////
2346 /// randomize matrix element values
2347 
2348 template<class Element>
2349 TMatrixTBase<Element> &TMatrixTSparse<Element>::Randomize(Element alpha,Element beta,Double_t &seed)
2350 {
2351  R__ASSERT(this->IsValid());
2352 
2353  const Element scale = beta-alpha;
2354  const Element shift = alpha/scale;
2355 
2356  Int_t * const pRowIndex = GetRowIndexArray();
2357  Int_t * const pColIndex = GetColIndexArray();
2358  Element * const ep = GetMatrixArray();
2359 
2360  const Int_t m = this->GetNrows();
2361  const Int_t n = this->GetNcols();
2362 
2363  // Knuth's algorithm for choosing "length" elements out of nn .
2364  const Int_t nn = this->GetNrows()*this->GetNcols();
2365  const Int_t length = (this->GetNoElements() <= nn) ? this->GetNoElements() : nn;
2366  Int_t chosen = 0;
2367  Int_t icurrent = 0;
2368  pRowIndex[0] = 0;
2369  for (Int_t k = 0; k < nn; k++) {
2370  const Element r = Drand(seed);
2371 
2372  if ((nn-k)*r < length-chosen) {
2373  pColIndex[chosen] = k%n;
2374  const Int_t irow = k/n;
2375 
2376  if (irow > icurrent) {
2377  for ( ; icurrent < irow; icurrent++)
2378  pRowIndex[icurrent+1] = chosen;
2379  }
2380  ep[chosen] = scale*(Drand(seed)+shift);
2381  chosen++;
2382  }
2383  }
2384  for ( ; icurrent < m; icurrent++)
2385  pRowIndex[icurrent+1] = length;
2386 
2387  R__ASSERT(chosen == length);
2388 
2389  return *this;
2390 }
2391 
2392 ////////////////////////////////////////////////////////////////////////////////
2393 /// randomize matrix element values but keep matrix symmetric positive definite
2394 
2395 template<class Element>
2396 TMatrixTSparse<Element> &TMatrixTSparse<Element>::RandomizePD(Element alpha,Element beta,Double_t &seed)
2397 {
2398  const Element scale = beta-alpha;
2399  const Element shift = alpha/scale;
2400 
2401  if (gMatrixCheck) {
2402  R__ASSERT(this->IsValid());
2403 
2404  if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
2405  Error("RandomizePD(Element &","matrix should be square");
2406  return *this;
2407  }
2408  }
2409 
2410  const Int_t n = this->fNcols;
2411 
2412  Int_t * const pRowIndex = this->GetRowIndexArray();
2413  Int_t * const pColIndex = this->GetColIndexArray();
2414  Element * const ep = GetMatrixArray();
2415 
2416  // We will always have non-zeros on the diagonal, so there
2417  // is no randomness there. In fact, choose the (0,0) element now
2418  pRowIndex[0] = 0;
2419  pColIndex[0] = 0;
2420  pRowIndex[1] = 1;
2421  ep[0] = 1e-8+scale*(Drand(seed)+shift);
2422 
2423  // Knuth's algorithm for choosing length elements out of nn .
2424  // nn here is the number of elements in the strict lower triangle.
2425  const Int_t nn = n*(n-1)/2;
2426 
2427  // length is the number of elements that can be stored, minus the number
2428  // of elements in the diagonal, which will always be in the matrix.
2429  // Int_t length = (this->fNelems-n)/2;
2430  Int_t length = this->fNelems-n;
2431  length = (length <= nn) ? length : nn;
2432 
2433  // chosen : the number of elements that have already been chosen (now 0)
2434  // nnz : the number of non-zeros in the matrix (now 1, because the
2435  // (0,0) element is already in the matrix.
2436  // icurrent : the index of the last row whose start has been stored in pRowIndex;
2437 
2438  Int_t chosen = 0;
2439  Int_t icurrent = 1;
2440  Int_t nnz = 1;
2441  for (Int_t k = 0; k < nn; k++ ) {
2442  const Element r = Drand(seed);
2443 
2444  if( (nn-k)*r < length-chosen) {
2445  // Element k is chosen. What row is it in?
2446  // In a lower triangular matrix (including a diagonal), it will be in
2447  // the largest row such that row*(row+1)/2 < k. In other words
2448 
2449  Int_t row = (int) TMath::Floor((-1+TMath::Sqrt(1.0+8.0*k))/2);
2450  // and its column will be the remainder
2451  Int_t col = k-row*(row+1)/2;
2452  // but since we are only filling in the *strict* lower triangle of
2453  // the matrix, we shift the row by 1
2454  row++;
2455 
2456  if (row > icurrent) {
2457  // We have chosen a row beyond the current row.
2458  // Choose a diagonal element for each intermediate row and fix the
2459  // data structure.
2460  for ( ; icurrent < row; icurrent++) {
2461  // Choose the diagonal
2462  ep[nnz] = 0.0;
2463  for (Int_t ll = pRowIndex[icurrent]; ll < nnz; ll++)
2464  ep[nnz] += TMath::Abs(ep[ll]);
2465  ep[nnz] += 1e-8+scale*(Drand(seed)+shift);
2466  pColIndex[nnz] = icurrent;
2467 
2468  nnz++;
2469  pRowIndex[icurrent+1] = nnz;
2470  }
2471  } // end if we have chosen a row beyond the current row;
2472  ep[nnz] = scale*(Drand(seed)+shift);
2473  pColIndex[nnz] = col;
2474  // add the value of this element (which occurs symmetrically in the
2475  // upper triangle) to the appropriate diagonal element
2476  ep[pRowIndex[col+1]-1] += TMath::Abs(ep[nnz]);
2477 
2478  nnz++; // We have added another element to the matrix
2479  chosen++; // And finished choosing another element.
2480  }
2481  }
2482 
2483  R__ASSERT(chosen == length);
2484 
2485  // and of course, we must choose all remaining diagonal elements .
2486  for ( ; icurrent < n; icurrent++) {
2487  // Choose the diagonal
2488  ep[nnz] = 0.0;
2489  for(Int_t ll = pRowIndex[icurrent]; ll < nnz; ll++)
2490  ep[nnz] += TMath::Abs(ep[ll]);
2491  ep[nnz] += 1e-8+scale*(Drand(seed)+shift);
2492  pColIndex[nnz] = icurrent;
2493 
2494  nnz++;
2495  pRowIndex[icurrent+1] = nnz;
2496  }
2497  this->fNelems = nnz;
2498 
2499  TMatrixTSparse<Element> tmp(TMatrixTSparse<Element>::kTransposed,*this);
2500  *this += tmp;
2501 
2502  // make sure to divide the diagonal by 2 becuase the operation
2503  // *this += tmp; adds the diagonal again
2504  {
2505  const Int_t * const pR = this->GetRowIndexArray();
2506  const Int_t * const pC = this->GetColIndexArray();
2507  Element * const pD = GetMatrixArray();
2508  for (Int_t irow = 0; irow < this->fNrows+1; irow++) {
2509  const Int_t sIndex = pR[irow];
2510  const Int_t eIndex = pR[irow+1];
2511  for (Int_t index = sIndex; index < eIndex; index++) {
2512  const Int_t icol = pC[index];
2513  if (irow == icol)
2514  pD[index] /= 2.;
2515  }
2516  }
2517  }
2518 
2519  return *this;
2520 }
2521 
2522 ////////////////////////////////////////////////////////////////////////////////
2523 
2524 template<class Element>
2525 TMatrixTSparse<Element> operator+(const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2)
2526 {
2527  TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kPlus,source2);
2528  return target;
2529 }
2530 
2531 ////////////////////////////////////////////////////////////////////////////////
2532 
2533 template<class Element>
2534 TMatrixTSparse<Element> operator+(const TMatrixTSparse<Element> &source1,const TMatrixT<Element> &source2)
2535 {
2536  TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kPlus,source2);
2537  return target;
2538 }
2539 
2540 ////////////////////////////////////////////////////////////////////////////////
2541 
2542 template<class Element>
2543 TMatrixTSparse<Element> operator+(const TMatrixT<Element> &source1,const TMatrixTSparse<Element> &source2)
2544 {
2545  TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kPlus,source2);
2546  return target;
2547 }
2548 
2549 ////////////////////////////////////////////////////////////////////////////////
2550 
2551 template<class Element>
2552 TMatrixTSparse<Element> operator+(const TMatrixTSparse<Element> &source,Element val)
2553 {
2554  TMatrixTSparse<Element> target(source);
2555  target += val;
2556  return target;
2557 }
2558 
2559 ////////////////////////////////////////////////////////////////////////////////
2560 
2561 template<class Element>
2562 TMatrixTSparse<Element> operator+(Element val,const TMatrixTSparse<Element> &source)
2563 {
2564  TMatrixTSparse<Element> target(source);
2565  target += val;
2566  return target;
2567 }
2568 
2569 ////////////////////////////////////////////////////////////////////////////////
2570 
2571 template<class Element>
2572 TMatrixTSparse<Element> operator-(const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2)
2573 {
2574  TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMinus,source2);
2575  return target;
2576 }
2577 
2578 ////////////////////////////////////////////////////////////////////////////////
2579 
2580 template<class Element>
2581 TMatrixTSparse<Element> operator-(const TMatrixTSparse<Element> &source1,const TMatrixT<Element> &source2)
2582 {
2583  TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMinus,source2);
2584  return target;
2585 }
2586 
2587 ////////////////////////////////////////////////////////////////////////////////
2588 
2589 template<class Element>
2590 TMatrixTSparse<Element> operator-(const TMatrixT<Element> &source1,const TMatrixTSparse<Element> &source2)
2591 {
2592  TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMinus,source2);
2593  return target;
2594 }
2595 
2596 ////////////////////////////////////////////////////////////////////////////////
2597 
2598 template<class Element>
2599 TMatrixTSparse<Element> operator-(const TMatrixTSparse<Element> &source,Element val)
2600 {
2601  TMatrixTSparse<Element> target(source);
2602  target -= val;
2603  return target;
2604 }
2605 
2606 ////////////////////////////////////////////////////////////////////////////////
2607 
2608 template<class Element>
2609 TMatrixTSparse<Element> operator-(Element val,const TMatrixTSparse<Element> &source)
2610 {
2611  TMatrixTSparse<Element> target(source);
2612  target -= val;
2613  return target;
2614 }
2615 
2616 ////////////////////////////////////////////////////////////////////////////////
2617 
2618 template<class Element>
2619 TMatrixTSparse<Element> operator*(const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2)
2620 {
2621  TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMult,source2);
2622  return target;
2623 }
2624 
2625 ////////////////////////////////////////////////////////////////////////////////
2626 
2627 template<class Element>
2628 TMatrixTSparse<Element> operator*(const TMatrixTSparse<Element> &source1,const TMatrixT<Element> &source2)
2629 {
2630  TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMult,source2);
2631  return target;
2632 }
2633 
2634 ////////////////////////////////////////////////////////////////////////////////
2635 
2636 template<class Element>
2637 TMatrixTSparse<Element> operator*(const TMatrixT<Element> &source1,const TMatrixTSparse<Element> &source2)
2638 {
2639  TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMult,source2);
2640  return target;
2641 }
2642 
2643 ////////////////////////////////////////////////////////////////////////////////
2644 
2645 template<class Element>
2646 TMatrixTSparse<Element> operator*(Element val,const TMatrixTSparse<Element> &source)
2647 {
2648  TMatrixTSparse<Element> target(source);
2649  target *= val;
2650  return target;
2651 }
2652 
2653 ////////////////////////////////////////////////////////////////////////////////
2654 
2655 template<class Element>
2656 TMatrixTSparse<Element> operator*(const TMatrixTSparse<Element> &source,Element val)
2657 {
2658  TMatrixTSparse<Element> target(source);
2659  target *= val;
2660  return target;
2661 }
2662 
2663 ////////////////////////////////////////////////////////////////////////////////
2664 /// Modify addition: target += scalar * source.
2665 
2666 template<class Element>
2667 TMatrixTSparse<Element> &Add(TMatrixTSparse<Element> &target,Element scalar,const TMatrixTSparse<Element> &source)
2668 {
2669  target += scalar * source;
2670  return target;
2671 }
2672 
2673 ////////////////////////////////////////////////////////////////////////////////
2674 /// Multiply target by the source, element-by-element.
2675 
2676 template<class Element>
2677 TMatrixTSparse<Element> &ElementMult(TMatrixTSparse<Element> &target,const TMatrixTSparse<Element> &source)
2678 {
2679  if (gMatrixCheck && !AreCompatible(target,source)) {
2680  ::Error("ElementMult(TMatrixTSparse &,const TMatrixTSparse &)","matrices not compatible");
2681  return target;
2682  }
2683 
2684  const Element *sp = source.GetMatrixArray();
2685  Element *tp = target.GetMatrixArray();
2686  const Element *ftp = tp+target.GetNoElements();
2687  while ( tp < ftp )
2688  *tp++ *= *sp++;
2689 
2690  return target;
2691 }
2692 
2693 ////////////////////////////////////////////////////////////////////////////////
2694 /// Divide target by the source, element-by-element.
2695 
2696 template<class Element>
2697 TMatrixTSparse<Element> &ElementDiv(TMatrixTSparse<Element> &target,const TMatrixTSparse<Element> &source)
2698 {
2699  if (gMatrixCheck && !AreCompatible(target,source)) {
2700  ::Error("ElementDiv(TMatrixT &,const TMatrixT &)","matrices not compatible");
2701  return target;
2702  }
2703 
2704  const Element *sp = source.GetMatrixArray();
2705  Element *tp = target.GetMatrixArray();
2706  const Element *ftp = tp+target.GetNoElements();
2707  while ( tp < ftp ) {
2708  if (*sp != 0.0)
2709  *tp++ /= *sp++;
2710  else {
2711  Error("ElementDiv","source element is zero");
2712  tp++;
2713  }
2714  }
2715 
2716  return target;
2717 }
2718 
2719 ////////////////////////////////////////////////////////////////////////////////
2720 
2721 template<class Element>
2722 Bool_t AreCompatible(const TMatrixTSparse<Element> &m1,const TMatrixTSparse<Element> &m2,Int_t verbose)
2723 {
2724  if (!m1.IsValid()) {
2725  if (verbose)
2726  ::Error("AreCompatible", "matrix 1 not valid");
2727  return kFALSE;
2728  }
2729  if (!m2.IsValid()) {
2730  if (verbose)
2731  ::Error("AreCompatible", "matrix 2 not valid");
2732  return kFALSE;
2733  }
2734 
2735  if (m1.GetNrows() != m2.GetNrows() || m1.GetNcols() != m2.GetNcols() ||
2736  m1.GetRowLwb() != m2.GetRowLwb() || m1.GetColLwb() != m2.GetColLwb()) {
2737  if (verbose)
2738  ::Error("AreCompatible", "matrices 1 and 2 not compatible");
2739  return kFALSE;
2740  }
2741 
2742  const Int_t *pR1 = m1.GetRowIndexArray();
2743  const Int_t *pR2 = m2.GetRowIndexArray();
2744  const Int_t nRows = m1.GetNrows();
2745  if (memcmp(pR1,pR2,(nRows+1)*sizeof(Int_t))) {
2746  if (verbose)
2747  ::Error("AreCompatible", "matrices 1 and 2 have different rowIndex");
2748  for (Int_t i = 0; i < nRows+1; i++)
2749  printf("%d: %d %d\n",i,pR1[i],pR2[i]);
2750  return kFALSE;
2751  }
2752  const Int_t *pD1 = m1.GetColIndexArray();
2753  const Int_t *pD2 = m2.GetColIndexArray();
2754  const Int_t nData = m1.GetNoElements();
2755  if (memcmp(pD1,pD2,nData*sizeof(Int_t))) {
2756  if (verbose)
2757  ::Error("AreCompatible", "matrices 1 and 2 have different colIndex");
2758  for (Int_t i = 0; i < nData; i++)
2759  printf("%d: %d %d\n",i,pD1[i],pD2[i]);
2760  return kFALSE;
2761  }
2762 
2763  return kTRUE;
2764 }
2765 
2766 ////////////////////////////////////////////////////////////////////////////////
2767 /// Stream an object of class TMatrixTSparse.
2768 
2769 template<class Element>
2770 void TMatrixTSparse<Element>::Streamer(TBuffer &R__b)
2771 {
2772  if (R__b.IsReading()) {
2773  UInt_t R__s, R__c;
2774  Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
2775  Clear();
2776  R__b.ReadClassBuffer(TMatrixTSparse<Element>::Class(),this,R__v,R__s,R__c);
2777  if (this->fNelems < 0)
2778  this->Invalidate();
2779  } else {
2780  R__b.WriteClassBuffer(TMatrixTSparse<Element>::Class(),this);
2781  }
2782 }
2783 
2784 template class TMatrixTSparse<Float_t>;
2785 
2786 #ifndef ROOT_TMatrixFSparsefwd
2787 #include "TMatrixFSparsefwd.h"
2788 #endif
2789 #ifndef ROOT_TMatrixFfwd
2790 #include "TMatrixFfwd.h"
2791 #endif
2792 
2793 template TMatrixFSparse operator+ <Float_t>(const TMatrixFSparse &source1,const TMatrixFSparse &source2);
2794 template TMatrixFSparse operator+ <Float_t>(const TMatrixFSparse &source1,const TMatrixF &source2);
2795 template TMatrixFSparse operator+ <Float_t>(const TMatrixF &source1,const TMatrixFSparse &source2);
2796 template TMatrixFSparse operator+ <Float_t>(const TMatrixFSparse &source , Float_t val );
2797 template TMatrixFSparse operator+ <Float_t>( Float_t val ,const TMatrixFSparse &source );
2798 template TMatrixFSparse operator- <Float_t>(const TMatrixFSparse &source1,const TMatrixFSparse &source2);
2799 template TMatrixFSparse operator- <Float_t>(const TMatrixFSparse &source1,const TMatrixF &source2);
2800 template TMatrixFSparse operator- <Float_t>(const TMatrixF &source1,const TMatrixFSparse &source2);
2801 template TMatrixFSparse operator- <Float_t>(const TMatrixFSparse &source , Float_t val );
2802 template TMatrixFSparse operator- <Float_t>( Float_t val ,const TMatrixFSparse &source );
2803 template TMatrixFSparse operator* <Float_t>(const TMatrixFSparse &source1,const TMatrixFSparse &source2);
2804 template TMatrixFSparse operator* <Float_t>(const TMatrixFSparse &source1,const TMatrixF &source2);
2805 template TMatrixFSparse operator* <Float_t>(const TMatrixF &source1,const TMatrixFSparse &source2);
2806 template TMatrixFSparse operator* <Float_t>( Float_t val ,const TMatrixFSparse &source );
2807 template TMatrixFSparse operator* <Float_t>(const TMatrixFSparse &source, Float_t val );
2808 
2809 template TMatrixFSparse &Add <Float_t>(TMatrixFSparse &target, Float_t scalar,const TMatrixFSparse &source);
2810 template TMatrixFSparse &ElementMult <Float_t>(TMatrixFSparse &target,const TMatrixFSparse &source);
2811 template TMatrixFSparse &ElementDiv <Float_t>(TMatrixFSparse &target,const TMatrixFSparse &source);
2812 
2813 template Bool_t AreCompatible<Float_t>(const TMatrixFSparse &m1,const TMatrixFSparse &m2,Int_t verbose);
2814 
2815 #ifndef ROOT_TMatrixDSparsefwd
2816 #include "TMatrixDSparsefwd.h"
2817 #endif
2818 #ifndef ROOT_TMatrixDfwd
2819 #include "TMatrixDfwd.h"
2820 #endif
2821 
2822 template class TMatrixTSparse<Double_t>;
2823 
2824 template TMatrixDSparse operator+ <Double_t>(const TMatrixDSparse &source1,const TMatrixDSparse &source2);
2825 template TMatrixDSparse operator+ <Double_t>(const TMatrixDSparse &source1,const TMatrixD &source2);
2826 template TMatrixDSparse operator+ <Double_t>(const TMatrixD &source1,const TMatrixDSparse &source2);
2827 template TMatrixDSparse operator+ <Double_t>(const TMatrixDSparse &source , Double_t val );
2828 template TMatrixDSparse operator+ <Double_t>( Double_t val ,const TMatrixDSparse &source );
2829 template TMatrixDSparse operator- <Double_t>(const TMatrixDSparse &source1,const TMatrixDSparse &source2);
2830 template TMatrixDSparse operator- <Double_t>(const TMatrixDSparse &source1,const TMatrixD &source2);
2831 template TMatrixDSparse operator- <Double_t>(const TMatrixD &source1,const TMatrixDSparse &source2);
2832 template TMatrixDSparse operator- <Double_t>(const TMatrixDSparse &source , Double_t val );
2833 template TMatrixDSparse operator- <Double_t>( Double_t val ,const TMatrixDSparse &source );
2834 template TMatrixDSparse operator* <Double_t>(const TMatrixDSparse &source1,const TMatrixDSparse &source2);
2835 template TMatrixDSparse operator* <Double_t>(const TMatrixDSparse &source1,const TMatrixD &source2);
2836 template TMatrixDSparse operator* <Double_t>(const TMatrixD &source1,const TMatrixDSparse &source2);
2837 template TMatrixDSparse operator* <Double_t>( Double_t val ,const TMatrixDSparse &source );
2838 template TMatrixDSparse operator* <Double_t>(const TMatrixDSparse &source, Double_t val );
2839 
2840 template TMatrixDSparse &Add <Double_t>(TMatrixDSparse &target, Double_t scalar,const TMatrixDSparse &source);
2841 template TMatrixDSparse &ElementMult <Double_t>(TMatrixDSparse &target,const TMatrixDSparse &source);
2842 template TMatrixDSparse &ElementDiv <Double_t>(TMatrixDSparse &target,const TMatrixDSparse &source);
2843 
2844 template Bool_t AreCompatible<Double_t>(const TMatrixDSparse &m1,const TMatrixDSparse &m2,Int_t verbose);
TMatrixTBase::GetRowLwb
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:132
TMatrixTSparse::AMinusB
void AMinusB(const TMatrixTSparse< Element > &a, const TMatrixTSparse< Element > &b, Int_t constr=0)
General matrix subtraction.
Definition: TMatrixTSparse.cxx:928
mrt.row
tuple row
Definition: mrt.py:26
TMatrixTBase::GetMatrixArray
virtual const Element * GetMatrixArray() const =0
TMatrixTSparse::InsertRow
virtual TMatrixTBase< Element > & InsertRow(Int_t row, Int_t col, const Element *v, Int_t n=-1)
Insert in row rown, n elements of array v at column coln.
Definition: TMatrixTSparse.cxx:360
TMatrixTSparse::GetMatrix2Array
virtual void GetMatrix2Array(Element *data, Option_t *="") const
Copy matrix data to array . It is assumed that array is of size >= fNelems.
Definition: TMatrixTSparse.cxx:1151
TBuffer::WriteClassBuffer
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
TMath::Floor
Double_t Floor(Double_t x)
Definition: TMath.h:473
template
Small helper to encapsulate whether to return the value pointed to by the iterator or its address...
TMath::LocMax
Long64_t LocMax(Long64_t n, const T *a)
Definition: TMath.h:724
TBuffer::IsReading
Bool_t IsReading() const
Definition: TBuffer.h:83
Version_t
short Version_t
Definition: RtypesCore.h:61
TMatrixTSparse::Transpose
TMatrixTSparse< Element > & Transpose(const TMatrixTSparse< Element > &source)
Transpose a matrix.
Definition: TMatrixTSparse.cxx:1980
Float_t
float Float_t
Definition: RtypesCore.h:53
TMatrixTSparse::NonZeros
virtual Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TMatrixTSparse.h:161
TMatrixTSparse::GetMatrixArray
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSparse.h:215
Option_t
const char Option_t
Definition: RtypesCore.h:62
TMatrixTSparse::UnitMatrix
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
Definition: TMatrixTSparse.cxx:2051
Add< Double_t >
template TMatrixDSparse & Add< Double_t >(TMatrixDSparse &target, Double_t scalar, const TMatrixDSparse &source)
TString::ToUpper
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1088
TBuffer
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
ElementMult< Double_t >
template TMatrixDSparse & ElementMult< Double_t >(TMatrixDSparse &target, const TMatrixDSparse &source)
R__ASSERT
#define R__ASSERT(e)
Definition: TError.h:98
TMatrixTBase::GetNoElements
Int_t GetNoElements() const
Definition: TMatrixTBase.h:138
ROOT::Math::Transpose
Expr< TransposeOp< SMatrix< T, D1, D2, R >, T, D1, D2 >, T, D2, D1, typename TranspPolicy< T, D1, D2, R >::RepType > Transpose(const SMatrix< T, D1, D2, R > &rhs)
Matrix Transpose B(i,j) = A(j,i) returning a matrix expression.
Definition: MatrixFunctions.h:546
TString
Basic string class.
Definition: TString.h:137
ElementDiv< Float_t >
template TMatrixFSparse & ElementDiv< Float_t >(TMatrixFSparse &target, const TMatrixFSparse &source)
TMatrixTSparse::SetMatrixArray
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *="")
Copy array data to matrix .
Definition: TMatrixTSparse.h:109
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:170
Int_t
int Int_t
Definition: RtypesCore.h:41
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
a
TArc * a
Definition: textangle.C:12
kFALSE
const Bool_t kFALSE
Definition: Rtypes.h:92
TMatrixTBase::NonZeros
virtual Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TMatrixTBase.cxx:741
TMatrixTBase::ResizeTo
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1)=0
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
at
TObject Int_t at
Definition: TStreamerInfo.cxx:84
Add< Float_t >
template TMatrixFSparse & Add< Float_t >(TMatrixFSparse &target, Float_t scalar, const TMatrixFSparse &source)
TMatrixTSparse::RowNorm
virtual Element RowNorm() const
Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
Definition: TMatrixTSparse.cxx:2093
ROOT::Math::beta
double beta(double x, double y)
Calculates the beta function.
Definition: SpecFuncMathCore.cxx:111
ElementMult< Float_t >
template TMatrixFSparse & ElementMult< Float_t >(TMatrixFSparse &target, const TMatrixFSparse &source)
TMatrixTBase::IsValid
Bool_t IsValid() const
Definition: TMatrixTBase.h:157
TMatrixTBase::GetTol
Element GetTol() const
Definition: TMatrixTBase.h:139
AreCompatible< Float_t >
template Bool_t AreCompatible< Float_t >(const TMatrixFSparse &m1, const TMatrixFSparse &m2, Int_t verbose)
ROOT::Math::Blas::AMultB
void AMultB(int n, int m, int k, const double *A, const double *B, double *C)
Definition: cblas.cxx:11
Info
void Info(const char *location, const char *msgfmt,...)
AreCompatible
Bool_t AreCompatible(const TMatrixTSparse< Element > &m1, const TMatrixTSparse< Element > &m2, Int_t verbose)
Definition: TMatrixTSparse.cxx:2722
TMatrixTSparse::Randomize
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
randomize matrix element values
Definition: TMatrixTSparse.cxx:2349
TMatrixTSparse::AMultBt
void AMultBt(const TMatrixTSparse< Element > &a, const TMatrixTSparse< Element > &b, Int_t constr=0)
General matrix multiplication.
Definition: TMatrixTSparse.cxx:484
ElementMult
TMatrixTSparse< Element > & ElementMult(TMatrixTSparse< Element > &target, const TMatrixTSparse< Element > &source)
Multiply target by the source, element-by-element.
Definition: TMatrixTSparse.cxx:2677
TObject::operator=
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.cxx:102
TMatrixTSparse::RandomizePD
virtual TMatrixTSparse< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
Definition: TMatrixTSparse.cxx:2396
TMatrixTBase::fNrowIndex
Int_t fNrowIndex
Definition: TMatrixTBase.h:106
Add
TMatrixTSparse< Element > & Add(TMatrixTSparse< Element > &target, Element scalar, const TMatrixTSparse< Element > &source)
Modify addition: target += scalar * source.
Definition: TMatrixTSparse.cxx:2667
Error
void Error(const char *location, const char *msgfmt,...)
TMatrixTBase::GetNrows
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
TMatrixTSparse::operator-=
TMatrixTSparse< Element > & operator-=(Element val)
Subtract val from every element of the matrix.
Definition: TMatrixTSparse.cxx:2317
TMatrixTSparse::Zero
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
Definition: TMatrixTSparse.cxx:2035
r
ROOT::R::TRInterface & r
Definition: Object.C:4
TMatrixTSparse::GetRowIndexArray
virtual const Int_t * GetRowIndexArray() const
Definition: TMatrixTSparse.h:217
RootCsg::length
Double_t length(const TVector2 &v)
Definition: CsgOps.cxx:347
TMatrixTBase::DoubleLexSort
static void DoubleLexSort(Int_t n, Int_t *first, Int_t *second, Element *data)
default kTRUE, when Use array kFALSE
TMatrixTSparse::ColNorm
virtual Element ColNorm() const
Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
Definition: TMatrixTSparse.cxx:2122
ElementDiv< Double_t >
template TMatrixDSparse & ElementDiv< Double_t >(TMatrixDSparse &target, const TMatrixDSparse &source)
v
SVector< double, 2 > v
Definition: Dict.h:5
TMatrixTSparse::operator*=
TMatrixTSparse< Element > & operator*=(Element val)
Multiply every element of the matrix with val.
Definition: TMatrixTSparse.cxx:2333
UInt_t
unsigned int UInt_t
Definition: RtypesCore.h:42
m
TMarker * m
Definition: textangle.C:8
verbose
bool verbose
Definition: binarySearchTime.cxx:27
TMatrixTSparse::SetSparseIndex
TMatrixTSparse< Element > & SetSparseIndex(Int_t nelem_new)
Increase/decrease the number of non-zero elements to nelems_new.
Definition: TMatrixTSparse.cxx:1255
operator+
TMatrixTSparse< Element > operator+(const TMatrixTSparse< Element > &source1, const TMatrixTSparse< Element > &source2)
Definition: TMatrixTSparse.cxx:2525
TMatrixT::GetMatrixArray
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
gMatrixCheck
R__EXTERN Int_t gMatrixCheck
Definition: TMatrixTBase.h:91
TMatrixTBase::GetColLwb
Int_t GetColLwb() const
Definition: TMatrixTBase.h:135
epsilon
REAL epsilon
Definition: triangle.c:617
TMatrixTSparse::Allocate
void Allocate(Int_t nrows, Int_t ncols, Int_t row_lwb=0, Int_t col_lwb=0, Int_t init=0, Int_t nr_nonzeros=0)
Allocate new matrix.
Definition: TMatrixTSparse.cxx:318
TBuffer::ReadClassBuffer
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
ElementDiv
TMatrixTSparse< Element > & ElementDiv(TMatrixTSparse< Element > &target, const TMatrixTSparse< Element > &source)
Divide target by the source, element-by-element.
Definition: TMatrixTSparse.cxx:2697
init
static Int_t init()
Definition: RooClassFactory.cxx:54
templateClassImp
templateClassImp(TMatrixTSparse) template< class Element > TMatrixTSparse< Element >
Space is allocated for row/column indices and data, but the sparse structure information has still to...
Definition: TMatrixTSparse.cxx:78
Double_t
double Double_t
Definition: RtypesCore.h:55
TMatrixTSparse::Use
TMatrixTSparse< Element > & Use(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Int_t nr_nonzeros, Int_t *pRowIndex, Int_t *pColIndex, Element *pData)
Definition: TMatrixTSparse.cxx:1704
TMatrixTSparse::operator()
Element operator()(Int_t rown, Int_t coln) const
Definition: TMatrixTSparse.cxx:2194
TMatrixTSparse::SetSparseIndexAB
TMatrixTSparse< Element > & SetSparseIndexAB(const TMatrixTSparse< Element > &a, const TMatrixTSparse< Element > &b)
Set the row/column indices to the "sum" of matrices a and b It is checked that enough space has been ...
Definition: TMatrixTSparse.cxx:1332
printf
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
TMatrixTBase::GetNcols
Int_t GetNcols() const
Definition: TMatrixTBase.h:137
AreCompatible< Double_t >
template Bool_t AreCompatible< Double_t >(const TMatrixDSparse &m1, const TMatrixDSparse &m2, Int_t verbose)
Cppyy::Allocate
TCppObject_t Allocate(TCppType_t type)
Definition: Cppyy.cxx:228
TMatrixTSparse::GetColIndexArray
virtual const Int_t * GetColIndexArray() const
Definition: TMatrixTSparse.h:219
TMatrixTSparse::ResizeTo
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1)
Set size of the matrix to nrows x ncols with nr_nonzeros non-zero entries if nr_nonzeros > 0 ...
Definition: TMatrixTSparse.cxx:1511
TMatrixTSparse::operator+=
TMatrixTSparse< Element > & operator+=(Element val)
Add val to every element of the matrix.
Definition: TMatrixTSparse.cxx:2301
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
AMultBt
void AMultBt(const Element *const ap, Int_t na, Int_t ncolsa, const Element *const bp, Int_t nb, Int_t ncolsb, Element *cp)
Elementary routine to calculate matrix multiplication A*B^T.
Definition: TMatrixT.cxx:3116
TMatrixTSparse::ExtractRow
virtual void ExtractRow(Int_t row, Int_t col, Element *v, Int_t n=-1) const
Store in array v, n matrix elements of row rown starting at column coln.
Definition: TMatrixTSparse.cxx:443
TMatrixTSparse::GetSub
virtual TMatrixTBase< Element > & GetSub(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, TMatrixTBase< Element > &target, Option_t *option="S") const
Get submatrix [row_lwb..row_upb][col_lwb..col_upb]; The indexing range of the returned matrix depends...
Definition: TMatrixTSparse.cxx:1744
TMatrixTSparse::operator=
TMatrixTSparse< Element > & operator=(const TMatrixT< Element > &source)
Notice that the sparsity of the matrix is NOT changed : its fRowIndex/fColIndex are used ! ...
Definition: TMatrixTSparse.cxx:2249
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
TMath::LocMin
Long64_t LocMin(Long64_t n, const T *a)
Definition: TMath.h:695
kTRUE
const Bool_t kTRUE
Definition: Rtypes.h:91
operator-
TMatrixTSparse< Element > operator-(const TMatrixTSparse< Element > &source1, const TMatrixTSparse< Element > &source2)
Definition: TMatrixTSparse.cxx:2572
TMatrixTSparse::APlusB
void APlusB(const TMatrixTSparse< Element > &a, const TMatrixTSparse< Element > &b, Int_t constr=0)
General matrix addition.
Definition: TMatrixTSparse.cxx:775
norm
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
Drand
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
Definition: TMatrixTUtils.cxx:1852
operator*
TMatrixTSparse< Element > operator*(const TMatrixTSparse< Element > &source1, const TMatrixTSparse< Element > &source2)
Definition: TMatrixTSparse.cxx:2619
n
const Int_t n
Definition: legend1.C:16
TMath::BinarySearch
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMath.h:944
TMatrixTSparse::SetSub
virtual TMatrixTBase< Element > & SetSub(Int_t row_lwb, Int_t col_lwb, const TMatrixTBase< Element > &source)
Insert matrix source starting at [row_lwb][col_lwb], thereby overwriting the part [row_lwb...
Definition: TMatrixTSparse.cxx:1844
TMatrixTBase::GetColIndexArray
virtual const Int_t * GetColIndexArray() const =0
TBuffer::ReadVersion
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
TMatrixTBase::GetRowIndexArray
virtual const Int_t * GetRowIndexArray() const =0