Logo ROOT  

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