// @(#)root/matrix:$Id$ // Authors: Fons Rademakers, Eddy Offermann Feb 2004 /************************************************************************* * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. * * All rights reserved. * * * * For the licensing terms see $ROOTSYS/LICENSE. * * For the list of contributors see $ROOTSYS/README/CREDITS. * *************************************************************************/ ////////////////////////////////////////////////////////////////////////// // // // TMatrixTSparse // // // // Template class of a general sparse matrix in the Harwell-Boeing // // format // // // // Besides the usual shape/size decsriptors of a matrix like fNrows, // // fRowLwb,fNcols and fColLwb, we also store a row index, fRowIndex and // // column index, fColIndex only for those elements unequal zero: // // // // fRowIndex[0,..,fNrows]: Stores for each row the index range of // // the elements in the data and column array // // fColIndex[0,..,fNelems-1]: Stores the column number for each data // // element != 0 // // // // As an example how to access all sparse data elements: // // // // for (Int_t irow = 0; irow < this->fNrows; irow++) { // // 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]; // // const Element data = fElements[index]; // // printf("data(%d,%d) = %.4e\n",irow+this->fRowLwb,icol+ // // this->fColLwb,data); // // } // // } // // // // When checking whether sparse matrices are compatible (like in an // // assigment !), not only the shape parameters are compared but also // // the sparse structure through fRowIndex and fColIndex . // // // // Several methods exist to fill a sparse matrix with data entries. // // Most are the same like for dense matrices but some care has to be // // taken with regard to performance. In the constructor, always the // // shape of the matrix has to be specified in some form . Data can be // // entered through the following methods : // // 1. constructor // // TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t dol_lwb, // // Int_t col_upb,Int_t nr_nonzeros, // // Int_t *row, Int_t *col,Element *data); // // It uses SetMatrixArray(..), see below // // 2. copy constructors // // 3. SetMatrixArray(Int_t nr,Int_t *irow,Int_t *icol,Element *data) // // where it is expected that the irow,icol and data array contain // // nr entries . Only the entries with non-zero data[i] value are // // inserted. Be aware that the input data array will be modified // // inside the routine for doing the necessary sorting of indices ! // // 4. TMatrixTSparse a(n,m); for(....) { a(i,j) = .... // // This is a very flexible method but expensive : // // - if no entry for slot (i,j) is found in the sparse index table // // it will be entered, which involves some memory management ! // // - before invoking this method in a loop it is smart to first // // set the index table through a call to SetSparseIndex(..) // // 5. SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixTBase &source) // // the matrix to be inserted at position (row_lwb,col_lwb) can be // // both dense or sparse . // // // ////////////////////////////////////////////////////////////////////////// #include "TMatrixTSparse.h" #include "TMatrixT.h" #include "TMath.h" #ifndef R__ALPHA templateClassImp(TMatrixTSparse) #endif //______________________________________________________________________________ template TMatrixTSparse::TMatrixTSparse(Int_t no_rows,Int_t no_cols) { // Space is allocated for row/column indices and data, but the sparse structure // information has still to be set ! Allocate(no_rows,no_cols,0,0,1); } //______________________________________________________________________________ template TMatrixTSparse::TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb) { // Space is allocated for row/column indices and data, but the sparse structure // information has still to be set ! Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb,1); } //______________________________________________________________________________ template TMatrixTSparse::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) { // Space is allocated for row/column indices and data. Sparse row/column index // structure together with data is coming from the arrays, row, col and data, resp . 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 TMatrixTSparse::TMatrixTSparse(const TMatrixTSparse &another) : TMatrixTBase(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 TMatrixTSparse::TMatrixTSparse(const TMatrixT &another) : TMatrixTBase(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 TMatrixTSparse::TMatrixTSparse(EMatrixCreatorsOp1 op,const TMatrixTSparse &prototype) { // Create a matrix applying a specific operation to the prototype. // Supported operations are: kZero, kUnit, kTransposed and kAtA 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 at(TMatrixTSparse::kTransposed,prototype); AMultBt(at,at,1); break; } default: Error("TMatrixTSparse(EMatrixCreatorOp1)","operation %d not yet implemented", op); } } //______________________________________________________________________________ template TMatrixTSparse::TMatrixTSparse(const TMatrixTSparse &a,EMatrixCreatorsOp2 op,const TMatrixTSparse &b) { // Create a matrix applying a specific operation to two prototypes. // Supported operations are: kMult (a*b), kMultTranspose (a*b'), kPlus (a+b), kMinus (a-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 TMatrixTSparse::TMatrixTSparse(const TMatrixTSparse &a,EMatrixCreatorsOp2 op,const TMatrixT &b) { // Create a matrix applying a specific operation to two prototypes. // Supported operations are: kMult (a*b), kMultTranspose (a*b'), kPlus (a+b), kMinus (a-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 TMatrixTSparse::TMatrixTSparse(const TMatrixT &a,EMatrixCreatorsOp2 op,const TMatrixTSparse &b) { // Create a matrix applying a specific operation to two prototypes. // Supported operations are: kMult (a*b), kMultTranspose (a*b'), kPlus (a+b), kMinus (a-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 void TMatrixTSparse::Allocate(Int_t no_rows,Int_t no_cols,Int_t row_lwb,Int_t col_lwb, Int_t init,Int_t nr_nonzeros) { // Allocate new matrix. Arguments are number of rows, columns, row lowerbound (0 default) // and column lowerbound (0 default), 0 initialization flag and number of non-zero // elements (only relevant for sparse format). 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::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 TMatrixTBase &TMatrixTSparse::InsertRow(Int_t rown,Int_t coln,const Element *v,Int_t n) { // Insert in row rown, n elements of array v at column coln 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]; // check first how many slots are available from [acoln,..,acoln+nr-1] // also note lIndex and rIndex so that [sIndex..lIndex] and [rIndex..eIndex-1] // contain the row entries except for the region to be inserted 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 void TMatrixTSparse::ExtractRow(Int_t rown, Int_t coln, Element *v,Int_t n) const { // Store in array v, n matrix elements of row rown starting at column coln 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 void TMatrixTSparse::AMultBt(const TMatrixTSparse &a,const TMatrixTSparse &b,Int_t constr) { // General matrix multiplication. Create a matrix C such that C = A * B'. // Note, matrix C is allocated for constr=1. 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) { // make a best guess of the sparse structure; it will guarantee // enough allocated space ! 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; // best guess 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 void TMatrixTSparse::AMultBt(const TMatrixTSparse &a,const TMatrixT &b,Int_t constr) { // General matrix multiplication. Create a matrix C such that C = A * B'. // Note, matrix C is allocated for constr=1. 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) { // make a best guess of the sparse structure; it will guarantee // enough allocated space ! 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; // best guess 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 void TMatrixTSparse::AMultBt(const TMatrixT &a,const TMatrixTSparse &b,Int_t constr) { // General matrix multiplication. Create a matrix C such that C = A * B'. // Note, matrix C is allocated for constr=1. 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) { // make a best guess of the sparse structure; it will guarantee // enough allocated space ! 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; // best guess 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 void TMatrixTSparse::APlusB(const TMatrixTSparse &a,const TMatrixTSparse &b,Int_t constr) { // General matrix addition. Create a matrix C such that C = A + B. // Note, matrix C is allocated for constr=1. 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 void TMatrixTSparse::APlusB(const TMatrixTSparse &a,const TMatrixT &b,Int_t constr) { // General matrix addition. Create a matrix C such that C = A + B. // Note, matrix C is allocated for constr=1. 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 void TMatrixTSparse::AMinusB(const TMatrixTSparse &a,const TMatrixTSparse &b,Int_t constr) { // General matrix subtraction. Create a matrix C such that C = A - B. // Note, matrix C is allocated for constr=1. 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 void TMatrixTSparse::AMinusB(const TMatrixTSparse &a,const TMatrixT &b,Int_t constr) { // General matrix subtraction. Create a matrix C such that C = A - B. // Note, matrix C is allocated for constr=1. 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 void TMatrixTSparse::AMinusB(const TMatrixT &a,const TMatrixTSparse &b,Int_t constr) { // General matrix subtraction. Create a matrix C such that C = A - B. // Note, matrix C is allocated for constr=1. 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 void TMatrixTSparse::GetMatrix2Array(Element *data,Option_t * /*option*/) const { // Copy matrix data to array . It is assumed that array is of size >= fNelems R__ASSERT(this->IsValid()); const Element * const elem = GetMatrixArray(); memcpy(data,elem,this->fNelems*sizeof(Element)); } //______________________________________________________________________________ template TMatrixTBase &TMatrixTSparse::SetMatrixArray(Int_t nr,Int_t *row,Int_t *col,Element *data) { // Copy nr elements from row/col index and data array to matrix . It is assumed // that arrays are of size >= nr // Note that the input arrays are not passed as const since they will be modified ! 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; } } TMatrixTBase::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 TMatrixTSparse &TMatrixTSparse::SetSparseIndex(Int_t nelems_new) { // Increase/decrease the number of non-zero elements to 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]; if (oIp) { memmove(fColIndex,oIp,nr*sizeof(Int_t)); delete [] oIp; } Element *oDp = fElements; fElements = new Element[nelems_new]; if (oDp) { memmove(fElements,oDp,nr*sizeof(Element)); 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 TMatrixTSparse &TMatrixTSparse::SetSparseIndex(const TMatrixTBase &source) { // Use non-zero data of matrix source to set the sparse structure 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 TMatrixTSparse &TMatrixTSparse::SetSparseIndexAB(const TMatrixTSparse &a,const TMatrixTSparse &b) { // Set the row/column indices to the "sum" of matrices a and b // It is checked that enough space has been allocated 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(); // Count first the number of non-zero slots that are needed 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 > sIndexa && eIndexa > 0) ? pColIndexa[eIndexa-1] : -1; if (pColIndexb[indexb++] > icola) nc++; } } // Allocate the necessary space in fRowIndex and fColIndex 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 TMatrixTSparse &TMatrixTSparse::SetSparseIndexAB(const TMatrixT &a,const TMatrixTSparse &b) { // Set the row/column indices to the "sum" of matrices a and b // It is checked that enough space has been allocated 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(); // Count first the number of non-zero slots that are needed 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; } } } } // Allocate thre necessary space in fRowIndex and fColIndex 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 TMatrixTBase &TMatrixTSparse::ResizeTo(Int_t nrows,Int_t ncols,Int_t nr_nonzeros) { // Set size of the matrix to nrows x ncols with nr_nonzeros non-zero entries // if nr_nonzeros > 0 . // New dynamic elements are created, the overlapping part of the old ones are // copied to the new structures, then the old elements are deleted. 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 TMatrixTBase &TMatrixTSparse::ResizeTo(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb, Int_t nr_nonzeros) { // Set size of the matrix to [row_lwb:row_upb] x [col_lwb:col_upb] with nr_nonzeros // non-zero entries if nr_nonzeros > 0 . // New dynamic elemenst are created, the overlapping part of the old ones are // copied to the new structures, then the old elements are deleted. 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; 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) { 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 TMatrixTSparse &TMatrixTSparse::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::epsilon(); fElements = pData; fRowIndex = pRowIndex; fColIndex = pColIndex; return *this; } //______________________________________________________________________________ template TMatrixTBase &TMatrixTSparse::GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb, TMatrixTBase &target,Option_t *option) const { // Get submatrix [row_lwb..row_upb][col_lwb..col_upb]; The indexing range of the // returned matrix depends on the argument option: // // option == "S" : return [0..row_upb-row_lwb+1][0..col_upb-col_lwb+1] (default) // else : return [row_lwb..row_upb][col_lwb..col_upb] 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 TMatrixTBase &TMatrixTSparse::SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixTBase &source) { // Insert matrix source starting at [row_lwb][col_lwb], thereby overwriting the part // [row_lwb..row_lwb+nrows_source-1][col_lwb..col_lwb+ncols_source-1]; 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(); // Determine how many non-zero's are already available in // [row_lwb..row_lwb+nrows_source-1][col_lwb..col_lwb+ncols_source-1] 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 TMatrixTSparse &TMatrixTSparse::Transpose(const TMatrixTSparse &source) { // Transpose a matrix. 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); TMatrixTBase::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 TMatrixTBase &TMatrixTSparse::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 TMatrixTBase &TMatrixTSparse::UnitMatrix() { // Make a unit matrix (matrix need not be a square one). 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 Element TMatrixTSparse::RowNorm() const { // Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}. // The norm is induced by the infinity vector norm. R__ASSERT(this->IsValid()); const Element * ep = GetMatrixArray(); const Element * const fp = ep+this->fNelems; const Int_t * const pR = GetRowIndexArray(); Element norm = 0; // Scan the matrix row-after-row 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 Element TMatrixTSparse::ColNorm() const { // Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}. // The norm is induced by the 1 vector norm. R__ASSERT(this->IsValid()); const TMatrixTSparse mt(kTransposed,*this); const Element * ep = mt.GetMatrixArray(); const Element * const fp = ep+this->fNcols; Element norm = 0; // Scan the matrix col-after-col while (ep < fp) { Element sum = 0; // Scan a col to compute the sum for (Int_t i = 0; i < this->fNrows; i++,ep += this->fNcols) sum += TMath::Abs(*ep); ep -= this->fNelems-1; // Point ep to the beginning of the next col norm = TMath::Max(norm,sum); } R__ASSERT(ep == fp); return norm; } //______________________________________________________________________________ template Element &TMatrixTSparse::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 Element TMatrixTSparse::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 TMatrixTSparse &TMatrixTSparse::operator=(const TMatrixTSparse &source) { // Notice that the sparsity of the matrix is NOT changed : its fRowIndex/fColIndex // are used ! 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 TMatrixTSparse &TMatrixTSparse::operator=(const TMatrixT &source) { // Notice that the sparsity of the matrix is NOT changed : its fRowIndex/fColIndex // are used ! if (gMatrixCheck && !AreCompatible(*this,(TMatrixTBase &)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 TMatrixTSparse &TMatrixTSparse::operator=(Element val) { // Assign val to every element of the matrix. Check that the row/col // indices are set ! 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 TMatrixTSparse &TMatrixTSparse::operator+=(Element val) { // Add val to every element of the matrix. 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 TMatrixTSparse &TMatrixTSparse::operator-=(Element val) { // Subtract val from every element of the matrix. 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 TMatrixTSparse &TMatrixTSparse::operator*=(Element val) { // Multiply every element of the matrix with 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 TMatrixTBase &TMatrixTSparse::Randomize(Element alpha,Element beta,Double_t &seed) { // randomize matrix element values 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(); // Knuth's algorithm for choosing "length" elements out of nn . 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 TMatrixTSparse &TMatrixTSparse::RandomizePD(Element alpha,Element beta,Double_t &seed) { // randomize matrix element values but keep matrix symmetric positive definite 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(); // We will always have non-zeros on the diagonal, so there // is no randomness there. In fact, choose the (0,0) element now pRowIndex[0] = 0; pColIndex[0] = 0; pRowIndex[1] = 1; ep[0] = 1e-8+scale*(Drand(seed)+shift); // Knuth's algorithm for choosing length elements out of nn . // nn here is the number of elements in the strict lower triangle. const Int_t nn = n*(n-1)/2; // length is the number of elements that can be stored, minus the number // of elements in the diagonal, which will always be in the matrix. // Int_t length = (this->fNelems-n)/2; Int_t length = this->fNelems-n; length = (length <= nn) ? length : nn; // chosen : the number of elements that have already been chosen (now 0) // nnz : the number of non-zeros in the matrix (now 1, because the // (0,0) element is already in the matrix. // icurrent : the index of the last row whose start has been stored in pRowIndex; 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) { // Element k is chosen. What row is it in? // In a lower triangular matrix (including a diagonal), it will be in // the largest row such that row*(row+1)/2 < k. In other words Int_t row = (int) TMath::Floor((-1+TMath::Sqrt(1.0+8.0*k))/2); // and its column will be the remainder Int_t col = k-row*(row+1)/2; // but since we are only filling in the *strict* lower triangle of // the matrix, we shift the row by 1 row++; if (row > icurrent) { // We have chosen a row beyond the current row. // Choose a diagonal element for each intermediate row and fix the // data structure. for ( ; icurrent < row; icurrent++) { // Choose the diagonal 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; } } // end if we have chosen a row beyond the current row; ep[nnz] = scale*(Drand(seed)+shift); pColIndex[nnz] = col; // add the value of this element (which occurs symmetrically in the // upper triangle) to the appropriate diagonal element ep[pRowIndex[col+1]-1] += TMath::Abs(ep[nnz]); nnz++; // We have added another element to the matrix chosen++; // And finished choosing another element. } } R__ASSERT(chosen == length); // and of course, we must choose all remaining diagonal elements . for ( ; icurrent < n; icurrent++) { // Choose the diagonal 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 tmp(TMatrixTSparse::kTransposed,*this); *this += tmp; // make sure to divide the diagonal by 2 becuase the operation // *this += tmp; adds the diagonal again { 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 TMatrixTSparse operator+(const TMatrixTSparse &source1,const TMatrixTSparse &source2) { TMatrixTSparse target(source1,TMatrixTSparse::kPlus,source2); return target; } //______________________________________________________________________________ template TMatrixTSparse operator+(const TMatrixTSparse &source1,const TMatrixT &source2) { TMatrixTSparse target(source1,TMatrixTSparse::kPlus,source2); return target; } //______________________________________________________________________________ template TMatrixTSparse operator+(const TMatrixT &source1,const TMatrixTSparse &source2) { TMatrixTSparse target(source1,TMatrixTSparse::kPlus,source2); return target; } //______________________________________________________________________________ template TMatrixTSparse operator+(const TMatrixTSparse &source,Element val) { TMatrixTSparse target(source); target += val; return target; } //______________________________________________________________________________ template TMatrixTSparse operator+(Element val,const TMatrixTSparse &source) { TMatrixTSparse target(source); target += val; return target; } //______________________________________________________________________________ template TMatrixTSparse operator-(const TMatrixTSparse &source1,const TMatrixTSparse &source2) { TMatrixTSparse target(source1,TMatrixTSparse::kMinus,source2); return target; } //______________________________________________________________________________ template TMatrixTSparse operator-(const TMatrixTSparse &source1,const TMatrixT &source2) { TMatrixTSparse target(source1,TMatrixTSparse::kMinus,source2); return target; } //______________________________________________________________________________ template TMatrixTSparse operator-(const TMatrixT &source1,const TMatrixTSparse &source2) { TMatrixTSparse target(source1,TMatrixTSparse::kMinus,source2); return target; } //______________________________________________________________________________ template TMatrixTSparse operator-(const TMatrixTSparse &source,Element val) { TMatrixTSparse target(source); target -= val; return target; } //______________________________________________________________________________ template TMatrixTSparse operator-(Element val,const TMatrixTSparse &source) { TMatrixTSparse target(source); target -= val; return target; } //______________________________________________________________________________ template TMatrixTSparse operator*(const TMatrixTSparse &source1,const TMatrixTSparse &source2) { TMatrixTSparse target(source1,TMatrixTSparse::kMult,source2); return target; } //______________________________________________________________________________ template TMatrixTSparse operator*(const TMatrixTSparse &source1,const TMatrixT &source2) { TMatrixTSparse target(source1,TMatrixTSparse::kMult,source2); return target; } //______________________________________________________________________________ template TMatrixTSparse operator*(const TMatrixT &source1,const TMatrixTSparse &source2) { TMatrixTSparse target(source1,TMatrixTSparse::kMult,source2); return target; } //______________________________________________________________________________ template TMatrixTSparse operator*(Element val,const TMatrixTSparse &source) { TMatrixTSparse target(source); target *= val; return target; } //______________________________________________________________________________ template TMatrixTSparse operator*(const TMatrixTSparse &source,Element val) { TMatrixTSparse target(source); target *= val; return target; } //______________________________________________________________________________ template TMatrixTSparse &Add(TMatrixTSparse &target,Element scalar,const TMatrixTSparse &source) { // Modify addition: target += scalar * source. target += scalar * source; return target; } //______________________________________________________________________________ template TMatrixTSparse &ElementMult(TMatrixTSparse &target,const TMatrixTSparse &source) { // Multiply target by the source, element-by-element. 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 TMatrixTSparse &ElementDiv(TMatrixTSparse &target,const TMatrixTSparse &source) { // Divide target by the source, element-by-element. 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 Bool_t AreCompatible(const TMatrixTSparse &m1,const TMatrixTSparse &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 void TMatrixTSparse::Streamer(TBuffer &R__b) { // Stream an object of class TMatrixTSparse. 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::Class(),this,R__v,R__s,R__c); if (this->fNelems < 0) this->Invalidate(); } else { R__b.WriteClassBuffer(TMatrixTSparse::Class(),this); } } template class TMatrixTSparse; #ifndef ROOT_TMatrixFSparsefwd #include "TMatrixFSparsefwd.h" #endif #ifndef ROOT_TMatrixFfwd #include "TMatrixFfwd.h" #endif template TMatrixFSparse operator+ (const TMatrixFSparse &source1,const TMatrixFSparse &source2); template TMatrixFSparse operator+ (const TMatrixFSparse &source1,const TMatrixF &source2); template TMatrixFSparse operator+ (const TMatrixF &source1,const TMatrixFSparse &source2); template TMatrixFSparse operator+ (const TMatrixFSparse &source , Float_t val ); template TMatrixFSparse operator+ ( Float_t val ,const TMatrixFSparse &source ); template TMatrixFSparse operator- (const TMatrixFSparse &source1,const TMatrixFSparse &source2); template TMatrixFSparse operator- (const TMatrixFSparse &source1,const TMatrixF &source2); template TMatrixFSparse operator- (const TMatrixF &source1,const TMatrixFSparse &source2); template TMatrixFSparse operator- (const TMatrixFSparse &source , Float_t val ); template TMatrixFSparse operator- ( Float_t val ,const TMatrixFSparse &source ); template TMatrixFSparse operator* (const TMatrixFSparse &source1,const TMatrixFSparse &source2); template TMatrixFSparse operator* (const TMatrixFSparse &source1,const TMatrixF &source2); template TMatrixFSparse operator* (const TMatrixF &source1,const TMatrixFSparse &source2); template TMatrixFSparse operator* ( Float_t val ,const TMatrixFSparse &source ); template TMatrixFSparse operator* (const TMatrixFSparse &source, Float_t val ); template TMatrixFSparse &Add (TMatrixFSparse &target, Float_t scalar,const TMatrixFSparse &source); template TMatrixFSparse &ElementMult (TMatrixFSparse &target,const TMatrixFSparse &source); template TMatrixFSparse &ElementDiv (TMatrixFSparse &target,const TMatrixFSparse &source); template Bool_t AreCompatible(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; template TMatrixDSparse operator+ (const TMatrixDSparse &source1,const TMatrixDSparse &source2); template TMatrixDSparse operator+ (const TMatrixDSparse &source1,const TMatrixD &source2); template TMatrixDSparse operator+ (const TMatrixD &source1,const TMatrixDSparse &source2); template TMatrixDSparse operator+ (const TMatrixDSparse &source , Double_t val ); template TMatrixDSparse operator+ ( Double_t val ,const TMatrixDSparse &source ); template TMatrixDSparse operator- (const TMatrixDSparse &source1,const TMatrixDSparse &source2); template TMatrixDSparse operator- (const TMatrixDSparse &source1,const TMatrixD &source2); template TMatrixDSparse operator- (const TMatrixD &source1,const TMatrixDSparse &source2); template TMatrixDSparse operator- (const TMatrixDSparse &source , Double_t val ); template TMatrixDSparse operator- ( Double_t val ,const TMatrixDSparse &source ); template TMatrixDSparse operator* (const TMatrixDSparse &source1,const TMatrixDSparse &source2); template TMatrixDSparse operator* (const TMatrixDSparse &source1,const TMatrixD &source2); template TMatrixDSparse operator* (const TMatrixD &source1,const TMatrixDSparse &source2); template TMatrixDSparse operator* ( Double_t val ,const TMatrixDSparse &source ); template TMatrixDSparse operator* (const TMatrixDSparse &source, Double_t val ); template TMatrixDSparse &Add (TMatrixDSparse &target, Double_t scalar,const TMatrixDSparse &source); template TMatrixDSparse &ElementMult (TMatrixDSparse &target,const TMatrixDSparse &source); template TMatrixDSparse &ElementDiv (TMatrixDSparse &target,const TMatrixDSparse &source); template Bool_t AreCompatible(const TMatrixDSparse &m1,const TMatrixDSparse &m2,Int_t verbose);