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