92 template<
class Element>
102 template<
class Element>
105 Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb,1);
112 template<
class Element>
121 if (row[irowmin] < row_lwb || row[irowmax] > row_upb) {
122 Error(
"TMatrixTSparse",
"Inconsistency between row index and its range");
123 if (row[irowmin] < row_lwb) {
124 Info(
"TMatrixTSparse",
"row index lower bound adjusted to %d",row[irowmin]);
125 row_lwb = row[irowmin];
127 if (row[irowmax] > row_upb) {
128 Info(
"TMatrixTSparse",
"row index upper bound adjusted to %d",row[irowmax]);
129 col_lwb = col[irowmax];
132 if (col[icolmin] < col_lwb || col[icolmax] > col_upb) {
133 Error(
"TMatrixTSparse",
"Inconsistency between column index and its range");
134 if (col[icolmin] < col_lwb) {
135 Info(
"TMatrixTSparse",
"column index lower bound adjusted to %d",col[icolmin]);
136 col_lwb = col[icolmin];
138 if (col[icolmax] > col_upb) {
139 Info(
"TMatrixTSparse",
"column index upper bound adjusted to %d",col[icolmax]);
140 col_upb = col[icolmax];
144 Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb,1,nr);
146 SetMatrixArray(nr,row,col,data);
151 template<
class Element>
164 template<
class Element>
177 template<
class Element>
182 Int_t nr_nonzeros = 0;
197 for (
Int_t i = rowLwb; i <= rowLwb+nrows-1; i++)
198 for (
Int_t j = colLwb; j <= colLwb+ncols-1; j++)
199 if (i==j) nr_nonzeros++;
200 Allocate(nrows,ncols,rowLwb,colLwb,1,nr_nonzeros);
219 Error(
"TMatrixTSparse(EMatrixCreatorOp1)",
"operation %d not yet implemented", op);
227 template<
class Element>
251 Error(
"TMatrixTSparse(EMatrixCreatorOp2)",
"operation %d not yet implemented",op);
259 template<
class Element>
283 Error(
"TMatrixTSparse(EMatrixCreatorOp2)",
"operation %d not yet implemented",op);
291 template<
class Element>
315 Error(
"TMatrixTSparse(EMatrixCreatorOp2)",
"operation %d not yet implemented",op);
324 template<
class Element>
328 if ( (nr_nonzeros > 0 && (no_rows == 0 || no_cols == 0)) ||
329 (no_rows < 0 || no_cols < 0 || nr_nonzeros < 0) )
331 Error(
"Allocate",
"no_rows=%d no_cols=%d non_zeros=%d",no_rows,no_cols,nr_nonzeros);
366 template<
class Element>
374 if (arown >= this->
fNrows || arown < 0) {
375 Error(
"InsertRow",
"row %d out of matrix range",rown);
379 if (acoln >= this->
fNcols || acoln < 0) {
380 Error(
"InsertRow",
"column %d out of matrix range",coln);
384 if (acoln+nr > this->
fNcols || nr < 0) {
385 Error(
"InsertRow",
"row length %d out of range",nr);
398 Int_t lIndex = sIndex-1;
399 Int_t rIndex = sIndex-1;
401 for (index = sIndex; index < eIndex; index++) {
404 if (icol >= acoln+nr)
break;
405 if (icol >= acoln) nslots++;
413 const Int_t ndiff = nr-nslots;
418 for (
Int_t irow = arown+1; irow < this->
fNrows+1; irow++)
423 memmove(
fElements,elements_old,(lIndex+1)*
sizeof(Element));
426 if (nelems_old > 0 && nelems_old-rIndex > 0) {
427 memmove(
fColIndex+rIndex+ndiff,colIndex_old+rIndex,(nelems_old-rIndex)*
sizeof(
Int_t));
428 memmove(
fElements+rIndex+ndiff,elements_old+rIndex,(nelems_old-rIndex)*
sizeof(Element));
432 for (
Int_t i = 0; i < nr; i++) {
438 if (colIndex_old)
delete [] (
Int_t*) colIndex_old;
439 if (elements_old)
delete [] (Element*) elements_old;
449 template<
class Element>
457 if (arown >= this->
fNrows || arown < 0) {
458 Error(
"ExtractRow",
"row %d out of matrix range",rown);
462 if (acoln >= this->
fNcols || acoln < 0) {
463 Error(
"ExtractRow",
"column %d out of matrix range",coln);
467 if (acoln+nr > this->
fNcols || nr < 0) {
468 Error(
"ExtractRow",
"row length %d out of range",nr);
476 memset(v,0,nr*
sizeof(Element));
479 for (
Int_t index = sIndex; index < eIndex; index++) {
480 const Int_t icol = pColIndex[index];
481 if (icol < acoln || icol >= acoln+nr)
continue;
482 v[icol-acoln] = pData[index];
490 template<
class Element>
498 Error(
"AMultBt",
"A and B columns incompatible");
503 Error(
"AMultB",
"this = &a");
508 Error(
"AMultB",
"this = &b");
524 Int_t nr_nonzero_rowa = 0;
527 if (pRowIndexa[irowa] < pRowIndexa[irowa+1])
530 Int_t nr_nonzero_rowb = 0;
533 if (pRowIndexb[irowb] < pRowIndexb[irowb+1])
537 const Int_t nc = nr_nonzero_rowa*nr_nonzero_rowb;
546 pRowIndexc[irowa+1] = pRowIndexc[irowa];
547 if (pRowIndexa[irowa] >= pRowIndexa[irowa+1])
continue;
549 if (pRowIndexb[irowb] >= pRowIndexb[irowb+1])
continue;
550 pRowIndexc[irowa+1]++;
551 pColIndexc[ielem++] = irowb;
564 const Int_t sIndexa = pRowIndexa[irowc];
565 const Int_t eIndexa = pRowIndexa[irowc+1];
567 const Int_t sIndexb = pRowIndexb[icolc];
568 const Int_t eIndexb = pRowIndexb[icolc+1];
570 Int_t indexb = sIndexb;
571 for (
Int_t indexa = sIndexa; indexa < eIndexa && indexb < eIndexb; indexa++) {
572 const Int_t icola = pColIndexa[indexa];
573 while (indexb < eIndexb && pColIndexb[indexb] <= icola) {
574 if (icola == pColIndexb[indexb]) {
575 sum += pDataa[indexa]*pDatab[indexb];
582 pColIndexc[indexc_r] = icolc;
583 pDatac[indexc_r] =
sum;
587 pRowIndexc[irowc+1] = indexc_r;
598 template<
class Element>
606 Error(
"AMultBt",
"A and B columns incompatible");
611 Error(
"AMultB",
"this = &a");
616 Error(
"AMultB",
"this = &b");
630 Int_t nr_nonzero_rowa = 0;
633 if (pRowIndexa[irowa] < pRowIndexa[irowa+1])
638 const Int_t nc = nr_nonzero_rowa*nr_nonzero_rowb;
647 pRowIndexc[irowa+1] = pRowIndexc[irowa];
649 pRowIndexc[irowa+1]++;
650 pColIndexc[ielem++] = irowb;
663 const Int_t sIndexa = pRowIndexa[irowc];
664 const Int_t eIndexa = pRowIndexa[irowc+1];
668 for (
Int_t indexa = sIndexa; indexa < eIndexa; indexa++) {
669 const Int_t icola = pColIndexa[indexa];
670 sum += pDataa[indexa]*pDatab[off+icola];
673 pColIndexc[indexc_r] = icolc;
674 pDatac[indexc_r] =
sum;
678 pRowIndexc[irowc+1]= indexc_r;
689 template<
class Element>
697 Error(
"AMultBt",
"A and B columns incompatible");
702 Error(
"AMultB",
"this = &a");
707 Error(
"AMultB",
"this = &b");
722 Int_t nr_nonzero_rowb = 0;
725 if (pRowIndexb[irowb] < pRowIndexb[irowb+1])
729 const Int_t nc = nr_nonzero_rowa*nr_nonzero_rowb;
738 pRowIndexc[irowa+1] = pRowIndexc[irowa];
740 if (pRowIndexb[irowb] >= pRowIndexb[irowb+1])
continue;
741 pRowIndexc[irowa+1]++;
742 pColIndexc[ielem++] = irowb;
757 const Int_t sIndexb = pRowIndexb[icolc];
758 const Int_t eIndexb = pRowIndexb[icolc+1];
760 for (
Int_t indexb = sIndexb; indexb < eIndexb; indexb++) {
761 const Int_t icolb = pColIndexb[indexb];
762 sum += pDataa[off+icolb]*pDatab[indexb];
765 pColIndexc[indexc_r] = icolc;
766 pDatac[indexc_r] =
sum;
770 pRowIndexc[irowc+1] = indexc_r;
781 template<
class Element>
790 Error(
"APlusB(const TMatrixTSparse &,const TMatrixTSparse &",
"matrices not compatible");
795 Error(
"APlusB",
"this = &a");
800 Error(
"APlusB",
"this = &b");
823 const Int_t sIndexa = pRowIndexa[irowc];
824 const Int_t eIndexa = pRowIndexa[irowc+1];
825 const Int_t sIndexb = pRowIndexb[irowc];
826 const Int_t eIndexb = pRowIndexb[irowc+1];
827 Int_t indexa = sIndexa;
828 Int_t indexb = sIndexb;
831 while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
832 if (icolc == pColIndexa[indexa]) {
833 sum += pDataa[indexa];
838 while (indexb < eIndexb && pColIndexb[indexb] <= icolc) {
839 if (icolc == pColIndexb[indexb]) {
840 sum += pDatab[indexb];
847 pColIndexc[indexc_r] = icolc;
848 pDatac[indexc_r] =
sum;
852 pRowIndexc[irowc+1] = indexc_r;
863 template<
class Element>
872 Error(
"APlusB(const TMatrixTSparse &,const TMatrixT &",
"matrices not compatible");
877 Error(
"APlusB",
"this = &a");
882 Error(
"APlusB",
"this = &b");
903 const Int_t sIndexa = pRowIndexa[irowc];
904 const Int_t eIndexa = pRowIndexa[irowc+1];
906 Int_t indexa = sIndexa;
908 Element
sum = pDatab[off+icolc];
909 while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
910 if (icolc == pColIndexa[indexa]) {
911 sum += pDataa[indexa];
918 pColIndexc[indexc_r] = icolc;
919 pDatac[indexc_r] =
sum;
923 pRowIndexc[irowc+1] = indexc_r;
934 template<
class Element>
943 Error(
"AMinusB(const TMatrixTSparse &,const TMatrixTSparse &",
"matrices not compatible");
948 Error(
"AMinusB",
"this = &a");
953 Error(
"AMinusB",
"this = &b");
976 const Int_t sIndexa = pRowIndexa[irowc];
977 const Int_t eIndexa = pRowIndexa[irowc+1];
978 const Int_t sIndexb = pRowIndexb[irowc];
979 const Int_t eIndexb = pRowIndexb[irowc+1];
980 Int_t indexa = sIndexa;
981 Int_t indexb = sIndexb;
984 while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
985 if (icolc == pColIndexa[indexa]) {
986 sum += pDataa[indexa];
991 while (indexb < eIndexb && pColIndexb[indexb] <= icolc) {
992 if (icolc == pColIndexb[indexb]) {
993 sum -= pDatab[indexb];
1000 pColIndexc[indexc_r] = icolc;
1001 pDatac[indexc_r] =
sum;
1005 pRowIndexc[irowc+1] = indexc_r;
1016 template<
class Element>
1025 Error(
"AMinusB(const TMatrixTSparse &,const TMatrixT &",
"matrices not compatible");
1030 Error(
"AMinusB",
"this = &a");
1035 Error(
"AMinusB",
"this = &b");
1056 const Int_t sIndexa = pRowIndexa[irowc];
1057 const Int_t eIndexa = pRowIndexa[irowc+1];
1059 Int_t indexa = sIndexa;
1061 Element
sum = -pDatab[off+icolc];
1062 while (indexa < eIndexa && pColIndexa[indexa] <= icolc) {
1063 if (icolc == pColIndexa[indexa]) {
1064 sum += pDataa[indexa];
1071 pColIndexc[indexc_r] = icolc;
1072 pDatac[indexc_r] =
sum;
1076 pRowIndexc[irowc+1] = indexc_r;
1087 template<
class Element>
1096 Error(
"AMinusB(const TMatrixT &,const TMatrixTSparse &",
"matrices not compatible");
1101 Error(
"AMinusB",
"this = &a");
1106 Error(
"AMinusB",
"this = &b");
1127 const Int_t sIndexb = pRowIndexb[irowc];
1128 const Int_t eIndexb = pRowIndexb[irowc+1];
1130 Int_t indexb = sIndexb;
1132 Element
sum = pDataa[off+icolc];
1133 while (indexb < eIndexb && pColIndexb[indexb] <= icolc) {
1134 if (icolc == pColIndexb[indexb]) {
1135 sum -= pDatab[indexb];
1142 pColIndexc[indexc_r] = icolc;
1143 pDatac[indexc_r] =
sum;
1147 pRowIndexc[irowc+1] = indexc_r;
1157 template<
class Element>
1163 memcpy(data,elem,this->
fNelems*
sizeof(Element));
1171 template<
class Element>
1176 Error(
"SetMatrixArray(Int_t,Int_t*,Int_t*,Element*",
"nr <= 0");
1189 Error(
"SetMatrixArray",
"Inconsistency between row index and its range");
1190 if (row[irowmin] < this->
fRowLwb) {
1191 Info(
"SetMatrixArray",
"row index lower bound adjusted to %d",row[irowmin]);
1195 Info(
"SetMatrixArray",
"row index upper bound adjusted to %d",row[irowmax]);
1200 Error(
"SetMatrixArray",
"Inconsistency between column index and its range");
1201 if (col[icolmin] < this->
fColLwb) {
1202 Info(
"SetMatrixArray",
"column index lower bound adjusted to %d",col[icolmin]);
1206 Info(
"SetMatrixArray",
"column index upper bound adjusted to %d",col[icolmax]);
1213 Int_t nr_nonzeros = 0;
1214 const Element *ep =
data;
1215 const Element *
const fp = data+nr;
1218 if (*ep++ != 0.0) nr_nonzeros++;
1220 if (nr_nonzeros != this->
fNelems) {
1239 for (
Int_t irow = 1; irow < this->
fNrows+1; irow++) {
1240 if (ielem < nr && row[ielem] < irow) {
1241 while (ielem < nr) {
1242 if (data[ielem] != 0.0) {
1248 if (ielem >= nr || row[ielem] != row[ielem-1])
1261 template<
class Element>
1264 if (nelems_new != this->
fNelems) {
1275 memmove(
fElements,oDp,nr*
sizeof(Element));
1279 if (nelems_new > nr) {
1280 memset(
fElements+nr,0,(nelems_new-nr)*
sizeof(Element));
1295 template<
class Element>
1302 Error(
"SetSparseIndex",
"matrices not compatible");
1309 if (nr_nonzeros != this->
fNelems)
1318 for (
Int_t irow = 0; irow < this->
fNrows; irow++) {
1320 for (
Int_t icol = 0; icol < this->
fNcols; icol++) {
1338 template<
class Element>
1347 Error(
"SetSparseIndexAB",
"source matrices not compatible");
1353 Error(
"SetSparseIndexAB",
"matrix not compatible with source matrices");
1366 const Int_t sIndexa = pRowIndexa[irowc];
1367 const Int_t eIndexa = pRowIndexa[irowc+1];
1368 const Int_t sIndexb = pRowIndexb[irowc];
1369 const Int_t eIndexb = pRowIndexb[irowc+1];
1370 nc += eIndexa-sIndexa;
1371 Int_t indexb = sIndexb;
1372 for (
Int_t indexa = sIndexa; indexa < eIndexa; indexa++) {
1373 const Int_t icola = pColIndexa[indexa];
1374 for (; indexb < eIndexb; indexb++) {
1375 if (pColIndexb[indexb] >= icola) {
1376 if (pColIndexb[indexb] == icola)
1383 while (indexb < eIndexb) {
1384 const Int_t icola = (eIndexa > sIndexa && eIndexa > 0) ? pColIndexa[eIndexa-1] : -1;
1385 if (pColIndexb[indexb++] > icola)
1400 const Int_t sIndexa = pRowIndexa[irowc];
1401 const Int_t eIndexa = pRowIndexa[irowc+1];
1402 const Int_t sIndexb = pRowIndexb[irowc];
1403 const Int_t eIndexb = pRowIndexb[irowc+1];
1404 Int_t indexb = sIndexb;
1405 for (
Int_t indexa = sIndexa; indexa < eIndexa; indexa++) {
1406 const Int_t icola = pColIndexa[indexa];
1407 for (; indexb < eIndexb; indexb++) {
1408 if (pColIndexb[indexb] >= icola) {
1409 if (pColIndexb[indexb] == icola)
1413 pColIndexc[nc++] = pColIndexb[indexb];
1415 pColIndexc[nc++] = pColIndexa[indexa];
1417 while (indexb < eIndexb) {
1418 const Int_t icola = (eIndexa > 0) ? pColIndexa[eIndexa-1] : -1;
1419 if (pColIndexb[indexb++] > icola)
1420 pColIndexc[nc++] = pColIndexb[indexb-1];
1422 pRowIndexc[irowc+1] = nc;
1432 template<
class Element>
1441 Error(
"SetSparseIndexAB",
"source matrices not compatible");
1447 Error(
"SetSparseIndexAB",
"matrix not compatible with source matrices");
1460 const Int_t sIndexb = pRowIndexb[irowc];
1461 const Int_t eIndexb = pRowIndexb[irowc+1];
1463 Int_t indexb = sIndexb;
1465 if (pDataa[off+icolc] != 0.0 || pColIndexb[indexb] > icolc)
continue;
1466 for (; indexb < eIndexb; indexb++) {
1467 if (pColIndexb[indexb] >= icolc) {
1468 if (pColIndexb[indexb] == icolc) {
1488 const Int_t sIndexb = pRowIndexb[irowc];
1489 const Int_t eIndexb = pRowIndexb[irowc+1];
1491 Int_t indexb = sIndexb;
1493 if (pDataa[off+icolc] != 0.0)
1494 pColIndexc[nc++] = icolc;
1495 else if (pColIndexb[indexb] <= icolc) {
1496 for (; indexb < eIndexb; indexb++) {
1497 if (pColIndexb[indexb] >= icolc) {
1498 if (pColIndexb[indexb] == icolc)
1499 pColIndexc[nc++] = pColIndexb[indexb++];
1505 pRowIndexc[irowc+1] = nc;
1517 template<
class Element>
1522 Error(
"ResizeTo(Int_t,Int_t,Int_t)",
"Not owner of data array,cannot resize");
1528 (this->
fNelems == nr_nonzeros || nr_nonzeros < 0))
1530 else if (nrows == 0 || ncols == 0 || nr_nonzeros == 0) {
1544 if (nr_nonzeros > 0)
1545 nelems_new = nr_nonzeros;
1548 for (
Int_t irow = 0; irow < nrows_old; irow++) {
1549 if (irow >= nrows)
continue;
1550 const Int_t sIndex = rowIndex_old[irow];
1551 const Int_t eIndex = rowIndex_old[irow+1];
1552 for (
Int_t index = sIndex; index < eIndex; index++) {
1553 const Int_t icol = colIndex_old[index];
1560 Allocate(nrows,ncols,0,0,1,nelems_new);
1567 Int_t nelems_copy = 0;
1568 rowIndex_new[0] = 0;
1570 for (
Int_t irow = 0; irow < nrows_old && cont; irow++) {
1571 if (irow >= nrows)
continue;
1572 const Int_t sIndex = rowIndex_old[irow];
1573 const Int_t eIndex = rowIndex_old[irow+1];
1574 for (
Int_t index = sIndex; index < eIndex; index++) {
1575 const Int_t icol = colIndex_old[index];
1577 rowIndex_new[irow+1] = nelems_copy+1;
1578 colIndex_new[nelems_copy] = icol;
1579 elements_new[nelems_copy] = elements_old[index];
1582 if (nelems_copy >= nelems_new) {
1589 if (rowIndex_old)
delete [] (
Int_t*) rowIndex_old;
1590 if (colIndex_old)
delete [] (
Int_t*) colIndex_old;
1591 if (elements_old)
delete [] (Element*) elements_old;
1595 rowIndex_new[irow] = rowIndex_new[nrowIndex_old-1];
1598 const Int_t nelems_new = (nr_nonzeros >= 0) ? nr_nonzeros : 0;
1599 Allocate(nrows,ncols,0,0,1,nelems_new);
1611 template<
class Element>
1617 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t,Int_t)",
"Not owner of data array,cannot resize");
1621 const Int_t new_nrows = row_upb-row_lwb+1;
1622 const Int_t new_ncols = col_upb-col_lwb+1;
1625 if (this->
fNrows == new_nrows && this->
fNcols == new_ncols &&
1627 (this->
fNelems == nr_nonzeros || nr_nonzeros < 0))
1629 else if (new_nrows == 0 || new_ncols == 0 || nr_nonzeros == 0) {
1647 if (nr_nonzeros > 0)
1648 nelems_new = nr_nonzeros;
1651 for (
Int_t irow = 0; irow < nrows_old; irow++) {
1652 if (irow+rowLwb_old > row_upb || irow+rowLwb_old < row_lwb)
continue;
1653 const Int_t sIndex = rowIndex_old[irow];
1654 const Int_t eIndex = rowIndex_old[irow+1];
1655 for (
Int_t index = sIndex; index < eIndex; index++) {
1656 const Int_t icol = colIndex_old[index];
1657 if (icol+colLwb_old <= col_upb && icol+colLwb_old >= col_lwb)
1663 Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1,nelems_new);
1670 Int_t nelems_copy = 0;
1671 rowIndex_new[0] = 0;
1672 const Int_t row_off = rowLwb_old-row_lwb;
1673 const Int_t col_off = colLwb_old-col_lwb;
1674 for (
Int_t irow = 0; irow < nrows_old; irow++) {
1675 if (irow+rowLwb_old > row_upb || irow+rowLwb_old < row_lwb)
continue;
1676 const Int_t sIndex = rowIndex_old[irow];
1677 const Int_t eIndex = rowIndex_old[irow+1];
1678 for (
Int_t index = sIndex; index < eIndex; index++) {
1679 const Int_t icol = colIndex_old[index];
1680 if (icol+colLwb_old <= col_upb && icol+colLwb_old >= col_lwb) {
1681 rowIndex_new[irow+row_off+1] = nelems_copy+1;
1682 colIndex_new[nelems_copy] = icol+col_off;
1683 elements_new[nelems_copy] = elements_old[index];
1686 if (nelems_copy >= nelems_new) {
1692 if (rowIndex_old)
delete [] (
Int_t*) rowIndex_old;
1693 if (colIndex_old)
delete [] (
Int_t*) colIndex_old;
1694 if (elements_old)
delete [] (Element*) elements_old;
1698 rowIndex_new[irow] = rowIndex_new[nrowIndex_old-1];
1701 const Int_t nelems_new = (nr_nonzeros >= 0) ? nr_nonzeros : 0;
1702 Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1,nelems_new);
1710 template<
class Element>
1715 if (row_upb < row_lwb) {
1716 Error(
"Use",
"row_upb=%d < row_lwb=%d",row_upb,row_lwb);
1719 if (col_upb < col_lwb) {
1720 Error(
"Use",
"col_upb=%d < col_lwb=%d",col_upb,col_lwb);
1727 this->
fNrows = row_upb-row_lwb+1;
1728 this->
fNcols = col_upb-col_lwb+1;
1750 template<
class Element>
1757 Error(
"GetSub",
"row_lwb out-of-bounds");
1761 Error(
"GetSub",
"col_lwb out-of-bounds");
1765 Error(
"GetSub",
"row_upb out-of-bounds");
1769 Error(
"GetSub",
"col_upb out-of-bounds");
1772 if (row_upb < row_lwb || col_upb < col_lwb) {
1773 Error(
"GetSub",
"row_upb < row_lwb || col_upb < col_lwb");
1782 const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
1783 const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
1784 const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
1785 const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
1787 Int_t nr_nonzeros = 0;
1788 for (
Int_t irow = 0; irow < this->
fNrows; irow++) {
1789 if (irow+this->
fRowLwb > row_upb || irow+this->
fRowLwb < row_lwb)
continue;
1792 for (
Int_t index = sIndex; index < eIndex; index++) {
1794 if (icol+this->fColLwb <= col_upb && icol+this->
fColLwb >= col_lwb)
1799 target.
ResizeTo(row_lwb_sub,row_upb_sub,col_lwb_sub,col_upb_sub,nr_nonzeros);
1808 Int_t nelems_copy = 0;
1809 rowIndex_sub[0] = 0;
1812 for (
Int_t irow = 0; irow < this->
fNrows; irow++) {
1813 if (irow+this->
fRowLwb > row_upb || irow+this->
fRowLwb < row_lwb)
continue;
1816 for (
Int_t index = sIndex; index < eIndex; index++) {
1818 if (icol+this->fColLwb <= col_upb && icol+this->
fColLwb >= col_lwb) {
1819 rowIndex_sub[irow+row_off+1] = nelems_copy+1;
1820 colIndex_sub[nelems_copy] = icol+col_off;
1821 ep_sub[nelems_copy] = ep[index];
1829 const Int_t ncols_sub = col_upb_sub-col_lwb_sub+1;
1830 for (
Int_t irow = 0; irow < this->
fNrows; irow++) {
1831 if (irow+this->
fRowLwb > row_upb || irow+this->
fRowLwb < row_lwb)
continue;
1834 const Int_t off = (irow+row_off)*ncols_sub;
1835 for (
Int_t index = sIndex; index < eIndex; index++) {
1837 if (icol+this->fColLwb <= col_upb && icol+this->
fColLwb >= col_lwb)
1838 ep_sub[off+icol+col_off] = ep[index];
1850 template<
class Element>
1858 Error(
"SetSub",
"row_lwb out-of-bounds");
1862 Error(
"SetSub",
"col_lwb out-of-bounds");
1866 Error(
"SetSub",
"source matrix too large");
1876 Int_t nr_nonzeros = 0;
1877 for (
Int_t irow = 0; irow < this->
fNrows; irow++) {
1878 if (irow+this->
fRowLwb >= row_lwb+nRows_source || irow+this->
fRowLwb < row_lwb)
continue;
1881 for (
Int_t index = sIndex; index < eIndex; index++) {
1883 if (icol+this->fColLwb < col_lwb+nCols_source && icol+this->
fColLwb >= col_lwb)
1897 const Int_t nelems_new = nelems_old+source.
NonZeros()-nr_nonzeros;
1911 rowIndex_new[0] = 0;
1912 for (
Int_t irow = 0; irow < this->
fNrows; irow++) {
1913 rowIndex_new[irow+1] = rowIndex_new[irow];
1915 if (irow+this->fRowLwb < row_lwb+nRows_source && irow+this->fRowLwb >= row_lwb)
1918 const Int_t sIndex_o = rowIndex_old[irow];
1919 const Int_t eIndex_o = rowIndex_old[irow+1];
1922 const Int_t icol_left = col_off-1;
1924 for (
Int_t index = sIndex_o; index <= left; index++) {
1925 rowIndex_new[irow+1]++;
1926 colIndex_new[nr] = colIndex_old[index];
1927 elements_new[nr] = elements_old[index];
1931 if (rowIndex_s && colIndex_s) {
1932 const Int_t sIndex_s = rowIndex_s[irow-row_off];
1933 const Int_t eIndex_s = rowIndex_s[irow-row_off+1];
1934 for (
Int_t index = sIndex_s; index < eIndex_s; index++) {
1935 rowIndex_new[irow+1]++;
1936 colIndex_new[nr] = colIndex_s[index]+col_off;
1937 elements_new[nr] = elements_s[index];
1941 const Int_t off = (irow-row_off)*nCols_source;
1942 for (
Int_t icol = 0; icol < nCols_source; icol++) {
1943 rowIndex_new[irow+1]++;
1944 colIndex_new[nr] = icol+col_off;
1945 elements_new[nr] = elements_s[off+icol];
1950 const Int_t icol_right = col_off+nCols_source-1;
1953 while (right < eIndex_o-1 && colIndex_old[right+1] <= icol_right)
1957 for (
Int_t index = right; index < eIndex_o; index++) {
1958 rowIndex_new[irow+1]++;
1959 colIndex_new[nr] = colIndex_old[index];
1960 elements_new[nr] = elements_old[index];
1965 for (
Int_t index = sIndex_o; index < eIndex_o; index++) {
1966 rowIndex_new[irow+1]++;
1967 colIndex_new[nr] = colIndex_old[index];
1968 elements_new[nr] = elements_old[index];
1976 if (rowIndex_old)
delete [] (
Int_t*) rowIndex_old;
1977 if (colIndex_old)
delete [] (
Int_t*) colIndex_old;
1978 if (elements_old)
delete [] (Element*) elements_old;
1986 template<
class Element>
1995 Error(
"Transpose",
"matrix has wrong shape");
2011 Element *pData_t =
new Element[nr_nonzeros];
2014 for (
Int_t irow_s = 0; irow_s < source.
GetNrows(); irow_s++) {
2015 const Int_t sIndex = pRowIndex_s[irow_s];
2016 const Int_t eIndex = pRowIndex_s[irow_s+1];
2017 for (
Int_t index = sIndex; index < eIndex; index++) {
2018 rownr[ielem] = pColIndex_s[index]+this->
fRowLwb;
2019 colnr[ielem] = irow_s+this->
fColLwb;
2020 pData_t[ielem] = pData_s[index];
2032 if (pData_t)
delete [] pData_t;
2033 if (rownr)
delete [] rownr;
2034 if (colnr)
delete [] colnr;
2041 template<
class Element>
2057 template<
class Element>
2064 Int_t nr_nonzeros = 0;
2067 if (i == j) nr_nonzeros++;
2069 if (nr_nonzeros != this->
fNelems) {
2073 if (oIp)
delete [] oIp;
2076 if (oDp)
delete [] oDp;
2081 for (i = this->
fRowLwb; i <= this->
fRowLwb+this->fNrows-1; i++) {
2099 template<
class Element>
2105 const Element *
const fp = ep+this->
fNelems;
2110 for (
Int_t irow = 0; irow < this->
fNrows; irow++) {
2111 const Int_t sIndex = pR[irow];
2112 const Int_t eIndex = pR[irow+1];
2114 for (
Int_t index = sIndex; index < eIndex; index++)
2128 template<
class Element>
2136 const Element *
const fp = ep+this->
fNcols;
2156 template<
class Element>
2163 if (arown >= this->
fNrows || arown < 0) {
2164 Error(
"operator()",
"Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->
fNrows);
2167 if (acoln >= this->
fNcols || acoln < 0) {
2168 Error(
"operator()",
"Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->
fNcols);
2181 if (index >= sIndex &&
fColIndex[index] == acoln)
2189 if (index >= sIndex &&
fColIndex[index] == acoln)
2192 Error(
"operator()(Int_t,Int_t",
"Insert row failed");
2200 template <
class Element>
2205 Error(
"operator()(Int_t,Int_t) const",
"row/col indices are not set");
2212 if (arown >= this->
fNrows || arown < 0) {
2213 Error(
"operator()",
"Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->
fNrows);
2216 if (acoln >= this->
fNcols || acoln < 0) {
2217 Error(
"operator()",
"Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->
fNcols);
2224 if (index < sIndex ||
fColIndex[index] != acoln)
return 0.0;
2232 template<
class Element>
2236 Error(
"operator=(const TMatrixTSparse &)",
"matrices not compatible");
2245 memcpy(tp,sp,this->
fNelems*
sizeof(Element));
2255 template<
class Element>
2259 Error(
"operator=(const TMatrixT &)",
"matrices not compatible");
2268 for (
Int_t irow = 0; irow < this->
fNrows; irow++) {
2272 for (
Int_t index = sIndex; index < eIndex; index++) {
2274 tp[index] = sp[off+icol];
2286 template<
class Element>
2292 Error(
"operator=(Element",
"row/col indices are not set");
2297 const Element *
const ep_last = ep+this->
fNelems;
2298 while (ep < ep_last)
2307 template<
class Element>
2313 const Element *
const ep_last = ep+this->
fNelems;
2314 while (ep < ep_last)
2323 template<
class Element>
2329 const Element *
const ep_last = ep+this->
fNelems;
2330 while (ep < ep_last)
2339 template<
class Element>
2345 const Element *
const ep_last = ep+this->
fNelems;
2346 while (ep < ep_last)
2355 template<
class Element>
2360 const Element scale = beta-alpha;
2361 const Element shift = alpha/scale;
2376 for (
Int_t k = 0; k <
nn; k++) {
2377 const Element
r =
Drand(seed);
2379 if ((nn-k)*r < length-chosen) {
2380 pColIndex[chosen] = k%
n;
2383 if (irow > icurrent) {
2384 for ( ; icurrent < irow; icurrent++)
2385 pRowIndex[icurrent+1] = chosen;
2387 ep[chosen] = scale*(
Drand(seed)+shift);
2391 for ( ; icurrent <
m; icurrent++)
2392 pRowIndex[icurrent+1] = length;
2402 template<
class Element>
2405 const Element scale = beta-alpha;
2406 const Element shift = alpha/scale;
2412 Error(
"RandomizePD(Element &",
"matrix should be square");
2428 ep[0] = 1
e-8+scale*(
Drand(seed)+shift);
2438 length = (length <=
nn) ? length : nn;
2448 for (
Int_t k = 0; k <
nn; k++ ) {
2449 const Element
r =
Drand(seed);
2451 if( (nn-k)*r < length-chosen) {
2458 Int_t col = k-row*(row+1)/2;
2463 if (row > icurrent) {
2467 for ( ; icurrent < row; icurrent++) {
2470 for (
Int_t ll = pRowIndex[icurrent]; ll < nnz; ll++)
2472 ep[nnz] += 1
e-8+scale*(
Drand(seed)+shift);
2473 pColIndex[nnz] = icurrent;
2476 pRowIndex[icurrent+1] = nnz;
2479 ep[nnz] = scale*(
Drand(seed)+shift);
2480 pColIndex[nnz] = col;
2483 ep[pRowIndex[col+1]-1] +=
TMath::Abs(ep[nnz]);
2493 for ( ; icurrent <
n; icurrent++) {
2496 for(
Int_t ll = pRowIndex[icurrent]; ll < nnz; ll++)
2498 ep[nnz] += 1
e-8+scale*(
Drand(seed)+shift);
2499 pColIndex[nnz] = icurrent;
2502 pRowIndex[icurrent+1] = nnz;
2515 for (
Int_t irow = 0; irow < this->
fNrows+1; irow++) {
2516 const Int_t sIndex = pR[irow];
2517 const Int_t eIndex = pR[irow+1];
2518 for (
Int_t index = sIndex; index < eIndex; index++) {
2519 const Int_t icol = pC[index];
2531 template<
class Element>
2540 template<
class Element>
2549 template<
class Element>
2558 template<
class Element>
2568 template<
class Element>
2578 template<
class Element>
2587 template<
class Element>
2596 template<
class Element>
2605 template<
class Element>
2615 template<
class Element>
2625 template<
class Element>
2634 template<
class Element>
2643 template<
class Element>
2652 template<
class Element>
2662 template<
class Element>
2673 template<
class Element>
2676 target += scalar * source;
2683 template<
class Element>
2687 ::Error(
"ElementMult(TMatrixTSparse &,const TMatrixTSparse &)",
"matrices not compatible");
2703 template<
class Element>
2707 ::Error(
"ElementDiv(TMatrixT &,const TMatrixT &)",
"matrices not compatible");
2714 while ( tp < ftp ) {
2718 Error(
"ElementDiv",
"source element is zero");
2728 template<
class Element>
2733 ::Error(
"AreCompatible",
"matrix 1 not valid");
2738 ::Error(
"AreCompatible",
"matrix 2 not valid");
2745 ::Error(
"AreCompatible",
"matrices 1 and 2 not compatible");
2752 if (memcmp(pR1,pR2,(nRows+1)*
sizeof(
Int_t))) {
2754 ::Error(
"AreCompatible",
"matrices 1 and 2 have different rowIndex");
2755 for (
Int_t i = 0; i < nRows+1; i++)
2756 printf(
"%d: %d %d\n",i,pR1[i],pR2[i]);
2762 if (memcmp(pD1,pD2,nData*
sizeof(
Int_t))) {
2764 ::Error(
"AreCompatible",
"matrices 1 and 2 have different colIndex");
2765 for (
Int_t i = 0; i < nData; i++)
2766 printf(
"%d: %d %d\n",i,pD1[i],pD2[i]);
2776 template<
class Element>
2793 #ifndef ROOT_TMatrixFSparsefwd 2796 #ifndef ROOT_TMatrixFfwd 2822 #ifndef ROOT_TMatrixDSparsefwd 2825 #ifndef ROOT_TMatrixDfwd void AMinusB(const TMatrixTSparse< Element > &a, const TMatrixTSparse< Element > &b, Int_t constr=0)
General matrix subtraction.
virtual const Element * GetMatrixArray() const =0
virtual const Element * GetMatrixArray() const
virtual Int_t NonZeros() const
Compute the number of elements != 0.0.
virtual Int_t NonZeros() const
Compute the number of elements != 0.0.
static long int sum(long int i)
virtual TMatrixTBase< Element > & InsertRow(Int_t row, Int_t col, const Element *v, Int_t n=-1)
Insert in row rown, n elements of array v at column coln.
virtual void GetMatrix2Array(Element *data, Option_t *="") const
Copy matrix data to array . It is assumed that array is of size >= fNelems.
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Double_t Floor(Double_t x)
Long64_t LocMax(Long64_t n, const T *a)
virtual const Int_t * GetRowIndexArray() const
TMatrixTSparse< Element > & Transpose(const TMatrixTSparse< Element > &source)
Transpose a matrix.
virtual const Element * GetMatrixArray() const
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
template TMatrixDSparse & Add< Double_t >(TMatrixDSparse &target, Double_t scalar, const TMatrixDSparse &source)
void ToUpper()
Change string to upper case.
Buffer base class used for serializing objects.
template TMatrixDSparse & ElementMult< Double_t >(TMatrixDSparse &target, const TMatrixDSparse &source)
template TMatrixFSparse & ElementDiv< Float_t >(TMatrixFSparse &target, const TMatrixFSparse &source)
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *="")
Copy array data to matrix .
Short_t Min(Short_t a, Short_t b)
Element operator()(Int_t rown, Int_t coln) const
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1)=0
template TMatrixFSparse & Add< Float_t >(TMatrixFSparse &target, Float_t scalar, const TMatrixFSparse &source)
double beta(double x, double y)
Calculates the beta function.
template TMatrixFSparse & ElementMult< Float_t >(TMatrixFSparse &target, const TMatrixFSparse &source)
Int_t GetNoElements() const
static void DoubleLexSort(Int_t n, Int_t *first, Int_t *second, Element *data)
default kTRUE, when Use array kFALSE
#define templateClassImp(name)
template Bool_t AreCompatible< Float_t >(const TMatrixFSparse &m1, const TMatrixFSparse &m2, Int_t verbose)
void Info(const char *location, const char *msgfmt,...)
Bool_t AreCompatible(const TMatrixTSparse< Element > &m1, const TMatrixTSparse< Element > &m2, Int_t verbose)
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
randomize matrix element values
void AMultBt(const TMatrixTSparse< Element > &a, const TMatrixTSparse< Element > &b, Int_t constr=0)
General matrix multiplication.
TMatrixTSparse< Element > & ElementMult(TMatrixTSparse< Element > &target, const TMatrixTSparse< Element > &source)
Multiply target by the source, element-by-element.
TObject & operator=(const TObject &rhs)
TObject assignment operator.
virtual TMatrixTSparse< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
TMatrixTSparse< Element > & Add(TMatrixTSparse< Element > &target, Element scalar, const TMatrixTSparse< Element > &source)
Modify addition: target += scalar * source.
void Error(const char *location, const char *msgfmt,...)
virtual Element ColNorm() const
Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
TMatrixTSparse< Element > & operator-=(Element val)
Subtract val from every element of the matrix.
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
template TMatrixDSparse & ElementDiv< Double_t >(TMatrixDSparse &target, const TMatrixDSparse &source)
virtual void Clear(Option_t *="")
TMatrixTSparse< Element > & operator*=(Element val)
Multiply every element of the matrix with val.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
TMatrixTSparse< Element > & SetSparseIndex(Int_t nelem_new)
Increase/decrease the number of non-zero elements to nelems_new.
TMatrixTSparse< Element > operator+(const TMatrixTSparse< Element > &source1, const TMatrixTSparse< Element > &source2)
R__EXTERN Int_t gMatrixCheck
virtual const Int_t * GetColIndexArray() const
void Allocate(Int_t nrows, Int_t ncols, Int_t row_lwb=0, Int_t col_lwb=0, Int_t init=0, Int_t nr_nonzeros=0)
Allocate new matrix.
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
TMatrixTSparse< Element > & ElementDiv(TMatrixTSparse< Element > &target, const TMatrixTSparse< Element > &source)
Divide target by the source, element-by-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)
TMatrixTSparse< Element > & SetSparseIndexAB(const TMatrixTSparse< Element > &a, const TMatrixTSparse< Element > &b)
Set the row/column indices to the "sum" of matrices a and b It is checked that enough space has been ...
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
template Bool_t AreCompatible< Double_t >(const TMatrixDSparse &m1, const TMatrixDSparse &m2, Int_t verbose)
TCppObject_t Allocate(TCppType_t type)
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1)
Set size of the matrix to nrows x ncols with nr_nonzeros non-zero entries if nr_nonzeros > 0 ...
TMatrixTSparse< Element > & operator+=(Element val)
Add val to every element of the matrix.
Short_t Max(Short_t a, Short_t b)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
virtual Element RowNorm() const
Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
TMatrixTSparse< Element > & operator=(const TMatrixT< Element > &source)
Notice that the sparsity of the matrix is NOT changed : its fRowIndex/fColIndex are used ! ...
Double_t Sqrt(Double_t x)
void AMultB(const TMatrixTSparse< Element > &a, const TMatrixTSparse< Element > &b, Int_t constr=0)
virtual TMatrixTBase< Element > & GetSub(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, TMatrixTBase< Element > &target, Option_t *option="S") const
Get submatrix [row_lwb..row_upb][col_lwb..col_upb]; The indexing range of the returned matrix depends...
virtual void ExtractRow(Int_t row, Int_t col, Element *v, Int_t n=-1) const
Store in array v, n matrix elements of row rown starting at column coln.
Long64_t LocMin(Long64_t n, const T *a)
TMatrixTSparse< Element > operator-(const TMatrixTSparse< Element > &source1, const TMatrixTSparse< Element > &source2)
void APlusB(const TMatrixTSparse< Element > &a, const TMatrixTSparse< Element > &b, Int_t constr=0)
General matrix addition.
double norm(double *x, double *p)
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
TMatrixTSparse< Element > operator*(const TMatrixTSparse< Element > &source1, const TMatrixTSparse< Element > &source2)
Long64_t BinarySearch(Long64_t n, const T *array, T value)
virtual TMatrixTBase< Element > & SetSub(Int_t row_lwb, Int_t col_lwb, const TMatrixTBase< Element > &source)
Insert matrix source starting at [row_lwb][col_lwb], thereby overwriting the part [row_lwb...
virtual const Int_t * GetColIndexArray() const =0
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
virtual const Int_t * GetRowIndexArray() const =0