#include "TMatrixTSparse.h"
#include "TMatrixT.h"
#include "TMath.h"
#ifndef R__ALPHA
templateClassImp(TMatrixTSparse)
#endif
template<class Element>
TMatrixTSparse<Element>::TMatrixTSparse(Int_t no_rows,Int_t no_cols)
{
Allocate(no_rows,no_cols,0,0,1);
}
template<class Element>
TMatrixTSparse<Element>::TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb)
{
Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb,1);
}
template<class Element>
TMatrixTSparse<Element>::TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
Int_t nr,Int_t *row, Int_t *col,Element *data)
{
const Int_t irowmin = TMath::LocMin(nr,row);
const Int_t irowmax = TMath::LocMax(nr,row);
const Int_t icolmin = TMath::LocMin(nr,col);
const Int_t icolmax = TMath::LocMax(nr,col);
if (row[irowmin] < row_lwb || row[irowmax] > row_upb) {
Error("TMatrixTSparse","Inconsistency between row index and its range");
if (row[irowmin] < row_lwb) {
Info("TMatrixTSparse","row index lower bound adjusted to %d",row[irowmin]);
row_lwb = row[irowmin];
}
if (row[irowmax] > row_upb) {
Info("TMatrixTSparse","row index upper bound adjusted to %d",row[irowmax]);
col_lwb = col[irowmax];
}
}
if (col[icolmin] < col_lwb || col[icolmax] > col_upb) {
Error("TMatrixTSparse","Inconsistency between column index and its range");
if (col[icolmin] < col_lwb) {
Info("TMatrixTSparse","column index lower bound adjusted to %d",col[icolmin]);
col_lwb = col[icolmin];
}
if (col[icolmax] > col_upb) {
Info("TMatrixTSparse","column index upper bound adjusted to %d",col[icolmax]);
col_upb = col[icolmax];
}
}
Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb,1,nr);
SetMatrixArray(nr,row,col,data);
}
template<class Element>
TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixTSparse<Element> &another) : TMatrixTBase<Element>(another)
{
Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb(),1,
another.GetNoElements());
memcpy(fRowIndex,another.GetRowIndexArray(),this->fNrowIndex*sizeof(Int_t));
memcpy(fColIndex,another.GetColIndexArray(),this->fNelems*sizeof(Int_t));
*this = another;
}
template<class Element>
TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixT<Element> &another) : TMatrixTBase<Element>(another)
{
const Int_t nr_nonzeros = another.NonZeros();
Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb(),1,nr_nonzeros);
SetSparseIndex(another);
*this = another;
}
template<class Element>
TMatrixTSparse<Element>::TMatrixTSparse(EMatrixCreatorsOp1 op,const TMatrixTSparse<Element> &prototype)
{
R__ASSERT(prototype.IsValid());
Int_t nr_nonzeros = 0;
switch(op) {
case kZero:
{
Allocate(prototype.GetNrows(),prototype.GetNcols(),
prototype.GetRowLwb(),prototype.GetColLwb(),1,nr_nonzeros);
break;
}
case kUnit:
{
const Int_t nrows = prototype.GetNrows();
const Int_t ncols = prototype.GetNcols();
const Int_t rowLwb = prototype.GetRowLwb();
const Int_t colLwb = prototype.GetColLwb();
for (Int_t i = rowLwb; i <= rowLwb+nrows-1; i++)
for (Int_t j = colLwb; j <= colLwb+ncols-1; j++)
if (i==j) nr_nonzeros++;
Allocate(nrows,ncols,rowLwb,colLwb,1,nr_nonzeros);
UnitMatrix();
break;
}
case kTransposed:
{
Allocate(prototype.GetNcols(), prototype.GetNrows(),
prototype.GetColLwb(),prototype.GetRowLwb(),1,prototype.GetNoElements());
Transpose(prototype);
break;
}
case kAtA:
{
const TMatrixTSparse<Element> at(TMatrixTSparse<Element>::kTransposed,prototype);
AMultBt(at,at,1);
break;
}
default:
Error("TMatrixTSparse(EMatrixCreatorOp1)","operation %d not yet implemented", op);
}
}
template<class Element>
TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b)
{
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
switch(op) {
case kMult:
AMultB(a,b,1);
break;
case kMultTranspose:
AMultBt(a,b,1);
break;
case kPlus:
APlusB(a,b,1);
break;
case kMinus:
AMinusB(a,b,1);
break;
default:
Error("TMatrixTSparse(EMatrixCreatorOp2)", "operation %d not yet implemented",op);
}
}
template<class Element>
TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixT<Element> &b)
{
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
switch(op) {
case kMult:
AMultB(a,b,1);
break;
case kMultTranspose:
AMultBt(a,b,1);
break;
case kPlus:
APlusB(a,b,1);
break;
case kMinus:
AMinusB(a,b,1);
break;
default:
Error("TMatrixTSparse(EMatrixCreatorOp2)", "operation %d not yet implemented",op);
}
}
template<class Element>
TMatrixTSparse<Element>::TMatrixTSparse(const TMatrixT<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b)
{
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
switch(op) {
case kMult:
AMultB(a,b,1);
break;
case kMultTranspose:
AMultBt(a,b,1);
break;
case kPlus:
APlusB(a,b,1);
break;
case kMinus:
AMinusB(a,b,1);
break;
default:
Error("TMatrixTSparse(EMatrixCreatorOp2)", "operation %d not yet implemented",op);
}
}
template<class Element>
void TMatrixTSparse<Element>::Allocate(Int_t no_rows,Int_t no_cols,Int_t row_lwb,Int_t col_lwb,
Int_t init,Int_t nr_nonzeros)
{
if ( (nr_nonzeros > 0 && (no_rows == 0 || no_cols == 0)) ||
(no_rows < 0 || no_cols < 0 || nr_nonzeros < 0) )
{
Error("Allocate","no_rows=%d no_cols=%d non_zeros=%d",no_rows,no_cols,nr_nonzeros);
this->Invalidate();
return;
}
this->MakeValid();
this->fNrows = no_rows;
this->fNcols = no_cols;
this->fRowLwb = row_lwb;
this->fColLwb = col_lwb;
this->fNrowIndex = this->fNrows+1;
this->fNelems = nr_nonzeros;
this->fIsOwner = kTRUE;
this->fTol = std::numeric_limits<Element>::epsilon();
fRowIndex = new Int_t[this->fNrowIndex];
if (init)
memset(fRowIndex,0,this->fNrowIndex*sizeof(Int_t));
if (this->fNelems > 0) {
fElements = new Element[this->fNelems];
fColIndex = new Int_t [this->fNelems];
if (init) {
memset(fElements,0,this->fNelems*sizeof(Element));
memset(fColIndex,0,this->fNelems*sizeof(Int_t));
}
} else {
fElements = 0;
fColIndex = 0;
}
}
template<class Element>
TMatrixTBase<Element> &TMatrixTSparse<Element>::InsertRow(Int_t rown,Int_t coln,const Element *v,Int_t n)
{
const Int_t arown = rown-this->fRowLwb;
const Int_t acoln = coln-this->fColLwb;
const Int_t nr = (n > 0) ? n : this->fNcols;
if (gMatrixCheck) {
if (arown >= this->fNrows || arown < 0) {
Error("InsertRow","row %d out of matrix range",rown);
return *this;
}
if (acoln >= this->fNcols || acoln < 0) {
Error("InsertRow","column %d out of matrix range",coln);
return *this;
}
if (acoln+nr > this->fNcols || nr < 0) {
Error("InsertRow","row length %d out of range",nr);
return *this;
}
}
const Int_t sIndex = fRowIndex[arown];
const Int_t eIndex = fRowIndex[arown+1];
Int_t nslots = 0;
Int_t lIndex = sIndex-1;
Int_t rIndex = sIndex-1;
Int_t index;
for (index = sIndex; index < eIndex; index++) {
const Int_t icol = fColIndex[index];
rIndex++;
if (icol >= acoln+nr) break;
if (icol >= acoln) nslots++;
else lIndex++;
}
const Int_t nelems_old = this->fNelems;
const Int_t *colIndex_old = fColIndex;
const Element *elements_old = fElements;
const Int_t ndiff = nr-nslots;
this->fNelems += ndiff;
fColIndex = new Int_t[this->fNelems];
fElements = new Element[this->fNelems];
for (Int_t irow = arown+1; irow < this->fNrows+1; irow++)
fRowIndex[irow] += ndiff;
if (lIndex+1 > 0) {
memmove(fColIndex,colIndex_old,(lIndex+1)*sizeof(Int_t));
memmove(fElements,elements_old,(lIndex+1)*sizeof(Element));
}
if (nelems_old > 0 && nelems_old-rIndex > 0) {
memmove(fColIndex+rIndex+ndiff,colIndex_old+rIndex,(nelems_old-rIndex)*sizeof(Int_t));
memmove(fElements+rIndex+ndiff,elements_old+rIndex,(nelems_old-rIndex)*sizeof(Element));
}
index = lIndex+1;
for (Int_t i = 0; i < nr; i++) {
fColIndex[index] = acoln+i;
fElements[index] = v[i];
index++;
}
if (colIndex_old) delete [] (Int_t*) colIndex_old;
if (elements_old) delete [] (Element*) elements_old;
R__ASSERT(this->fNelems == fRowIndex[this->fNrowIndex-1]);
return *this;
}
template<class Element>
void TMatrixTSparse<Element>::ExtractRow(Int_t rown, Int_t coln, Element *v,Int_t n) const
{
const Int_t arown = rown-this->fRowLwb;
const Int_t acoln = coln-this->fColLwb;
const Int_t nr = (n > 0) ? n : this->fNcols;
if (gMatrixCheck) {
if (arown >= this->fNrows || arown < 0) {
Error("ExtractRow","row %d out of matrix range",rown);
return;
}
if (acoln >= this->fNcols || acoln < 0) {
Error("ExtractRow","column %d out of matrix range",coln);
return;
}
if (acoln+nr > this->fNcols || nr < 0) {
Error("ExtractRow","row length %d out of range",nr);
return;
}
}
const Int_t sIndex = fRowIndex[arown];
const Int_t eIndex = fRowIndex[arown+1];
memset(v,0,nr*sizeof(Element));
const Int_t * const pColIndex = GetColIndexArray();
const Element * const pData = GetMatrixArray();
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = pColIndex[index];
if (icol < acoln || icol >= acoln+nr) continue;
v[icol-acoln] = pData[index];
}
}
template<class Element>
void TMatrixTSparse<Element>::AMultBt(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
Error("AMultBt","A and B columns incompatible");
return;
}
if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
Error("AMultB","this = &a");
return;
}
if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
Error("AMultB","this = &b");
return;
}
}
const Int_t * const pRowIndexa = a.GetRowIndexArray();
const Int_t * const pColIndexa = a.GetColIndexArray();
const Int_t * const pRowIndexb = b.GetRowIndexArray();
const Int_t * const pColIndexb = b.GetColIndexArray();
Int_t *pRowIndexc;
Int_t *pColIndexc;
if (constr) {
Int_t nr_nonzero_rowa = 0;
{
for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++)
if (pRowIndexa[irowa] < pRowIndexa[irowa+1])
nr_nonzero_rowa++;
}
Int_t nr_nonzero_rowb = 0;
{
for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++)
if (pRowIndexb[irowb] < pRowIndexb[irowb+1])
nr_nonzero_rowb++;
}
const Int_t nc = nr_nonzero_rowa*nr_nonzero_rowb;
Allocate(a.GetNrows(),b.GetNrows(),a.GetRowLwb(),b.GetRowLwb(),1,nc);
pRowIndexc = this->GetRowIndexArray();
pColIndexc = this->GetColIndexArray();
pRowIndexc[0] = 0;
Int_t ielem = 0;
for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++) {
pRowIndexc[irowa+1] = pRowIndexc[irowa];
if (pRowIndexa[irowa] >= pRowIndexa[irowa+1]) continue;
for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++) {
if (pRowIndexb[irowb] >= pRowIndexb[irowb+1]) continue;
pRowIndexc[irowa+1]++;
pColIndexc[ielem++] = irowb;
}
}
} else {
pRowIndexc = this->GetRowIndexArray();
pColIndexc = this->GetColIndexArray();
}
const Element * const pDataa = a.GetMatrixArray();
const Element * const pDatab = b.GetMatrixArray();
Element * const pDatac = this->GetMatrixArray();
Int_t indexc_r = 0;
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t sIndexa = pRowIndexa[irowc];
const Int_t eIndexa = pRowIndexa[irowc+1];
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
const Int_t sIndexb = pRowIndexb[icolc];
const Int_t eIndexb = pRowIndexb[icolc+1];
Element sum = 0.0;
Int_t indexb = sIndexb;
for (Int_t indexa = sIndexa; indexa < eIndexa && indexb < eIndexb; indexa++) {
const Int_t icola = pColIndexa[indexa];
while (indexb < eIndexb && pColIndexb[indexb] <= icola) {
if (icola == pColIndexb[indexb]) {
sum += pDataa[indexa]*pDatab[indexb];
break;
}
indexb++;
}
}
if (sum != 0.0) {
pColIndexc[indexc_r] = icolc;
pDatac[indexc_r] = sum;
indexc_r++;
}
}
pRowIndexc[irowc+1] = indexc_r;
}
if (constr)
SetSparseIndex(indexc_r);
}
template<class Element>
void TMatrixTSparse<Element>::AMultBt(const TMatrixTSparse<Element> &a,const TMatrixT<Element> &b,Int_t constr)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
Error("AMultBt","A and B columns incompatible");
return;
}
if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
Error("AMultB","this = &a");
return;
}
if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
Error("AMultB","this = &b");
return;
}
}
const Int_t * const pRowIndexa = a.GetRowIndexArray();
const Int_t * const pColIndexa = a.GetColIndexArray();
Int_t *pRowIndexc;
Int_t *pColIndexc;
if (constr) {
Int_t nr_nonzero_rowa = 0;
{
for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++)
if (pRowIndexa[irowa] < pRowIndexa[irowa+1])
nr_nonzero_rowa++;
}
Int_t nr_nonzero_rowb = b.GetNrows();
const Int_t nc = nr_nonzero_rowa*nr_nonzero_rowb;
Allocate(a.GetNrows(),b.GetNrows(),a.GetRowLwb(),b.GetRowLwb(),1,nc);
pRowIndexc = this->GetRowIndexArray();
pColIndexc = this->GetColIndexArray();
pRowIndexc[0] = 0;
Int_t ielem = 0;
for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++) {
pRowIndexc[irowa+1] = pRowIndexc[irowa];
for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++) {
pRowIndexc[irowa+1]++;
pColIndexc[ielem++] = irowb;
}
}
} else {
pRowIndexc = this->GetRowIndexArray();
pColIndexc = this->GetColIndexArray();
}
const Element * const pDataa = a.GetMatrixArray();
const Element * const pDatab = b.GetMatrixArray();
Element * const pDatac = this->GetMatrixArray();
Int_t indexc_r = 0;
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t sIndexa = pRowIndexa[irowc];
const Int_t eIndexa = pRowIndexa[irowc+1];
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
const Int_t off = icolc*b.GetNcols();
Element sum = 0.0;
for (Int_t indexa = sIndexa; indexa < eIndexa; indexa++) {
const Int_t icola = pColIndexa[indexa];
sum += pDataa[indexa]*pDatab[off+icola];
}
if (sum != 0.0) {
pColIndexc[indexc_r] = icolc;
pDatac[indexc_r] = sum;
indexc_r++;
}
}
pRowIndexc[irowc+1]= indexc_r;
}
if (constr)
SetSparseIndex(indexc_r);
}
template<class Element>
void TMatrixTSparse<Element>::AMultBt(const TMatrixT<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
Error("AMultBt","A and B columns incompatible");
return;
}
if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
Error("AMultB","this = &a");
return;
}
if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
Error("AMultB","this = &b");
return;
}
}
const Int_t * const pRowIndexb = b.GetRowIndexArray();
const Int_t * const pColIndexb = b.GetColIndexArray();
Int_t *pRowIndexc;
Int_t *pColIndexc;
if (constr) {
Int_t nr_nonzero_rowa = a.GetNrows();
Int_t nr_nonzero_rowb = 0;
{
for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++)
if (pRowIndexb[irowb] < pRowIndexb[irowb+1])
nr_nonzero_rowb++;
}
const Int_t nc = nr_nonzero_rowa*nr_nonzero_rowb;
Allocate(a.GetNrows(),b.GetNrows(),a.GetRowLwb(),b.GetRowLwb(),1,nc);
pRowIndexc = this->GetRowIndexArray();
pColIndexc = this->GetColIndexArray();
pRowIndexc[0] = 0;
Int_t ielem = 0;
for (Int_t irowa = 0; irowa < a.GetNrows(); irowa++) {
pRowIndexc[irowa+1] = pRowIndexc[irowa];
for (Int_t irowb = 0; irowb < b.GetNrows(); irowb++) {
if (pRowIndexb[irowb] >= pRowIndexb[irowb+1]) continue;
pRowIndexc[irowa+1]++;
pColIndexc[ielem++] = irowb;
}
}
} else {
pRowIndexc = this->GetRowIndexArray();
pColIndexc = this->GetColIndexArray();
}
const Element * const pDataa = a.GetMatrixArray();
const Element * const pDatab = b.GetMatrixArray();
Element * const pDatac = this->GetMatrixArray();
Int_t indexc_r = 0;
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t off = irowc*a.GetNcols();
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
const Int_t sIndexb = pRowIndexb[icolc];
const Int_t eIndexb = pRowIndexb[icolc+1];
Element sum = 0.0;
for (Int_t indexb = sIndexb; indexb < eIndexb; indexb++) {
const Int_t icolb = pColIndexb[indexb];
sum += pDataa[off+icolb]*pDatab[indexb];
}
if (sum != 0.0) {
pColIndexc[indexc_r] = icolc;
pDatac[indexc_r] = sum;
indexc_r++;
}
}
pRowIndexc[irowc+1] = indexc_r;
}
if (constr)
SetSparseIndex(indexc_r);
}
template<class Element>
void TMatrixTSparse<Element>::APlusB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
Error("APlusB(const TMatrixTSparse &,const TMatrixTSparse &","matrices not compatible");
return;
}
if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
Error("APlusB","this = &a");
return;
}
if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
Error("APlusB","this = &b");
return;
}
}
const Int_t * const pRowIndexa = a.GetRowIndexArray();
const Int_t * const pRowIndexb = b.GetRowIndexArray();
const Int_t * const pColIndexa = a.GetColIndexArray();
const Int_t * const pColIndexb = b.GetColIndexArray();
if (constr) {
Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
SetSparseIndexAB(a,b);
}
Int_t * const pRowIndexc = this->GetRowIndexArray();
Int_t * const pColIndexc = this->GetColIndexArray();
const Element * const pDataa = a.GetMatrixArray();
const Element * const pDatab = b.GetMatrixArray();
Element * const pDatac = this->GetMatrixArray();
Int_t indexc_r = 0;
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t sIndexa = pRowIndexa[irowc];
const Int_t eIndexa = pRowIndexa[irowc+1];
const Int_t sIndexb = pRowIndexb[irowc];
const Int_t eIndexb = pRowIndexb[irowc+1];
Int_t indexa = sIndexa;
Int_t indexb = sIndexb;
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
Element sum = 0.0;
while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
if (icolc == pColIndexa[indexa]) {
sum += pDataa[indexa];
break;
}
indexa++;
}
while (indexb < eIndexb && pColIndexb[indexb] <= icolc) {
if (icolc == pColIndexb[indexb]) {
sum += pDatab[indexb];
break;
}
indexb++;
}
if (sum != 0.0) {
pColIndexc[indexc_r] = icolc;
pDatac[indexc_r] = sum;
indexc_r++;
}
}
pRowIndexc[irowc+1] = indexc_r;
}
if (constr)
SetSparseIndex(indexc_r);
}
template<class Element>
void TMatrixTSparse<Element>::APlusB(const TMatrixTSparse<Element> &a,const TMatrixT<Element> &b,Int_t constr)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
Error("APlusB(const TMatrixTSparse &,const TMatrixT &","matrices not compatible");
return;
}
if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
Error("APlusB","this = &a");
return;
}
if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
Error("APlusB","this = &b");
return;
}
}
if (constr) {
Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
SetSparseIndexAB(a,b);
}
Int_t * const pRowIndexc = this->GetRowIndexArray();
Int_t * const pColIndexc = this->GetColIndexArray();
const Int_t * const pRowIndexa = a.GetRowIndexArray();
const Int_t * const pColIndexa = a.GetColIndexArray();
const Element * const pDataa = a.GetMatrixArray();
const Element * const pDatab = b.GetMatrixArray();
Element * const pDatac = this->GetMatrixArray();
Int_t indexc_r = 0;
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t sIndexa = pRowIndexa[irowc];
const Int_t eIndexa = pRowIndexa[irowc+1];
const Int_t off = irowc*this->GetNcols();
Int_t indexa = sIndexa;
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
Element sum = pDatab[off+icolc];
while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
if (icolc == pColIndexa[indexa]) {
sum += pDataa[indexa];
break;
}
indexa++;
}
if (sum != 0.0) {
pColIndexc[indexc_r] = icolc;
pDatac[indexc_r] = sum;
indexc_r++;
}
}
pRowIndexc[irowc+1] = indexc_r;
}
if (constr)
SetSparseIndex(indexc_r);
}
template<class Element>
void TMatrixTSparse<Element>::AMinusB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
Error("AMinusB(const TMatrixTSparse &,const TMatrixTSparse &","matrices not compatible");
return;
}
if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
Error("AMinusB","this = &a");
return;
}
if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
Error("AMinusB","this = &b");
return;
}
}
const Int_t * const pRowIndexa = a.GetRowIndexArray();
const Int_t * const pRowIndexb = b.GetRowIndexArray();
const Int_t * const pColIndexa = a.GetColIndexArray();
const Int_t * const pColIndexb = b.GetColIndexArray();
if (constr) {
Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
SetSparseIndexAB(a,b);
}
Int_t * const pRowIndexc = this->GetRowIndexArray();
Int_t * const pColIndexc = this->GetColIndexArray();
const Element * const pDataa = a.GetMatrixArray();
const Element * const pDatab = b.GetMatrixArray();
Element * const pDatac = this->GetMatrixArray();
Int_t indexc_r = 0;
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t sIndexa = pRowIndexa[irowc];
const Int_t eIndexa = pRowIndexa[irowc+1];
const Int_t sIndexb = pRowIndexb[irowc];
const Int_t eIndexb = pRowIndexb[irowc+1];
Int_t indexa = sIndexa;
Int_t indexb = sIndexb;
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
Element sum = 0.0;
while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
if (icolc == pColIndexa[indexa]) {
sum += pDataa[indexa];
break;
}
indexa++;
}
while (indexb < eIndexb && pColIndexb[indexb] <= icolc) {
if (icolc == pColIndexb[indexb]) {
sum -= pDatab[indexb];
break;
}
indexb++;
}
if (sum != 0.0) {
pColIndexc[indexc_r] = icolc;
pDatac[indexc_r] = sum;
indexc_r++;
}
}
pRowIndexc[irowc+1] = indexc_r;
}
if (constr)
SetSparseIndex(indexc_r);
}
template<class Element>
void TMatrixTSparse<Element>::AMinusB(const TMatrixTSparse<Element> &a,const TMatrixT<Element> &b,Int_t constr)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
Error("AMinusB(const TMatrixTSparse &,const TMatrixT &","matrices not compatible");
return;
}
if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
Error("AMinusB","this = &a");
return;
}
if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
Error("AMinusB","this = &b");
return;
}
}
if (constr) {
Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
SetSparseIndexAB(a,b);
}
Int_t * const pRowIndexc = this->GetRowIndexArray();
Int_t * const pColIndexc = this->GetColIndexArray();
const Int_t * const pRowIndexa = a.GetRowIndexArray();
const Int_t * const pColIndexa = a.GetColIndexArray();
const Element * const pDataa = a.GetMatrixArray();
const Element * const pDatab = b.GetMatrixArray();
Element * const pDatac = this->GetMatrixArray();
Int_t indexc_r = 0;
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t sIndexa = pRowIndexa[irowc];
const Int_t eIndexa = pRowIndexa[irowc+1];
const Int_t off = irowc*this->GetNcols();
Int_t indexa = sIndexa;
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
Element sum = -pDatab[off+icolc];
while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
if (icolc == pColIndexa[indexa]) {
sum += pDataa[indexa];
break;
}
indexa++;
}
if (sum != 0.0) {
pColIndexc[indexc_r] = icolc;
pDatac[indexc_r] = sum;
indexc_r++;
}
}
pRowIndexc[irowc+1] = indexc_r;
}
if (constr)
SetSparseIndex(indexc_r);
}
template<class Element>
void TMatrixTSparse<Element>::AMinusB(const TMatrixT<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
Error("AMinusB(const TMatrixT &,const TMatrixTSparse &","matrices not compatible");
return;
}
if (!constr && this->GetMatrixArray() == a.GetMatrixArray()) {
Error("AMinusB","this = &a");
return;
}
if (!constr && this->GetMatrixArray() == b.GetMatrixArray()) {
Error("AMinusB","this = &b");
return;
}
}
if (constr) {
Allocate(a.GetNrows(),a.GetNcols(),a.GetRowLwb(),a.GetColLwb());
SetSparseIndexAB(a,b);
}
Int_t * const pRowIndexc = this->GetRowIndexArray();
Int_t * const pColIndexc = this->GetColIndexArray();
const Int_t * const pRowIndexb = b.GetRowIndexArray();
const Int_t * const pColIndexb = b.GetColIndexArray();
const Element * const pDataa = a.GetMatrixArray();
const Element * const pDatab = b.GetMatrixArray();
Element * const pDatac = this->GetMatrixArray();
Int_t indexc_r = 0;
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t sIndexb = pRowIndexb[irowc];
const Int_t eIndexb = pRowIndexb[irowc+1];
const Int_t off = irowc*this->GetNcols();
Int_t indexb = sIndexb;
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
Element sum = pDataa[off+icolc];
while (indexb < eIndexb && pColIndexb[indexb] <= icolc) {
if (icolc == pColIndexb[indexb]) {
sum -= pDatab[indexb];
break;
}
indexb++;
}
if (sum != 0.0) {
pColIndexc[indexc_r] = icolc;
pDatac[indexc_r] = sum;
indexc_r++;
}
}
pRowIndexc[irowc+1] = indexc_r;
}
if (constr)
SetSparseIndex(indexc_r);
}
template<class Element>
void TMatrixTSparse<Element>::GetMatrix2Array(Element *data,Option_t * ) const
{
R__ASSERT(this->IsValid());
const Element * const elem = GetMatrixArray();
memcpy(data,elem,this->fNelems*sizeof(Element));
}
template<class Element>
TMatrixTBase<Element> &TMatrixTSparse<Element>::SetMatrixArray(Int_t nr,Int_t *row,Int_t *col,Element *data)
{
R__ASSERT(this->IsValid());
if (nr <= 0) {
Error("SetMatrixArray(Int_t,Int_t*,Int_t*,Element*","nr <= 0");
return *this;
}
const Int_t irowmin = TMath::LocMin(nr,row);
const Int_t irowmax = TMath::LocMax(nr,row);
const Int_t icolmin = TMath::LocMin(nr,col);
const Int_t icolmax = TMath::LocMax(nr,col);
R__ASSERT(row[irowmin] >= this->fRowLwb && row[irowmax] <= this->fRowLwb+this->fNrows-1);
R__ASSERT(col[icolmin] >= this->fColLwb && col[icolmax] <= this->fColLwb+this->fNcols-1);
if (row[irowmin] < this->fRowLwb || row[irowmax] > this->fRowLwb+this->fNrows-1) {
Error("SetMatrixArray","Inconsistency between row index and its range");
if (row[irowmin] < this->fRowLwb) {
Info("SetMatrixArray","row index lower bound adjusted to %d",row[irowmin]);
this->fRowLwb = row[irowmin];
}
if (row[irowmax] > this->fRowLwb+this->fNrows-1) {
Info("SetMatrixArray","row index upper bound adjusted to %d",row[irowmax]);
this->fNrows = row[irowmax]-this->fRowLwb+1;
}
}
if (col[icolmin] < this->fColLwb || col[icolmax] > this->fColLwb+this->fNcols-1) {
Error("SetMatrixArray","Inconsistency between column index and its range");
if (col[icolmin] < this->fColLwb) {
Info("SetMatrixArray","column index lower bound adjusted to %d",col[icolmin]);
this->fColLwb = col[icolmin];
}
if (col[icolmax] > this->fColLwb+this->fNcols-1) {
Info("SetMatrixArray","column index upper bound adjusted to %d",col[icolmax]);
this->fNcols = col[icolmax]-this->fColLwb+1;
}
}
DoubleLexSort(nr,row,col,data);
Int_t nr_nonzeros = 0;
const Element *ep = data;
const Element * const fp = data+nr;
while (ep < fp)
if (*ep++ != 0.0) nr_nonzeros++;
if (nr_nonzeros != this->fNelems) {
if (fColIndex) { delete [] fColIndex; fColIndex = 0; }
if (fElements) { delete [] fElements; fElements = 0; }
this->fNelems = nr_nonzeros;
if (this->fNelems > 0) {
fColIndex = new Int_t[nr_nonzeros];
fElements = new Element[nr_nonzeros];
} else {
fColIndex = 0;
fElements = 0;
}
}
if (this->fNelems <= 0)
return *this;
fRowIndex[0] = 0;
Int_t ielem = 0;
nr_nonzeros = 0;
for (Int_t irow = 1; irow < this->fNrows+1; irow++) {
if (ielem < nr && row[ielem] < irow) {
while (ielem < nr) {
if (data[ielem] != 0.0) {
fColIndex[nr_nonzeros] = col[ielem]-this->fColLwb;
fElements[nr_nonzeros] = data[ielem];
nr_nonzeros++;
}
ielem++;
if (ielem >= nr || row[ielem] != row[ielem-1])
break;
}
}
fRowIndex[irow] = nr_nonzeros;
}
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::SetSparseIndex(Int_t nelems_new)
{
if (nelems_new != this->fNelems) {
Int_t nr = TMath::Min(nelems_new,this->fNelems);
Int_t *oIp = fColIndex;
fColIndex = new Int_t[nelems_new];
memmove(fColIndex,oIp,nr*sizeof(Int_t));
if (oIp) delete [] oIp;
Element *oDp = fElements;
fElements = new Element[nelems_new];
memmove(fElements,oDp,nr*sizeof(Element));
if (oDp) delete [] oDp;
this->fNelems = nelems_new;
if (nelems_new > nr) {
memset(fElements+nr,0,(nelems_new-nr)*sizeof(Element));
memset(fColIndex+nr,0,(nelems_new-nr)*sizeof(Int_t));
} else {
for (Int_t irow = 0; irow < this->fNrowIndex; irow++)
if (fRowIndex[irow] > nelems_new)
fRowIndex[irow] = nelems_new;
}
}
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::SetSparseIndex(const TMatrixTBase<Element> &source)
{
if (gMatrixCheck) {
R__ASSERT(source.IsValid());
if (this->GetNrows() != source.GetNrows() || this->GetNcols() != source.GetNcols() ||
this->GetRowLwb() != source.GetRowLwb() || this->GetColLwb() != source.GetColLwb()) {
Error("SetSparseIndex","matrices not compatible");
return *this;
}
}
const Int_t nr_nonzeros = source.NonZeros();
if (nr_nonzeros != this->fNelems)
SetSparseIndex(nr_nonzeros);
if (source.GetRowIndexArray() && source.GetColIndexArray()) {
memmove(fRowIndex,source.GetRowIndexArray(),this->fNrowIndex*sizeof(Int_t));
memmove(fColIndex,source.GetColIndexArray(),this->fNelems*sizeof(Int_t));
} else {
const Element *ep = source.GetMatrixArray();
Int_t nr = 0;
for (Int_t irow = 0; irow < this->fNrows; irow++) {
fRowIndex[irow] = nr;
for (Int_t icol = 0; icol < this->fNcols; icol++) {
if (*ep != 0.0) {
fColIndex[nr] = icol;
nr++;
}
ep++;
}
}
fRowIndex[this->fNrows] = nr;
}
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
Error("SetSparseIndexAB","source matrices not compatible");
return *this;
}
if (this->GetNrows() != a.GetNrows() || this->GetNcols() != a.GetNcols() ||
this->GetRowLwb() != a.GetRowLwb() || this->GetColLwb() != a.GetColLwb()) {
Error("SetSparseIndexAB","matrix not compatible with source matrices");
return *this;
}
}
const Int_t * const pRowIndexa = a.GetRowIndexArray();
const Int_t * const pRowIndexb = b.GetRowIndexArray();
const Int_t * const pColIndexa = a.GetColIndexArray();
const Int_t * const pColIndexb = b.GetColIndexArray();
Int_t nc = 0;
for (Int_t irowc = 0; irowc < a.GetNrows(); irowc++) {
const Int_t sIndexa = pRowIndexa[irowc];
const Int_t eIndexa = pRowIndexa[irowc+1];
const Int_t sIndexb = pRowIndexb[irowc];
const Int_t eIndexb = pRowIndexb[irowc+1];
nc += eIndexa-sIndexa;
Int_t indexb = sIndexb;
for (Int_t indexa = sIndexa; indexa < eIndexa; indexa++) {
const Int_t icola = pColIndexa[indexa];
for (; indexb < eIndexb; indexb++) {
if (pColIndexb[indexb] >= icola) {
if (pColIndexb[indexb] == icola)
indexb++;
break;
}
nc++;
}
}
while (indexb < eIndexb) {
const Int_t icola = (eIndexa > 0) ? pColIndexa[eIndexa-1] : -1;
if (pColIndexb[indexb++] > icola)
nc++;
}
}
if (this->NonZeros() != nc)
SetSparseIndex(nc);
Int_t * const pRowIndexc = this->GetRowIndexArray();
Int_t * const pColIndexc = this->GetColIndexArray();
nc = 0;
pRowIndexc[0] = 0;
for (Int_t irowc = 0; irowc < a.GetNrows(); irowc++) {
const Int_t sIndexa = pRowIndexa[irowc];
const Int_t eIndexa = pRowIndexa[irowc+1];
const Int_t sIndexb = pRowIndexb[irowc];
const Int_t eIndexb = pRowIndexb[irowc+1];
Int_t indexb = sIndexb;
for (Int_t indexa = sIndexa; indexa < eIndexa; indexa++) {
const Int_t icola = pColIndexa[indexa];
for (; indexb < eIndexb; indexb++) {
if (pColIndexb[indexb] >= icola) {
if (pColIndexb[indexb] == icola)
indexb++;
break;
}
pColIndexc[nc++] = pColIndexb[indexb];
}
pColIndexc[nc++] = pColIndexa[indexa];
}
while (indexb < eIndexb) {
const Int_t icola = (eIndexa > 0) ? pColIndexa[eIndexa-1] : -1;
if (pColIndexb[indexb++] > icola)
pColIndexc[nc++] = pColIndexb[indexb-1];
}
pRowIndexc[irowc+1] = nc;
}
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::SetSparseIndexAB(const TMatrixT<Element> &a,const TMatrixTSparse<Element> &b)
{
if (gMatrixCheck) {
R__ASSERT(a.IsValid());
R__ASSERT(b.IsValid());
if (a.GetNrows() != b.GetNrows() || a.GetNcols() != b.GetNcols() ||
a.GetRowLwb() != b.GetRowLwb() || a.GetColLwb() != b.GetColLwb()) {
Error("SetSparseIndexAB","source matrices not compatible");
return *this;
}
if (this->GetNrows() != a.GetNrows() || this->GetNcols() != a.GetNcols() ||
this->GetRowLwb() != a.GetRowLwb() || this->GetColLwb() != a.GetColLwb()) {
Error("SetSparseIndexAB","matrix not compatible with source matrices");
return *this;
}
}
const Element * const pDataa = a.GetMatrixArray();
const Int_t * const pRowIndexb = b.GetRowIndexArray();
const Int_t * const pColIndexb = b.GetColIndexArray();
Int_t nc = a.NonZeros();
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t sIndexb = pRowIndexb[irowc];
const Int_t eIndexb = pRowIndexb[irowc+1];
const Int_t off = irowc*this->GetNcols();
Int_t indexb = sIndexb;
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
if (pDataa[off+icolc] != 0.0 || pColIndexb[indexb] > icolc) continue;
for (; indexb < eIndexb; indexb++) {
if (pColIndexb[indexb] >= icolc) {
if (pColIndexb[indexb] == icolc) {
nc++;
indexb++;
}
break;
}
}
}
}
if (this->NonZeros() != nc)
SetSparseIndex(nc);
Int_t * const pRowIndexc = this->GetRowIndexArray();
Int_t * const pColIndexc = this->GetColIndexArray();
nc = 0;
pRowIndexc[0] = 0;
for (Int_t irowc = 0; irowc < this->GetNrows(); irowc++) {
const Int_t sIndexb = pRowIndexb[irowc];
const Int_t eIndexb = pRowIndexb[irowc+1];
const Int_t off = irowc*this->GetNcols();
Int_t indexb = sIndexb;
for (Int_t icolc = 0; icolc < this->GetNcols(); icolc++) {
if (pDataa[off+icolc] != 0.0)
pColIndexc[nc++] = icolc;
else if (pColIndexb[indexb] <= icolc) {
for (; indexb < eIndexb; indexb++) {
if (pColIndexb[indexb] >= icolc) {
if (pColIndexb[indexb] == icolc)
pColIndexc[nc++] = pColIndexb[indexb++];
break;
}
}
}
}
pRowIndexc[irowc+1] = nc;
}
return *this;
}
template<class Element>
TMatrixTBase<Element> &TMatrixTSparse<Element>::ResizeTo(Int_t nrows,Int_t ncols,Int_t nr_nonzeros)
{
R__ASSERT(this->IsValid());
if (!this->fIsOwner) {
Error("ResizeTo(Int_t,Int_t,Int_t)","Not owner of data array,cannot resize");
return *this;
}
if (this->fNelems > 0) {
if (this->fNrows == nrows && this->fNcols == ncols &&
(this->fNelems == nr_nonzeros || nr_nonzeros < 0))
return *this;
else if (nrows == 0 || ncols == 0 || nr_nonzeros == 0) {
this->fNrows = nrows; this->fNcols = ncols;
Clear();
return *this;
}
const Element *elements_old = GetMatrixArray();
const Int_t *rowIndex_old = GetRowIndexArray();
const Int_t *colIndex_old = GetColIndexArray();
Int_t nrows_old = this->fNrows;
Int_t nrowIndex_old = this->fNrowIndex;
Int_t nelems_new;
if (nr_nonzeros > 0)
nelems_new = nr_nonzeros;
else {
nelems_new = 0;
for (Int_t irow = 0; irow < nrows_old; irow++) {
if (irow >= nrows) continue;
const Int_t sIndex = rowIndex_old[irow];
const Int_t eIndex = rowIndex_old[irow+1];
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = colIndex_old[index];
if (icol < ncols)
nelems_new++;
}
}
}
Allocate(nrows,ncols,0,0,1,nelems_new);
R__ASSERT(this->IsValid());
Element *elements_new = GetMatrixArray();
Int_t *rowIndex_new = GetRowIndexArray();
Int_t *colIndex_new = GetColIndexArray();
Int_t nelems_copy = 0;
rowIndex_new[0] = 0;
Bool_t cont = kTRUE;
for (Int_t irow = 0; irow < nrows_old && cont; irow++) {
if (irow >= nrows) continue;
const Int_t sIndex = rowIndex_old[irow];
const Int_t eIndex = rowIndex_old[irow+1];
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = colIndex_old[index];
if (icol < ncols) {
rowIndex_new[irow+1] = nelems_copy+1;
colIndex_new[nelems_copy] = icol;
elements_new[nelems_copy] = elements_old[index];
nelems_copy++;
}
if (nelems_copy >= nelems_new) {
cont = kFALSE;
break;
}
}
}
if (rowIndex_old) delete [] (Int_t*) rowIndex_old;
if (colIndex_old) delete [] (Int_t*) colIndex_old;
if (elements_old) delete [] (Element*) elements_old;
if (nrowIndex_old < this->fNrowIndex) {
for (Int_t irow = nrowIndex_old; irow < this->fNrowIndex; irow++)
rowIndex_new[irow] = rowIndex_new[nrowIndex_old-1];
}
} else {
const Int_t nelems_new = (nr_nonzeros >= 0) ? nr_nonzeros : 0;
Allocate(nrows,ncols,0,0,1,nelems_new);
}
return *this;
}
template<class Element>
TMatrixTBase<Element> &TMatrixTSparse<Element>::ResizeTo(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
Int_t nr_nonzeros)
{
R__ASSERT(this->IsValid());
if (!this->fIsOwner) {
Error("ResizeTo(Int_t,Int_t,Int_t,Int_t,Int_t)","Not owner of data array,cannot resize");
return *this;
}
const Int_t new_nrows = row_upb-row_lwb+1;
const Int_t new_ncols = col_upb-col_lwb+1;
if (this->fNelems > 0) {
if (this->fNrows == new_nrows && this->fNcols == new_ncols &&
this->fRowLwb == row_lwb && this->fColLwb == col_lwb &&
(this->fNelems == nr_nonzeros || nr_nonzeros < 0))
return *this;
else if (new_nrows == 0 || new_ncols == 0 || nr_nonzeros == 0) {
this->fNrows = new_nrows; this->fNcols = new_ncols;
this->fRowLwb = row_lwb; this->fColLwb = col_lwb;
Clear();
return *this;
}
const Int_t *rowIndex_old = GetRowIndexArray();
const Int_t *colIndex_old = GetColIndexArray();
const Element *elements_old = GetMatrixArray();
const Int_t nrowIndex_old = this->fNrowIndex;
const Int_t nrows_old = this->fNrows;
const Int_t rowLwb_old = this->fRowLwb;
const Int_t colLwb_old = this->fColLwb;
Int_t nelems_new;
if (nr_nonzeros > 0)
nelems_new = nr_nonzeros;
else {
nelems_new = 0;
for (Int_t irow = 0; irow < nrows_old; irow++) {
if (irow+rowLwb_old > row_upb || irow+rowLwb_old < row_lwb) continue;
const Int_t sIndex = rowIndex_old[irow];
const Int_t eIndex = rowIndex_old[irow+1];
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = colIndex_old[index];
if (icol+colLwb_old <= col_upb && icol+colLwb_old >= col_lwb)
nelems_new++;
}
}
}
Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1,nelems_new);
R__ASSERT(this->IsValid());
Int_t *rowIndex_new = GetRowIndexArray();
Int_t *colIndex_new = GetColIndexArray();
Element *elements_new = GetMatrixArray();
Int_t nelems_copy = 0;
rowIndex_new[0] = 0;
Bool_t cont = kTRUE;
const Int_t row_off = rowLwb_old-row_lwb;
const Int_t col_off = colLwb_old-col_lwb;
for (Int_t irow = 0; irow < nrows_old; irow++) {
if (irow+rowLwb_old > row_upb || irow+rowLwb_old < row_lwb) continue;
const Int_t sIndex = rowIndex_old[irow];
const Int_t eIndex = rowIndex_old[irow+1];
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = colIndex_old[index];
if (icol+colLwb_old <= col_upb && icol+colLwb_old >= col_lwb) {
rowIndex_new[irow+row_off+1] = nelems_copy+1;
colIndex_new[nelems_copy] = icol+col_off;
elements_new[nelems_copy] = elements_old[index];
nelems_copy++;
}
if (nelems_copy >= nelems_new) {
cont = kFALSE;
break;
}
}
}
if (rowIndex_old) delete [] (Int_t*) rowIndex_old;
if (colIndex_old) delete [] (Int_t*) colIndex_old;
if (elements_old) delete [] (Element*) elements_old;
if (nrowIndex_old < this->fNrowIndex) {
for (Int_t irow = nrowIndex_old; irow < this->fNrowIndex; irow++)
rowIndex_new[irow] = rowIndex_new[nrowIndex_old-1];
}
} else {
const Int_t nelems_new = (nr_nonzeros >= 0) ? nr_nonzeros : 0;
Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1,nelems_new);
}
return *this;
}
template<class Element>
TMatrixTSparse<Element> &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)
{
if (gMatrixCheck) {
if (row_upb < row_lwb) {
Error("Use","row_upb=%d < row_lwb=%d",row_upb,row_lwb);
return *this;
}
if (col_upb < col_lwb) {
Error("Use","col_upb=%d < col_lwb=%d",col_upb,col_lwb);
return *this;
}
}
Clear();
this->fNrows = row_upb-row_lwb+1;
this->fNcols = col_upb-col_lwb+1;
this->fRowLwb = row_lwb;
this->fColLwb = col_lwb;
this->fNrowIndex = this->fNrows+1;
this->fNelems = nr_nonzeros;
this->fIsOwner = kFALSE;
this->fTol = std::numeric_limits<Element>::epsilon();
fElements = pData;
fRowIndex = pRowIndex;
fColIndex = pColIndex;
return *this;
}
template<class Element>
TMatrixTBase<Element> &TMatrixTSparse<Element>::GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
TMatrixTBase<Element> &target,Option_t *option) const
{
if (gMatrixCheck) {
R__ASSERT(this->IsValid());
if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
Error("GetSub","row_lwb out-of-bounds");
return target;
}
if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
Error("GetSub","col_lwb out-of-bounds");
return target;
}
if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
Error("GetSub","row_upb out-of-bounds");
return target;
}
if (col_upb < this->fColLwb || col_upb > this->fColLwb+this->fNcols-1) {
Error("GetSub","col_upb out-of-bounds");
return target;
}
if (row_upb < row_lwb || col_upb < col_lwb) {
Error("GetSub","row_upb < row_lwb || col_upb < col_lwb");
return target;
}
}
TString opt(option);
opt.ToUpper();
const Int_t shift = (opt.Contains("S")) ? 1 : 0;
const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
Int_t nr_nonzeros = 0;
for (Int_t irow = 0; irow < this->fNrows; irow++) {
if (irow+this->fRowLwb > row_upb || irow+this->fRowLwb < row_lwb) continue;
const Int_t sIndex = fRowIndex[irow];
const Int_t eIndex = fRowIndex[irow+1];
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = fColIndex[index];
if (icol+this->fColLwb <= col_upb && icol+this->fColLwb >= col_lwb)
nr_nonzeros++;
}
}
target.ResizeTo(row_lwb_sub,row_upb_sub,col_lwb_sub,col_upb_sub,nr_nonzeros);
const Element *ep = this->GetMatrixArray();
Int_t *rowIndex_sub = target.GetRowIndexArray();
Int_t *colIndex_sub = target.GetColIndexArray();
Element *ep_sub = target.GetMatrixArray();
if (target.GetRowIndexArray() && target.GetColIndexArray()) {
Int_t nelems_copy = 0;
rowIndex_sub[0] = 0;
const Int_t row_off = this->fRowLwb-row_lwb;
const Int_t col_off = this->fColLwb-col_lwb;
for (Int_t irow = 0; irow < this->fNrows; irow++) {
if (irow+this->fRowLwb > row_upb || irow+this->fRowLwb < row_lwb) continue;
const Int_t sIndex = fRowIndex[irow];
const Int_t eIndex = fRowIndex[irow+1];
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = fColIndex[index];
if (icol+this->fColLwb <= col_upb && icol+this->fColLwb >= col_lwb) {
rowIndex_sub[irow+row_off+1] = nelems_copy+1;
colIndex_sub[nelems_copy] = icol+col_off;
ep_sub[nelems_copy] = ep[index];
nelems_copy++;
}
}
}
} else {
const Int_t row_off = this->fRowLwb-row_lwb;
const Int_t col_off = this->fColLwb-col_lwb;
const Int_t ncols_sub = col_upb_sub-col_lwb_sub+1;
for (Int_t irow = 0; irow < this->fNrows; irow++) {
if (irow+this->fRowLwb > row_upb || irow+this->fRowLwb < row_lwb) continue;
const Int_t sIndex = fRowIndex[irow];
const Int_t eIndex = fRowIndex[irow+1];
const Int_t off = (irow+row_off)*ncols_sub;
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = fColIndex[index];
if (icol+this->fColLwb <= col_upb && icol+this->fColLwb >= col_lwb)
ep_sub[off+icol+col_off] = ep[index];
}
}
}
return target;
}
template<class Element>
TMatrixTBase<Element> &TMatrixTSparse<Element>::SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source)
{
if (gMatrixCheck) {
R__ASSERT(this->IsValid());
R__ASSERT(source.IsValid());
if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
Error("SetSub","row_lwb out-of-bounds");
return *this;
}
if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
Error("SetSub","col_lwb out-of-bounds");
return *this;
}
if (row_lwb+source.GetNrows() > this->fRowLwb+this->fNrows || col_lwb+source.GetNcols() > this->fColLwb+this->fNcols) {
Error("SetSub","source matrix too large");
return *this;
}
}
const Int_t nRows_source = source.GetNrows();
const Int_t nCols_source = source.GetNcols();
Int_t nr_nonzeros = 0;
for (Int_t irow = 0; irow < this->fNrows; irow++) {
if (irow+this->fRowLwb >= row_lwb+nRows_source || irow+this->fRowLwb < row_lwb) continue;
const Int_t sIndex = fRowIndex[irow];
const Int_t eIndex = fRowIndex[irow+1];
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = fColIndex[index];
if (icol+this->fColLwb < col_lwb+nCols_source && icol+this->fColLwb >= col_lwb)
nr_nonzeros++;
}
}
const Int_t *rowIndex_s = source.GetRowIndexArray();
const Int_t *colIndex_s = source.GetColIndexArray();
const Element *elements_s = source.GetMatrixArray();
const Int_t nelems_old = this->fNelems;
const Int_t *rowIndex_old = GetRowIndexArray();
const Int_t *colIndex_old = GetColIndexArray();
const Element *elements_old = GetMatrixArray();
const Int_t nelems_new = nelems_old+source.NonZeros()-nr_nonzeros;
fRowIndex = new Int_t[this->fNrowIndex];
fColIndex = new Int_t[nelems_new];
fElements = new Element[nelems_new];
this->fNelems = nelems_new;
Int_t *rowIndex_new = GetRowIndexArray();
Int_t *colIndex_new = GetColIndexArray();
Element *elements_new = GetMatrixArray();
const Int_t row_off = row_lwb-this->fRowLwb;
const Int_t col_off = col_lwb-this->fColLwb;
Int_t nr = 0;
rowIndex_new[0] = 0;
for (Int_t irow = 0; irow < this->fNrows; irow++) {
rowIndex_new[irow+1] = rowIndex_new[irow];
Bool_t flagRow = kFALSE;
if (irow+this->fRowLwb < row_lwb+nRows_source && irow+this->fRowLwb >= row_lwb)
flagRow = kTRUE;
const Int_t sIndex_o = rowIndex_old[irow];
const Int_t eIndex_o = rowIndex_old[irow+1];
if (flagRow) {
const Int_t icol_left = col_off-1;
const Int_t left = TMath::BinarySearch(eIndex_o-sIndex_o,colIndex_old+sIndex_o,icol_left)+sIndex_o;
for (Int_t index = sIndex_o; index <= left; index++) {
rowIndex_new[irow+1]++;
colIndex_new[nr] = colIndex_old[index];
elements_new[nr] = elements_old[index];
nr++;
}
if (rowIndex_s && colIndex_s) {
const Int_t sIndex_s = rowIndex_s[irow-row_off];
const Int_t eIndex_s = rowIndex_s[irow-row_off+1];
for (Int_t index = sIndex_s; index < eIndex_s; index++) {
rowIndex_new[irow+1]++;
colIndex_new[nr] = colIndex_s[index]+col_off;
elements_new[nr] = elements_s[index];
nr++;
}
} else {
const Int_t off = (irow-row_off)*nCols_source;
for (Int_t icol = 0; icol < nCols_source; icol++) {
rowIndex_new[irow+1]++;
colIndex_new[nr] = icol+col_off;
elements_new[nr] = elements_s[off+icol];
nr++;
}
}
const Int_t icol_right = col_off+nCols_source-1;
if (colIndex_old) {
Int_t right = TMath::BinarySearch(eIndex_o-sIndex_o,colIndex_old+sIndex_o,icol_right)+sIndex_o;
while (right < eIndex_o-1 && colIndex_old[right+1] <= icol_right)
right++;
right++;
for (Int_t index = right; index < eIndex_o; index++) {
rowIndex_new[irow+1]++;
colIndex_new[nr] = colIndex_old[index];
elements_new[nr] = elements_old[index];
nr++;
}
}
} else {
for (Int_t index = sIndex_o; index < eIndex_o; index++) {
rowIndex_new[irow+1]++;
colIndex_new[nr] = colIndex_old[index];
elements_new[nr] = elements_old[index];
nr++;
}
}
}
R__ASSERT(this->fNelems == fRowIndex[this->fNrowIndex-1]);
if (rowIndex_old) delete [] (Int_t*) rowIndex_old;
if (colIndex_old) delete [] (Int_t*) colIndex_old;
if (elements_old) delete [] (Element*) elements_old;
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::Transpose(const TMatrixTSparse<Element> &source)
{
if (gMatrixCheck) {
R__ASSERT(this->IsValid());
R__ASSERT(source.IsValid());
if (this->fNrows != source.GetNcols() || this->fNcols != source.GetNrows() ||
this->fRowLwb != source.GetColLwb() || this->fColLwb != source.GetRowLwb()) {
Error("Transpose","matrix has wrong shape");
return *this;
}
if (source.NonZeros() <= 0)
return *this;
}
const Int_t nr_nonzeros = source.NonZeros();
const Int_t * const pRowIndex_s = source.GetRowIndexArray();
const Int_t * const pColIndex_s = source.GetColIndexArray();
const Element * const pData_s = source.GetMatrixArray();
Int_t *rownr = new Int_t [nr_nonzeros];
Int_t *colnr = new Int_t [nr_nonzeros];
Element *pData_t = new Element[nr_nonzeros];
Int_t ielem = 0;
for (Int_t irow_s = 0; irow_s < source.GetNrows(); irow_s++) {
const Int_t sIndex = pRowIndex_s[irow_s];
const Int_t eIndex = pRowIndex_s[irow_s+1];
for (Int_t index = sIndex; index < eIndex; index++) {
rownr[ielem] = pColIndex_s[index]+this->fRowLwb;
colnr[ielem] = irow_s+this->fColLwb;
pData_t[ielem] = pData_s[index];
ielem++;
}
}
R__ASSERT(nr_nonzeros >= ielem);
DoubleLexSort(nr_nonzeros,rownr,colnr,pData_t);
SetMatrixArray(nr_nonzeros,rownr,colnr,pData_t);
R__ASSERT(this->fNelems == fRowIndex[this->fNrowIndex-1]);
if (pData_t) delete [] pData_t;
if (rownr) delete [] rownr;
if (colnr) delete [] colnr;
return *this;
}
template<class Element>
TMatrixTBase<Element> &TMatrixTSparse<Element>::Zero()
{
R__ASSERT(this->IsValid());
if (fElements) { delete [] fElements; fElements = 0; }
if (fColIndex) { delete [] fColIndex; fColIndex = 0; }
this->fNelems = 0;
memset(this->GetRowIndexArray(),0,this->fNrowIndex*sizeof(Int_t));
return *this;
}
template<class Element>
TMatrixTBase<Element> &TMatrixTSparse<Element>::UnitMatrix()
{
R__ASSERT(this->IsValid());
Int_t i;
Int_t nr_nonzeros = 0;
for (i = this->fRowLwb; i <= this->fRowLwb+this->fNrows-1; i++)
for (Int_t j = this->fColLwb; j <= this->fColLwb+this->fNcols-1; j++)
if (i == j) nr_nonzeros++;
if (nr_nonzeros != this->fNelems) {
this->fNelems = nr_nonzeros;
Int_t *oIp = fColIndex;
fColIndex = new Int_t[nr_nonzeros];
if (oIp) delete [] oIp;
Element *oDp = fElements;
fElements = new Element[nr_nonzeros];
if (oDp) delete [] oDp;
}
Int_t ielem = 0;
fRowIndex[0] = 0;
for (i = this->fRowLwb; i <= this->fRowLwb+this->fNrows-1; i++) {
for (Int_t j = this->fColLwb; j <= this->fColLwb+this->fNcols-1; j++) {
if (i == j) {
const Int_t irow = i-this->fRowLwb;
fRowIndex[irow+1] = ielem+1;
fElements[ielem] = 1.0;
fColIndex[ielem++] = j-this->fColLwb;
}
}
}
return *this;
}
template<class Element>
Element TMatrixTSparse<Element>::RowNorm() const
{
R__ASSERT(this->IsValid());
const Element * ep = GetMatrixArray();
const Element * const fp = ep+this->fNelems;
const Int_t * const pR = GetRowIndexArray();
Element norm = 0;
for (Int_t irow = 0; irow < this->fNrows; irow++) {
const Int_t sIndex = pR[irow];
const Int_t eIndex = pR[irow+1];
Element sum = 0;
for (Int_t index = sIndex; index < eIndex; index++)
sum += TMath::Abs(*ep++);
norm = TMath::Max(norm,sum);
}
R__ASSERT(ep == fp);
return norm;
}
template<class Element>
Element TMatrixTSparse<Element>::ColNorm() const
{
R__ASSERT(this->IsValid());
const TMatrixTSparse<Element> mt(kTransposed,*this);
const Element * ep = mt.GetMatrixArray();
const Element * const fp = ep+this->fNcols;
Element norm = 0;
while (ep < fp) {
Element sum = 0;
for (Int_t i = 0; i < this->fNrows; i++,ep += this->fNcols)
sum += TMath::Abs(*ep);
ep -= this->fNelems-1;
norm = TMath::Max(norm,sum);
}
R__ASSERT(ep == fp);
return norm;
}
template<class Element>
Element &TMatrixTSparse<Element>::operator()(Int_t rown,Int_t coln)
{
R__ASSERT(this->IsValid());
const Int_t arown = rown-this->fRowLwb;
const Int_t acoln = coln-this->fColLwb;
if (arown >= this->fNrows || arown < 0) {
Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
return fElements[0];
}
if (acoln >= this->fNcols || acoln < 0) {
Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
return fElements[0];
}
Int_t index = -1;
Int_t sIndex = 0;
Int_t eIndex = 0;
if (this->fNrowIndex > 0 && fRowIndex[this->fNrowIndex-1] != 0) {
sIndex = fRowIndex[arown];
eIndex = fRowIndex[arown+1];
index = TMath::BinarySearch(eIndex-sIndex,fColIndex+sIndex,acoln)+sIndex;
}
if (index >= sIndex && fColIndex[index] == acoln)
return fElements[index];
else {
Element val = 0.;
InsertRow(rown,coln,&val,1);
sIndex = fRowIndex[arown];
eIndex = fRowIndex[arown+1];
index = TMath::BinarySearch(eIndex-sIndex,fColIndex+sIndex,acoln)+sIndex;
if (index >= sIndex && fColIndex[index] == acoln)
return fElements[index];
else {
Error("operator()(Int_t,Int_t","Insert row failed");
return fElements[0];
}
}
}
template <class Element>
Element TMatrixTSparse<Element>::operator()(Int_t rown,Int_t coln) const
{
R__ASSERT(this->IsValid());
if (this->fNrowIndex > 0 && this->fRowIndex[this->fNrowIndex-1] == 0) {
Error("operator()(Int_t,Int_t) const","row/col indices are not set");
Info("operator()","fNrowIndex = %d fRowIndex[fNrowIndex-1] = %d\n",this->fNrowIndex,this->fRowIndex[this->fNrowIndex-1]);
return 0.0;
}
const Int_t arown = rown-this->fRowLwb;
const Int_t acoln = coln-this->fColLwb;
if (arown >= this->fNrows || arown < 0) {
Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
return 0.0;
}
if (acoln >= this->fNcols || acoln < 0) {
Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
return 0.0;
}
const Int_t sIndex = fRowIndex[arown];
const Int_t eIndex = fRowIndex[arown+1];
const Int_t index = (Int_t)TMath::BinarySearch(eIndex-sIndex,fColIndex+sIndex,acoln)+sIndex;
if (index < sIndex || fColIndex[index] != acoln) return 0.0;
else return fElements[index];
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator=(const TMatrixTSparse<Element> &source)
{
if (gMatrixCheck && !AreCompatible(*this,source)) {
Error("operator=(const TMatrixTSparse &)","matrices not compatible");
return *this;
}
if (this->GetMatrixArray() != source.GetMatrixArray()) {
TObject::operator=(source);
const Element * const sp = source.GetMatrixArray();
Element * const tp = this->GetMatrixArray();
memcpy(tp,sp,this->fNelems*sizeof(Element));
this->fTol = source.GetTol();
}
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator=(const TMatrixT<Element> &source)
{
if (gMatrixCheck && !AreCompatible(*this,(TMatrixTBase<Element> &)source)) {
Error("operator=(const TMatrixT &)","matrices not compatible");
return *this;
}
if (this->GetMatrixArray() != source.GetMatrixArray()) {
TObject::operator=(source);
const Element * const sp = source.GetMatrixArray();
Element * const tp = this->GetMatrixArray();
for (Int_t irow = 0; irow < this->fNrows; irow++) {
const Int_t sIndex = fRowIndex[irow];
const Int_t eIndex = fRowIndex[irow+1];
const Int_t off = irow*this->fNcols;
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = fColIndex[index];
tp[index] = sp[off+icol];
}
}
this->fTol = source.GetTol();
}
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator=(Element val)
{
R__ASSERT(this->IsValid());
if (fRowIndex[this->fNrowIndex-1] == 0) {
Error("operator=(Element","row/col indices are not set");
return *this;
}
Element *ep = this->GetMatrixArray();
const Element * const ep_last = ep+this->fNelems;
while (ep < ep_last)
*ep++ = val;
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator+=(Element val)
{
R__ASSERT(this->IsValid());
Element *ep = this->GetMatrixArray();
const Element * const ep_last = ep+this->fNelems;
while (ep < ep_last)
*ep++ += val;
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator-=(Element val)
{
R__ASSERT(this->IsValid());
Element *ep = this->GetMatrixArray();
const Element * const ep_last = ep+this->fNelems;
while (ep < ep_last)
*ep++ -= val;
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::operator*=(Element val)
{
R__ASSERT(this->IsValid());
Element *ep = this->GetMatrixArray();
const Element * const ep_last = ep+this->fNelems;
while (ep < ep_last)
*ep++ *= val;
return *this;
}
template<class Element>
TMatrixTBase<Element> &TMatrixTSparse<Element>::Randomize(Element alpha,Element beta,Double_t &seed)
{
R__ASSERT(this->IsValid());
const Element scale = beta-alpha;
const Element shift = alpha/scale;
Int_t * const pRowIndex = GetRowIndexArray();
Int_t * const pColIndex = GetColIndexArray();
Element * const ep = GetMatrixArray();
const Int_t m = this->GetNrows();
const Int_t n = this->GetNcols();
const Int_t nn = this->GetNrows()*this->GetNcols();
const Int_t length = (this->GetNoElements() <= nn) ? this->GetNoElements() : nn;
Int_t chosen = 0;
Int_t icurrent = 0;
pRowIndex[0] = 0;
for (Int_t k = 0; k < nn; k++) {
const Element r = Drand(seed);
if ((nn-k)*r < length-chosen) {
pColIndex[chosen] = k%n;
const Int_t irow = k/n;
if (irow > icurrent) {
for ( ; icurrent < irow; icurrent++)
pRowIndex[icurrent+1] = chosen;
}
ep[chosen] = scale*(Drand(seed)+shift);
chosen++;
}
}
for ( ; icurrent < m; icurrent++)
pRowIndex[icurrent+1] = length;
R__ASSERT(chosen == length);
return *this;
}
template<class Element>
TMatrixTSparse<Element> &TMatrixTSparse<Element>::RandomizePD(Element alpha,Element beta,Double_t &seed)
{
const Element scale = beta-alpha;
const Element shift = alpha/scale;
if (gMatrixCheck) {
R__ASSERT(this->IsValid());
if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
Error("RandomizePD(Element &","matrix should be square");
return *this;
}
}
const Int_t n = this->fNcols;
Int_t * const pRowIndex = this->GetRowIndexArray();
Int_t * const pColIndex = this->GetColIndexArray();
Element * const ep = GetMatrixArray();
pRowIndex[0] = 0;
pColIndex[0] = 0;
pRowIndex[1] = 1;
ep[0] = 1e-8+scale*(Drand(seed)+shift);
const Int_t nn = n*(n-1)/2;
Int_t length = this->fNelems-n;
length = (length <= nn) ? length : nn;
Int_t chosen = 0;
Int_t icurrent = 1;
Int_t nnz = 1;
for (Int_t k = 0; k < nn; k++ ) {
const Element r = Drand(seed);
if( (nn-k)*r < length-chosen) {
Int_t row = (int) TMath::Floor((-1+TMath::Sqrt(1.0+8.0*k))/2);
Int_t col = k-row*(row+1)/2;
row++;
if (row > icurrent) {
for ( ; icurrent < row; icurrent++) {
ep[nnz] = 0.0;
for (Int_t ll = pRowIndex[icurrent]; ll < nnz; ll++)
ep[nnz] += TMath::Abs(ep[ll]);
ep[nnz] += 1e-8+scale*(Drand(seed)+shift);
pColIndex[nnz] = icurrent;
nnz++;
pRowIndex[icurrent+1] = nnz;
}
}
ep[nnz] = scale*(Drand(seed)+shift);
pColIndex[nnz] = col;
ep[pRowIndex[col+1]-1] += TMath::Abs(ep[nnz]);
nnz++;
chosen++;
}
}
R__ASSERT(chosen == length);
for ( ; icurrent < n; icurrent++) {
ep[nnz] = 0.0;
for(Int_t ll = pRowIndex[icurrent]; ll < nnz; ll++)
ep[nnz] += TMath::Abs(ep[ll]);
ep[nnz] += 1e-8+scale*(Drand(seed)+shift);
pColIndex[nnz] = icurrent;
nnz++;
pRowIndex[icurrent+1] = nnz;
}
this->fNelems = nnz;
TMatrixTSparse<Element> tmp(TMatrixTSparse<Element>::kTransposed,*this);
*this += tmp;
{
const Int_t * const pR = this->GetRowIndexArray();
const Int_t * const pC = this->GetColIndexArray();
Element * const pD = GetMatrixArray();
for (Int_t irow = 0; irow < this->fNrows+1; irow++) {
const Int_t sIndex = pR[irow];
const Int_t eIndex = pR[irow+1];
for (Int_t index = sIndex; index < eIndex; index++) {
const Int_t icol = pC[index];
if (irow == icol)
pD[index] /= 2.;
}
}
}
return *this;
}
template<class Element>
TMatrixTSparse<Element> operator+(const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2)
{
TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kPlus,source2);
return target;
}
template<class Element>
TMatrixTSparse<Element> operator+(const TMatrixTSparse<Element> &source1,const TMatrixT<Element> &source2)
{
TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kPlus,source2);
return target;
}
template<class Element>
TMatrixTSparse<Element> operator+(const TMatrixT<Element> &source1,const TMatrixTSparse<Element> &source2)
{
TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kPlus,source2);
return target;
}
template<class Element>
TMatrixTSparse<Element> operator+(const TMatrixTSparse<Element> &source,Element val)
{
TMatrixTSparse<Element> target(source);
target += val;
return target;
}
template<class Element>
TMatrixTSparse<Element> operator+(Element val,const TMatrixTSparse<Element> &source)
{
TMatrixTSparse<Element> target(source);
target += val;
return target;
}
template<class Element>
TMatrixTSparse<Element> operator-(const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2)
{
TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMinus,source2);
return target;
}
template<class Element>
TMatrixTSparse<Element> operator-(const TMatrixTSparse<Element> &source1,const TMatrixT<Element> &source2)
{
TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMinus,source2);
return target;
}
template<class Element>
TMatrixTSparse<Element> operator-(const TMatrixT<Element> &source1,const TMatrixTSparse<Element> &source2)
{
TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMinus,source2);
return target;
}
template<class Element>
TMatrixTSparse<Element> operator-(const TMatrixTSparse<Element> &source,Element val)
{
TMatrixTSparse<Element> target(source);
target -= val;
return target;
}
template<class Element>
TMatrixTSparse<Element> operator-(Element val,const TMatrixTSparse<Element> &source)
{
TMatrixTSparse<Element> target(source);
target -= val;
return target;
}
template<class Element>
TMatrixTSparse<Element> operator*(const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2)
{
TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMult,source2);
return target;
}
template<class Element>
TMatrixTSparse<Element> operator*(const TMatrixTSparse<Element> &source1,const TMatrixT<Element> &source2)
{
TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMult,source2);
return target;
}
template<class Element>
TMatrixTSparse<Element> operator*(const TMatrixT<Element> &source1,const TMatrixTSparse<Element> &source2)
{
TMatrixTSparse<Element> target(source1,TMatrixTSparse<Element>::kMult,source2);
return target;
}
template<class Element>
TMatrixTSparse<Element> operator*(Element val,const TMatrixTSparse<Element> &source)
{
TMatrixTSparse<Element> target(source);
target *= val;
return target;
}
template<class Element>
TMatrixTSparse<Element> operator*(const TMatrixTSparse<Element> &source,Element val)
{
TMatrixTSparse<Element> target(source);
target *= val;
return target;
}
template<class Element>
TMatrixTSparse<Element> &Add(TMatrixTSparse<Element> &target,Element scalar,const TMatrixTSparse<Element> &source)
{
target += scalar * source;
return target;
}
template<class Element>
TMatrixTSparse<Element> &ElementMult(TMatrixTSparse<Element> &target,const TMatrixTSparse<Element> &source)
{
if (gMatrixCheck && !AreCompatible(target,source)) {
::Error("ElementMult(TMatrixTSparse &,const TMatrixTSparse &)","matrices not compatible");
return target;
}
const Element *sp = source.GetMatrixArray();
Element *tp = target.GetMatrixArray();
const Element *ftp = tp+target.GetNoElements();
while ( tp < ftp )
*tp++ *= *sp++;
return target;
}
template<class Element>
TMatrixTSparse<Element> &ElementDiv(TMatrixTSparse<Element> &target,const TMatrixTSparse<Element> &source)
{
if (gMatrixCheck && !AreCompatible(target,source)) {
::Error("ElementDiv(TMatrixT &,const TMatrixT &)","matrices not compatible");
return target;
}
const Element *sp = source.GetMatrixArray();
Element *tp = target.GetMatrixArray();
const Element *ftp = tp+target.GetNoElements();
while ( tp < ftp ) {
if (*sp != 0.0)
*tp++ /= *sp++;
else {
Error("ElementDiv","source element is zero");
tp++;
}
}
return target;
}
template<class Element>
Bool_t AreCompatible(const TMatrixTSparse<Element> &m1,const TMatrixTSparse<Element> &m2,Int_t verbose)
{
if (!m1.IsValid()) {
if (verbose)
::Error("AreCompatible", "matrix 1 not valid");
return kFALSE;
}
if (!m2.IsValid()) {
if (verbose)
::Error("AreCompatible", "matrix 2 not valid");
return kFALSE;
}
if (m1.GetNrows() != m2.GetNrows() || m1.GetNcols() != m2.GetNcols() ||
m1.GetRowLwb() != m2.GetRowLwb() || m1.GetColLwb() != m2.GetColLwb()) {
if (verbose)
::Error("AreCompatible", "matrices 1 and 2 not compatible");
return kFALSE;
}
const Int_t *pR1 = m1.GetRowIndexArray();
const Int_t *pR2 = m2.GetRowIndexArray();
const Int_t nRows = m1.GetNrows();
if (memcmp(pR1,pR2,(nRows+1)*sizeof(Int_t))) {
if (verbose)
::Error("AreCompatible", "matrices 1 and 2 have different rowIndex");
for (Int_t i = 0; i < nRows+1; i++)
printf("%d: %d %d\n",i,pR1[i],pR2[i]);
return kFALSE;
}
const Int_t *pD1 = m1.GetColIndexArray();
const Int_t *pD2 = m2.GetColIndexArray();
const Int_t nData = m1.GetNoElements();
if (memcmp(pD1,pD2,nData*sizeof(Int_t))) {
if (verbose)
::Error("AreCompatible", "matrices 1 and 2 have different colIndex");
for (Int_t i = 0; i < nData; i++)
printf("%d: %d %d\n",i,pD1[i],pD2[i]);
return kFALSE;
}
return kTRUE;
}
template<class Element>
void TMatrixTSparse<Element>::Streamer(TBuffer &R__b)
{
if (R__b.IsReading()) {
UInt_t R__s, R__c;
Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
Clear();
R__b.ReadClassBuffer(TMatrixTSparse<Element>::Class(),this,R__v,R__s,R__c);
if (this->fNelems < 0)
this->Invalidate();
} else {
R__b.WriteClassBuffer(TMatrixTSparse<Element>::Class(),this);
}
}
template class TMatrixTSparse<Float_t>;
#ifndef ROOT_TMatrixFSparsefwd
#include "TMatrixFSparsefwd.h"
#endif
#ifndef ROOT_TMatrixFfwd
#include "TMatrixFfwd.h"
#endif
template TMatrixFSparse operator+ <Float_t>(const TMatrixFSparse &source1,const TMatrixFSparse &source2);
template TMatrixFSparse operator+ <Float_t>(const TMatrixFSparse &source1,const TMatrixF &source2);
template TMatrixFSparse operator+ <Float_t>(const TMatrixF &source1,const TMatrixFSparse &source2);
template TMatrixFSparse operator+ <Float_t>(const TMatrixFSparse &source , Float_t val );
template TMatrixFSparse operator+ <Float_t>( Float_t val ,const TMatrixFSparse &source );
template TMatrixFSparse operator- <Float_t>(const TMatrixFSparse &source1,const TMatrixFSparse &source2);
template TMatrixFSparse operator- <Float_t>(const TMatrixFSparse &source1,const TMatrixF &source2);
template TMatrixFSparse operator- <Float_t>(const TMatrixF &source1,const TMatrixFSparse &source2);
template TMatrixFSparse operator- <Float_t>(const TMatrixFSparse &source , Float_t val );
template TMatrixFSparse operator- <Float_t>( Float_t val ,const TMatrixFSparse &source );
template TMatrixFSparse operator* <Float_t>(const TMatrixFSparse &source1,const TMatrixFSparse &source2);
template TMatrixFSparse operator* <Float_t>(const TMatrixFSparse &source1,const TMatrixF &source2);
template TMatrixFSparse operator* <Float_t>(const TMatrixF &source1,const TMatrixFSparse &source2);
template TMatrixFSparse operator* <Float_t>( Float_t val ,const TMatrixFSparse &source );
template TMatrixFSparse operator* <Float_t>(const TMatrixFSparse &source, Float_t val );
template TMatrixFSparse &Add <Float_t>(TMatrixFSparse &target, Float_t scalar,const TMatrixFSparse &source);
template TMatrixFSparse &ElementMult <Float_t>(TMatrixFSparse &target,const TMatrixFSparse &source);
template TMatrixFSparse &ElementDiv <Float_t>(TMatrixFSparse &target,const TMatrixFSparse &source);
template Bool_t AreCompatible<Float_t>(const TMatrixFSparse &m1,const TMatrixFSparse &m2,Int_t verbose);
#ifndef ROOT_TMatrixDSparsefwd
#include "TMatrixDSparsefwd.h"
#endif
#ifndef ROOT_TMatrixDfwd
#include "TMatrixDfwd.h"
#endif
template class TMatrixTSparse<Double_t>;
template TMatrixDSparse operator+ <Double_t>(const TMatrixDSparse &source1,const TMatrixDSparse &source2);
template TMatrixDSparse operator+ <Double_t>(const TMatrixDSparse &source1,const TMatrixD &source2);
template TMatrixDSparse operator+ <Double_t>(const TMatrixD &source1,const TMatrixDSparse &source2);
template TMatrixDSparse operator+ <Double_t>(const TMatrixDSparse &source , Double_t val );
template TMatrixDSparse operator+ <Double_t>( Double_t val ,const TMatrixDSparse &source );
template TMatrixDSparse operator- <Double_t>(const TMatrixDSparse &source1,const TMatrixDSparse &source2);
template TMatrixDSparse operator- <Double_t>(const TMatrixDSparse &source1,const TMatrixD &source2);
template TMatrixDSparse operator- <Double_t>(const TMatrixD &source1,const TMatrixDSparse &source2);
template TMatrixDSparse operator- <Double_t>(const TMatrixDSparse &source , Double_t val );
template TMatrixDSparse operator- <Double_t>( Double_t val ,const TMatrixDSparse &source );
template TMatrixDSparse operator* <Double_t>(const TMatrixDSparse &source1,const TMatrixDSparse &source2);
template TMatrixDSparse operator* <Double_t>(const TMatrixDSparse &source1,const TMatrixD &source2);
template TMatrixDSparse operator* <Double_t>(const TMatrixD &source1,const TMatrixDSparse &source2);
template TMatrixDSparse operator* <Double_t>( Double_t val ,const TMatrixDSparse &source );
template TMatrixDSparse operator* <Double_t>(const TMatrixDSparse &source, Double_t val );
template TMatrixDSparse &Add <Double_t>(TMatrixDSparse &target, Double_t scalar,const TMatrixDSparse &source);
template TMatrixDSparse &ElementMult <Double_t>(TMatrixDSparse &target,const TMatrixDSparse &source);
template TMatrixDSparse &ElementDiv <Double_t>(TMatrixDSparse &target,const TMatrixDSparse &source);
template Bool_t AreCompatible<Double_t>(const TMatrixDSparse &m1,const TMatrixDSparse &m2,Int_t verbose);