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