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