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