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