Logo ROOT  

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