34 memset(fInfo,0,21*
sizeof(
Int_t));
74 fNrows = row_upb-row_lwb+1;
123 Int_t nr_nonzeros = 0;
124 for (
Int_t irow = 0; irow < nrows; irow++ ) {
125 const Int_t rown = irow+rowLwb;
126 for (
Int_t index = pRowIndex[irow]; index < pRowIndex[irow+1]; index++ ) {
127 const Int_t coln = pColIndex[index]+colLwb;
128 if (coln >= rown) nr_nonzeros++;
149 for (
Int_t irow = 0; irow < nrows; irow++ ) {
150 const Int_t rown = irow+rowLwb;
151 for (
Int_t index = pRowIndex[irow]; index < pRowIndex[irow+1]; index++ ) {
152 const Int_t coln = pColIndex[index]+colLwb;
153 if (coln >= rown) b[nr++] = pData[index];
165 fA.
Use(*const_cast<TMatrixDSparse *>(&a));
180 for (
Int_t index = rowIndex[irow]; index < rowIndex[irow+1]; index++ ) {
203 Error(
"SetMatrix(const TMatrixDSparse &",
"nRows = %d out of range",fNrows);
206 Error(
"SetMatrix(const TMatrixDSparse &",
"nr_nonzeros = %d out of range",
fNnonZeros);
209 Error(
"SetMatrix(const TMatrixDSparse &",
210 "insufficient space in fIw of %d suggest reset to %d",
fIw.
GetSize(),this->
IError());
213 Error(
"SetMatrix(const TMatrixDSparse &",
214 "detected %d entries out of rage in row/col indices; ignored",this->
IError());
239 Error(
"Decompose()",
"Matrix has not been set");
256 Error(
"Decompose()",
"nRows = %d out of range",
fNrows);
264 Info(
"Decompose()",
"insufficient space of fIw: %d",
fIw.
GetSize());
270 Info(
"Decompose()",
"resetting to fIw: %d",nIw);
284 Info(
"Decompose()",
"reseting to: %d",nFact);
290 Info(
"Decompose()",
"matrix apparently numerically singular");
291 Info(
"Decompose()",
"detected at stage %d",this->
IError());
292 Info(
"Decompose()",
"accept this factorization and hope for the best..");
298 Info(
"Decompose()",
"change of sign of pivots detected at stage %d",this->
IError());
299 Info(
"Decompose()",
"but who cares ");
304 Error(
"Decompose()",
"value of fNsteps out of range: %d",
fNsteps);
308 Info(
"Decompose()",
"detected %d entries out of range in row/column index",this->
IError());
309 Info(
"Decompose()",
"they are ignored");
315 Info(
"Decompose()",
"rank deficient matrix detected; apparent rank = %d",this->
IError());
323 }
while (!done && tries < 10);
326 if ( !done && tries >= 10) {
329 Error(
"Decompose()",
"did not get a factorization after 10 tries");
345 Error(
"Solve()",
"Matrix is singular");
350 Error(
"Solve()",
"Decomposition failed");
357 Error(
"Solve(TVectorD &",
"vector and matrix incompatible");
370 Int_t refactorizations = 0;
372 while (!done && refactorizations < 10) {
384 || refactorizations > 10) {
396 Info(
"Solve",
"Setting ThresholdPivoting parameter to %.4e for future factorizations",
419 const Int_t ifrlvl = 5;
434 fIcntl[ifrlvl+11] = 32639;
435 fIcntl[ifrlvl+12] = 32639;
436 fIcntl[ifrlvl+13] = 32639;
437 fIcntl[ifrlvl+14] = 32689;
475 Int_t i,iwfr,k,l1,l2,lliw;
484 for (i = 1; i < 16; i++)
487 if (icntl[3] > 0 && icntl[2] > 0) {
488 ::Info(
"TDecompSparse::InitPivot",
"Start with n = %d nz = %d liw = %d iflag = %d",n,nz,liw,iflag);
491 if (icntl[3] > 1) k = nz;
493 printf(
"matrix non-zeros:\n");
494 for (i = 1; i < k+1; i++) {
495 printf(
"%d %d ",irn[i],icn[i]);
496 if (i%5 == 0 || i == k) printf(
"\n");
501 if (icntl[3] > 1) k =
n;
502 if (iflag == 1 && k > 0) {
503 for (i = 1; i < k+1; i++) {
504 printf(
"%d ",ikeep[i]);
505 if (i%10 == 0 || i == k) printf(
"\n");
510 if (n >= 1 && n <= icntl[4]) {
514 ::Error(
"TDecompSparse::InitPivot",
"info[1]= %d; value of nz out of range .. = %d",info[1],nz);
521 if (liw < 2*nz+3*n+1) {
523 info[2] = 2*nz+3*n+1;
525 ::Error(
"TDecompSparse::InitPivot",
"info[1]= %d; liw too small, must be increased from %d to at least %d",info[1],liw,info[2]);
528 InitPivot_sub1(n,nz,irn,icn,iw,iw1,iw1+n+1,iw+l1-1,iwfr,icntl,info);
530 ikeep+2*(n+1),ikeep,icntl[4],info[11],cntl[2]);
532 if (liw < nz+3*n+1) {
536 ::Error(
"TDecompSparse::InitPivot",
"info[1]= %d; liw too small, must be increased from %d to at least %d",info[1],liw,info[2]);
539 InitPivot_sub3(n,nz,irn,icn,ikeep,iw,iw1,iw1+n+1,iw+l1-1,iwfr,icntl,info);
540 InitPivot_sub4(n,iw1,iw,lliw,iwfr,ikeep,ikeep+n+1,iw+l1-1,iw+l2-1,info[11]);
542 InitPivot_sub5(n,iw1,iw+l1-1,ikeep,ikeep+n+1,ikeep+2*(n+1),iw+l2-1,nsteps,icntl[5]);
543 if (nz >= 1) iw[1] = irn[1]+1;
544 InitPivot_sub6(n,nz,irn,icn,ikeep,ikeep+2*(n+1),ikeep+n+1,iw+l2-1,
545 nsteps,iw1,iw1+n+1,iw,info,ops);
549 ::Error(
"TDecompSparse::InitPivot",
"info[1]= %d; value of n out of range ... = %d",info[1],n);
553 if (icntl[3] <= 0 || icntl[2] <= 0)
return;
555 printf(
"Leaving with nsteps =%d info(1)=%d ops=%14.5e ierror=%d\n",nsteps,info[1],ops,info[2]);
556 printf(
"nrltot=%d nirtot=%d nrlnec=%d nirnec=%d nrladu=%d niradu=%d ncmpa=%d\n",
557 info[3],info[4],info[5],info[6],info[7],info[8],info[11]);
560 if (icntl[3] > 1) k =
n;
562 printf(
"ikeep[0][.]=\n");
563 for (i = 1; i < k+1; i++) {
564 printf(
"%d ",ikeep[i]);
565 if (k%10 == 0 || i == k) printf(
"\n");
570 printf(
"ikeep[2][.]=\n");
571 for (i = 1; i < k+1; i++) {
572 printf(
"%d ",ikeep[2*(n+1)+i]);
573 if (k%10 == 0 || i == k) printf(
"\n");
585 Int_t i,iapos,iblk,ipos,irows,j1,j2,jj,k,kblk,kz,len,ncols,nrows,nz1;
598 if (icntl[3] > 0 && icntl[2] > 0) {
599 printf(
"entering Factor with n=%d nz=%d la=%d liw=%d nsteps=%d u=%10.2e\n",
600 n,nz,la,liw,nsteps,cntl[1]);
602 if (icntl[3] > 1) kz = nz;
604 printf(
"matrix non-zeros:\n");
605 for (i = 1; i < kz+1; i++) {
606 printf(
"%16.3e %d %d ",a[i],irn[i],icn[i]);
607 if (i%2 == 0 || i==kz) printf(
"\n");
611 if (icntl[3] > 1) k =
n;
613 printf(
"ikeep(0,.)=\n");
614 for (i = 1; i < k+1; i++) {
615 printf(
"%d ",ikeep[i]);
616 if (i%10 == 0 || i == k) printf(
"\n");
621 printf(
"ikeep(1,.)=\n");
622 for (i = 1; i < k+1; i++) {
623 printf(
"%d ",ikeep[n+1+i]);
624 if (i%10 == 0 || i == k) printf(
"\n");
626 printf(
"ikeep(2,.)=\n");
627 for (i = 1; i < k+1; i++) {
628 printf(
"%d ",ikeep[2*(n+1)+i]);
629 if (i%10 == 0 || i == k) printf(
"\n");
634 if (n < 1 || n > icntl[4])
641 }
else if (la < nz+n) {
644 }
else if (nsteps < 1 || nsteps > n)
647 Factor_sub1(n,nz,nz1,a,la,irn,icn,iw,liw,ikeep,iw1,icntl,info);
648 if (info[1] != -3 && info[1] != -4) {
649 Factor_sub2(n,nz1,a,la,iw,liw,ikeep,ikeep+2*(n+1),nsteps,maxfrt,ikeep+n+1,iw1,icntl,cntl,info);
650 if (info[1] == 3 && icntl[2] > 0)
651 ::Warning(
"TDecompSparse::Factor",
"info[1]= %d; matrix is singular. rank=%d",info[1],info[2]);
658 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; value of n out of range ... =%d",info[1],n);
662 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; value of nz out of range ... =%d",info[1],nz);
666 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; liw too small, must be increased from %d to at least %d",info[1],liw,info[2]);
670 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; la too small, must be increased from %d to at least %d",info[1],la,info[2]);
674 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; zero pivot at stage %d zero pivot at stage",info[1],info[2]);
678 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; change in sign of pivot encountered when factoring allegedly definite matrix",info[1]);
682 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; nsteps is out of range",info[1]);
687 if (icntl[3] <= 0 || icntl[2] <= 0 || info[1] < 0)
690 ::Info(
"TDecompSparse::Factor",
"leaving Factor with maxfrt=%d info[1]=%d nrlbdu=%d nirbdu=%d ncmpbr=%d ncmpbi=%d ntwo=%d ierror=%d",
691 maxfrt,info[1],info[9],info[10],info[12],info[13],info[14],info[2]);
693 if (info[1] < 0)
return;
696 if (kblk == 0)
return;
697 if (icntl[3] == 1) kblk = 1;
701 for (iblk = 1; iblk < kblk+1; iblk++) {
710 ::Info(
"TDecompSparse::Factor",
"block pivot =%d nrows =%d ncols =%d",iblk,nrows,ncols);
714 printf(
" column indices =\n");
715 for (jj = j1; jj < j2+1; jj++) {
716 printf(
"%d ",iw[jj]);
717 if (jj%10 == 0 || jj == j2) printf(
"\n");
720 printf(
" real entries .. each row starts on a new line\n");
722 for (irows = 1; irows < nrows+1; irows++) {
725 for (jj = j1; jj < j2+1; jj++) {
726 printf(
"%13.4e ",a[jj]);
727 if (jj%5 == 0 || jj == j2) printf(
"\n");
742 Int_t i,iapos,iblk,ipos,irows,j1,j2,jj,k,kblk,latop,len,nblk,ncols,nrows;
756 if (icntl[3] > 0 && icntl[2] > 0) {
757 printf(
"nentering Solve with n=%d la=%d liw=%d maxfrt=%d nsteps=%d",n,la,liw,maxfrt,nsteps);
761 if (icntl[3] == 1) kblk = 1;
764 for (iblk = 1; iblk < kblk+1; iblk++) {
773 printf(
"block pivot=%d nrows=%d ncols=%d\n",iblk,nrows,ncols);
776 printf(
"column indices =\n");
777 for (jj = j1; jj < j2+1; jj++) {
778 printf(
"%d ",iw[jj]);
779 if (jj%10 == 0 || jj == j2) printf(
"\n");
781 printf(
"real entries .. each row starts on a new line\n");
783 for (irows = 1; irows < nrows+1; irows++) {
786 for (jj = j1; jj < j2+1; jj++) {
787 printf(
"%13.3e ",a[jj]);
788 if (jj%5 == 0 || jj == j2) printf(
"\n");
797 if (icntl[3] > 1) k =
n;
800 for (i = 1; i < k+1; i++) {
801 printf(
"%13.3e ",rhs[i]);
802 if (i%5 == 0 || i == k) printf(
"\n");
810 for (i = 1; i < n+1; i++)
813 nblk = (iw[1] <= 0) ? -iw[1] : iw[1];
814 Solve_sub1(n,a,iw+1,w,rhs,iw1,nblk,latop,icntl);
815 Solve_sub2(n,a,iw+1,w,rhs,iw1,nblk,latop,icntl);
818 if (icntl[3] > 0 && icntl[2] > 0) {
819 printf(
"leaving Solve with:\n");
822 for (i = 1; i < k+1; i++) {
823 printf(
"%13.3e ",rhs[i]);
824 if (i%5 == 0 || i == k) printf(
"\n");
840 Int_t i,id,j,jn,k,k1,k2,
l,last,lr,n1,ndup;
843 for (i = 1; i < n+1; i++)
848 for (k = 1; k < nz+1; k++) {
852 Bool_t outRange = (i < 1 || i > n || j < 1 || j >
n);
856 if (info[2] <= 1 && icntl[2]> 0)
857 ::Warning(
"TDecompSparse::InitPivot_sub1",
"info[1]= %d; %d th non-zero (in row=%d and column=%d) ignored",info[1],k,i,j);
860 if (outRange || i == j) {
876 for (i = 1; i < n1+1; i++) {
878 if (ipe[i] == 0) ipe[i] = -1;
879 iq[i+1] = ipe[i]+iq[i]+1;
888 for (k = k1; k < last+1; k++)
894 for (k = 1; k < nz+1; k++) {
896 if (j <= 0)
continue;
899 for (
id = 1;
id < nz+1;
id++) {
925 for (i = 1; i < n+1; i++) {
932 for (k = k1; k < k2+1; k++) {
947 iq[i] = iq[i]-ipe[i];
948 if (ndup == 0) iw[k1-1] = iq[i];
954 for (i = 1; i < n+1; i++) {
956 if (k1 == 1)
continue;
961 for (k = k1; k < k2+1; k++) {
962 if (iw[k] == 0)
continue;
980 Int_t i,id,idl,idn,ie,ip,is,jp,jp1,jp2,js,k,k1,k2,ke,kp,kp0,kp1,
981 kp2,ks,
l,len,limit,ln,
ls,lwfr,md,me,ml,ms,nel,nflg,np,
982 np0,ns,nvpiv,nvroot,root;
984 for (i = 1; i < n+1; i++) {
998 for (is = 1; is < n+1; is++) {
1003 if (ns > 0) lst[ns] = is;
1015 for (ml = 1; ml < n+1; ml++) {
1016 if (nel+nvroot+1 >= n)
break;
1017 for (
id = md;
id < n+1;
id++) {
1027 if (ns > 0) lst[ns] = -id;
1037 for (kp1 = 1; kp1 < len+1; kp1++) {
1040 if (flag[ke] > -2) {
1041 if (flag[ke] <= 0) {
1042 if (ipe[ke] != -root)
continue;
1044 if (flag[ke] <= 0)
continue;
1055 for (jp1 = 1; jp1 < ln+1; jp1++) {
1058 if (flag[is] <= 0) {
1059 if (ipe[is] == -root) {
1062 if (flag[is] <= 0)
continue;
1076 for (jp = ip; jp < jp2+1; jp++) {
1092 if (ns > 0) lst[ns] =
ls;
1112 for (k = k1; k < k2+1; k++) {
1114 if (is == root)
continue;
1116 for (i = 1; i < n+1; i++) {
1117 if (flag[i] > 0) flag[i] = iovflo;
1118 if (flag[i] <= -2) flag[i] = -iovflo;
1126 kp2 = iw[kp1-1]+kp1-1;
1129 for (kp = kp1; kp < kp2+1; kp++) {
1131 if (flag[ke] == -1) {
1132 if (ipe[ke] != -root)
continue;
1135 if (flag[ke] == -1)
continue;
1137 if (flag[ke] >= 0) {
1142 jp2 = iw[jp1-1]+jp1-1;
1144 for (jp = jp1; jp < jp2+1; jp++) {
1146 if (flag[js] <= nflg)
continue;
1152 for (jp = jp1; jp < jp2+1; jp++) {
1154 if (flag[js] != 0) {
1179 for (kp = kp0; kp < kp2+1; kp++) {
1181 if (flag[ks] <= nflg) {
1182 if (ipe[ks] == -root) {
1185 if (flag[ks] <= nflg)
continue;
1201 iw[kp1-1] = np-kp1+1;
1203 for (l = 1; l < n+1; l++) {
1209 if (iw[kp1] != me) {
1213 kp2 = kp1-1+iw[kp1-1];
1214 Int_t stayInLoop = 0;
1215 for (kp = kp1; kp < kp2+1; kp++) {
1232 nv[is] = nv[is]+nv[js];
1237 if (ns > 0) lst[ns] = is;
1238 if (ls > 0) nxt[
ls] = is;
1243 if (ipd[
id] == js) ipd[id] = is;
1244 }
else if (doit == 2) {
1251 nv[root] = nv[root]+nv[is];
1256 }
else if (doit == 3) {
1258 if (ns > 0) lst[ns] = is;
1266 for (k = k1; k < k2+1; k++) {
1268 if (nv[is] == 0)
continue;
1283 for (is = 1; is < n+1; is++) {
1284 if (nxt[is] != 0 || lst[is] != 0) {
1291 nvroot = nvroot+nv[is];
1296 for (ie = 1; ie < n+1; ie++)
1297 if (ipe[ie] > 0) ipe[ie] = -root;
1299 if (nvroot> 0) nv[root] = nvroot;
1308 Int_t i,ir,k,k1,k2,lwfr;
1311 for (i = 1; i < n+1; i++) {
1313 if (k1 <= 0)
continue;
1320 for (ir = 1; ir < n+1; ir++) {
1321 if (lwfr > lw)
break;
1323 for (k = lwfr; k < lw+1; k++) {
1337 for (k = k1; k < k2+1; k++) {
1353 Int_t i,id,in,j,jdummy,k,k1,k2,
l,lbig,len;
1357 for (i = 1; i < n+1; i++)
1361 for (k = 1; k < nz+1; k++) {
1366 Bool_t outRange = (i < 1 || i > n || j < 1 || j >
n);
1368 info[2] = info[2]+1;
1370 if (info[2] <= 1 && icntl[2] > 0)
1371 ::Warning(
"TDecompSparse::InitPivot_sub3",
"info[1]= %d; %d 'th non-zero (in row %d and column %d) ignored",info[1],k,i,j);
1374 if (outRange || i==j) {
1377 if (perm[j] <= perm[i])
1387 for (i = 1; i < n+1; i++) {
1395 for (k = 1; k < nz+1; k++) {
1397 if (i <= 0)
continue;
1400 for (
id = 1;
id < nz+1;
id++) {
1402 if (perm[i] >= perm[j]) {
1414 if (i <= 0)
continue;
1421 for (i = 1; i < n+1; i++) {
1426 for (jdummy = 1; jdummy < len+1; jdummy++) {
1436 if (lbig < icntl[4]) {
1437 for (i = 1; i < n+1; i++) {
1440 if (iq[i] == 0) ipe[i] = 0;
1444 for (i = 1; i < n+1; i++) {
1452 for (k = k1; k < k2+1; k++) {
1454 if (flag[j] == i)
continue;
1475 Int_t i,ie,ip,j,je,jp,jp1,jp2,js,kdummy,ln,lwfr,me,minjs,ml,ms;
1477 for (i = 1; i < n+1; i++) {
1485 for (ml = 1; ml < n+1; ml++) {
1493 for (kdummy = 1; kdummy < n+1; kdummy++) {
1498 for (jp1 = 1; jp1 < ln+1; jp1++) {
1501 if (flag[js] == me)
continue;
1510 for (jp = ip; jp < jp2+1; jp++) {
1551 Int_t i,ib,iff,il,is,ison,k,
l,nr;
1554 for (i = 1; i < n+1; i++) {
1558 for (i = 1; i < n+1; i++) {
1559 if (nv[i] > 0)
continue;
1562 if (is > 0) ipe[i] = is;
1567 for (i = 1; i < n+1; i++) {
1568 if (nv[i] <= 0)
continue;
1583 for (k = 1; k < n+1; k++) {
1591 for (l = 1; l < n+1; l++) {
1592 if (ips[i] >= 0)
break;
1603 if (il < n) na[il+1] = na[il+1]+1;
1607 Bool_t doit = (na[is] == 1 && (nd[is-1]-ne[is-1] == nd[is])) ||
1608 (na[is] != 1 && ne[is] < nemin && na[is] != 0 && ne[is-1] < nemin);
1611 na[is-1] = na[is-1]+na[is]-1;
1612 nd[is-1] = nd[is]+ne[is-1];
1613 ne[is-1] = ne[is]+ne[is-1];
1641 Int_t i,inew,iold,iorg,irow,istki,istkr,itop,itree,jold,jorg,k,lstk,nassr,nelim,nfr,nstk,
1642 numorg,nz1,nz2,nrladu,niradu,nirtot,nrltot,nirnec,nrlnec;
1645 if (nz != 0 && irn[1] == iw[1]) {
1648 for (iold = 1; iold < n+1; iold++) {
1650 lstki[inew] = lstkr[iold]+1;
1651 nz2 = nz2+lstkr[iold];
1656 for (i = 1; i < n+1; i++)
1660 for (i = 1; i < nz+1; i++) {
1663 if (iold < 1 || iold > n)
continue;
1664 if (jold < 1 || jold > n)
continue;
1665 if (iold == jold)
continue;
1667 irow =
TMath::Min(perm[iold]+0,perm[jold]+0);
1668 lstki[irow] = lstki[irow]+1;
1685 for (itree = 1; itree < nsteps+1; itree++) {
1690 nassr = nfr*(nfr+1)/2;
1691 if (nstk != 0) nassr = nassr-lstkr[itop]+1;
1692 nrltot =
TMath::Max(nrltot,nrladu+nassr+istkr+nz1);
1693 nirtot =
TMath::Max(nirtot,niradu+nfr+2+istki+nz1);
1694 nrlnec =
TMath::Max(nrlnec,nrladu+nassr+istkr+nz2);
1695 nirnec =
TMath::Max(nirnec,niradu+nfr+2+istki+nz2);
1696 for (iorg = 1; iorg < nelim+1; iorg++) {
1698 nz2 = nz2-lstki[jorg];
1700 numorg = numorg+nelim;
1702 for (k = 1; k < nstk+1; k++) {
1710 nrladu = nrladu+(nelim* (2*nfr-nelim+1))/2;
1711 niradu = niradu+2+nfr;
1712 if (nelim == 1) niradu = niradu-1;
1713 ops = ops+((nfr*delim*(nfr+1)-(2*nfr+1)*delim*(delim+1)/2+delim*(delim+1)*(2*delim+1)/6)/2);
1714 if (itree == nsteps || nfr == nelim)
continue;
1716 lstkr[itop] = (nfr-nelim)* (nfr-nelim+1)/2;
1717 lstki[itop] = nfr-nelim+1;
1718 istki = istki+lstki[itop];
1719 istkr = istkr+lstkr[itop];
1720 nirtot =
TMath::Max(nirtot,niradu+istki+nz1);
1721 nirnec =
TMath::Max(nirnec,niradu+istki+nz2);
1745 Int_t i,ia,ich,ii,iiw,inew,iold,ipos,j1,j2,jj,jnew,jold,jpos,k;
1751 for (iold = 1; iold < n+1; iold++) {
1760 for (k = 1; k < nz+1; k++) {
1763 Bool_t outRange = (iold < 1 || iold > n || jold < 1 || jold >
n);
1768 if (!outRange && inew == jnew) {
1775 iw2[inew] = iw2[inew]+1;
1780 info[2] = info[2]+1;
1781 if (info[2] <= 1 && icntl[2] > 0)
1782 ::Warning(
"TDecompSparse::Factor_sub1",
"info[1]= %d; %d 'th non-zero (in row %d and column %d) ignored",
1783 info[1],k,irn[k],icn[k]);
1790 if (nz >= nz1 || nz1 == n) {
1792 for (i = 1; i < n+1; i++) {
1798 for (i = 1; i < n+1; i++) {
1817 for (k = 1; k < nz+1; k++) {
1819 if (iold <= 0)
continue;
1823 for (ich = 1; ich < nz+1; ich++) {
1827 if (inew == perm[jold]) jold = iold;
1834 if (iold == 0)
break;
1843 for (ii = 1; ii < n+1; ii++) {
1848 for (jj = j1; jj < j2+1; jj++) {
1849 iw[ipos] = iw[jpos];
1861 for (iold = 1; iold < n+1; iold++) {
1871 for (i = 1; i < nz1+1; i++) {
1888 Double_t amax,amult,amult1,amult2,detpiv,rmax,swop,thresh,tmax,uu;
1889 Int_t ainput,apos,apos1,apos2,apos3,astk,astk2,azero,i,iass;
1890 Int_t ibeg,idummy,iell,iend,iexch,ifr,iinput,ioldps,iorg,ipiv;
1891 Int_t ipmnp,ipos,irow,isnpiv,istk,istk2,iswop,iwpos,j,j1;
1892 Int_t j2,jcol,jdummy,jfirst,jj,jj1,jjj,jlast,jmax,jmxmip,jnew;
1893 Int_t jnext,jpiv,jpos,k,k1,k2,kdummy,kk,kmax,krow,laell,lapos2;
1894 Int_t liell,lnass,lnpiv,lt,ltopst,nass,nblk,newel,nfront,npiv;
1895 Int_t npivp1,ntotpv,numass,numorg,numstk,pivsiz,posfac,pospv1,pospv2;
1896 Int_t ntwo,neig,ncmpbi,ncmpbr,nrlbdu,nirbdu;
1918 for (i = 1; i < n+1; i++)
1932 for (iass = 1; iass < nsteps+1; iass++) {
1937 numstk = nstk[iass];
1943 ltopst = ((iw[istk]+1)*iw[istk])/2;
1944 for (iell = 1; iell < numstk+1; iell++) {
1949 for (jj = j1; jj < j2+1; jj++) {
1951 if (iw2[j] > 0)
continue;
1953 if (jnew <= numass) nass = nass+1;
1954 for (idummy = 1; idummy < n+1; idummy++) {
1955 if (jnext == n+1)
break;
1956 if (perm[jnext] > jnew)
break;
1972 numorg = nelim[iass];
1974 for (iorg = 1; iorg < numorg+1; iorg++) {
1976 for (idummy = 1; idummy < liw+1; idummy++) {
1981 for (jdummy = 1; jdummy < n+1; jdummy++) {
1982 if (jnext == n+1)
break;
1983 if (perm[jnext] > jnew)
break;
1995 if (j1 > liw)
break;
2001 if (newel+nfront >= istk)
2002 Factor_sub3(a,iw,istk,istk2,iinput,2,ncmpbr,ncmpbi);
2003 if (newel+nfront >= istk) {
2005 info[2] = liw+1+newel+nfront-istk;
2010 for (ifr = 1; ifr < nfront+1; ifr++) {
2014 iw2[j] = newel-(iwpos+1);
2020 laell = ((nfront+1)*nfront)/2;
2021 apos2 = posfac+laell-1;
2022 if (numstk != 0) lnass = lnass*(2*nfront-lnass+1)/2;
2024 if (posfac+lnass-1 >= astk || apos2 >= astk+ltopst-1) {
2025 Factor_sub3(a,iw,astk,astk2,ainput,1,ncmpbr,ncmpbi);
2026 if (posfac+lnass-1 >= astk || apos2 >= astk+ltopst-1) {
2028 info[2] = la+
TMath::Max(posfac+lnass,apos2-ltopst+2)-astk;
2033 if (apos2 > azero) {
2036 if (lapos2 >= apos) {
2037 for (k= apos; k< lapos2+1; k++)
2044 for (iell = 1; iell < numstk+1; iell++) {
2047 for (jj = j1; jj < j2+1; jj++) {
2050 apos = posfac+
IDiag(nfront,irow);
2051 for (jjj = jj; jjj < j2+1; jjj++) {
2053 apos2 = apos+iw2[j]-irow;
2054 a[apos2] = a[apos2]+a[astk];
2063 for (iorg = 1; iorg < numorg+1; iorg++) {
2066 apos = posfac+
IDiag(nfront,irow);
2067 for (idummy = 1; idummy < nz+1; idummy++) {
2068 apos2 = apos+iw2[j]-irow;
2069 a[apos2] = a[apos2]+a[ainput];
2072 if (iinput > liw)
break;
2077 numass = numass+numorg;
2079 j2 = iwpos+nfront+1;
2080 for (k = j1; k < j2+1; k++) {
2087 for (kdummy = 1; kdummy < nass+1; kdummy++) {
2088 if (npiv == nass)
break;
2089 if (npiv == lnpiv)
break;
2093 for (ipiv = npivp1; ipiv < nass+1; ipiv++) {
2095 if (jpiv == 1)
continue;
2096 apos = posfac+
IDiag(nfront-npiv,ipiv-npiv);
2105 if (a[apos] > zero) isnpiv = 1;
2106 if (a[apos] < zero) isnpiv = -1;
2108 if ((a[apos] <= zero || isnpiv != 1) && (a[apos] >= zero || isnpiv != -1)) {
2109 if (info[1] != 2) info[2] = 0;
2110 info[2] = info[2]+1;
2113 if (icntl[2] > 0 && info[2] <= 10)
2114 ::Warning(
"TDecompSparse::Factor_sub2",
"info[1]= %d; pivot %d has different sign from the previous one",
2118 if ((a[apos] > zero && isnpiv == 1) || (a[apos] < zero && isnpiv == -1) || (uu == zero))
goto hack;
2127 j2 = apos+nass-ipiv;
2129 for (jj = j1; jj < j2+1; jj++) {
2131 jmax = ipiv+jj-j1+1;
2136 j2 = apos+nfront-ipiv;
2138 for (jj = j1; jj < j2+1; jj++)
2146 for (k = 1; k < lt+1; k++) {
2154 apos2 = posfac+
IDiag(nfront-npiv,jmax-npiv);
2155 detpiv = a[apos]*a[apos2]-amax*amax;
2158 if (thresh <= rmax)
continue;
2161 j2 = apos2+nfront-jmax;
2163 for (jj = j1; jj < j2+1; jj++)
2168 jmxmip = jmax-ipiv-1;
2170 for (k = 1; k < jmxmip+1; k++) {
2176 ipmnp = ipiv-npiv-1;
2180 for (k = 1; k < ipmnp+1; k++) {
2186 if (thresh <= rmax)
continue;
2193 for (krow = 1; krow < pivsiz+1; krow++) {
2196 j2 = posfac+nfront-(npiv+1);
2199 for (jj = j1; jj < j2+1; jj++) {
2209 kk = nfront-(irow+npiv);
2211 for (jjj = j1; jjj < j2+1; jjj++) {
2224 for (jj = 1; jj < npiv+1; jj++) {
2229 a[apos2] = a[apos1];
2234 a[apos] = a[posfac];
2236 ipos = iwpos+npiv+2;
2237 iexch = iwpos+irow+npiv+1;
2239 iw[ipos] = iw[iexch];
2242 if (pivsiz == 1)
continue;
2244 irow = jmax-(npiv+1);
2246 posfac = posfac+nfront-npiv;
2257 a[posfac] = one/a[posfac];
2258 if (a[posfac] < zero) neig = neig+1;
2260 j2 = posfac+nfront-(npiv+1);
2263 for (jj = j1; jj < j2+1; jj++) {
2264 amult = -a[jj]*a[posfac];
2265 iend = ibeg+nfront-(npiv+jj-j1+2);
2266 for (irow = ibeg; irow < iend+1; irow++) {
2267 jcol = jj+irow-ibeg;
2268 a[irow] = a[irow]+amult*a[jcol];
2277 posfac = posfac+nfront-npiv+1;
2279 ipos = iwpos+npiv+2;
2281 iw[ipos] = -iw[ipos];
2283 pospv2 = posfac+nfront-npiv;
2285 if (detpiv < zero) neig = neig+1;
2286 if (detpiv > zero && swop < zero) neig = neig+2;
2287 a[pospv2] = a[pospv1]/detpiv;
2288 a[pospv1] = swop/detpiv;
2289 a[pospv1+1] = -a[pospv1+1]/detpiv;
2291 j2 = pospv1+nfront-(npiv+1);
2294 ibeg = pospv2+nfront-(npiv+1);
2295 for (jj = j1; jj < j2+1; jj++) {
2297 amult1 =-(a[pospv1]*a[jj]+a[pospv1+1]*a[jj1]);
2298 amult2 =-(a[pospv1+1]*a[jj]+a[pospv2]*a[jj1]);
2299 iend = ibeg+nfront-(npiv+jj-j1+3);
2300 for (irow = ibeg; irow < iend+1; irow++) {
2303 a[irow] = a[irow]+amult1*a[k1]+amult2*a[k2];
2313 posfac = pospv2+nfront-npiv+1;
2318 if (npiv != 0) nblk = nblk+1;
2320 iwpos = iwpos+nfront+2;
2323 iw[ioldps] = -iw[ioldps];
2324 for (k = 1; k < nfront+1; k++) {
2330 iw[ioldps+1] = npiv;
2333 liell = nfront-npiv;
2335 if (liell != 0 && iass != nsteps) {
2336 if (iwpos+liell >= istk)
2337 Factor_sub3(a,iw,istk,istk2,iinput,2,ncmpbr,ncmpbi);
2338 istk = istk-liell-1;
2342 for (k = 1; k < liell+1; k++) {
2347 laell = ((liell+1)*liell)/2;
2353 for (k = 1; k < laell+1; k++) {
2359 for ( k = kk; k < kmax+1; k++)
2364 if (npiv == 0) iwpos = ioldps;
2368 if (ntwo > 0) iw[1] = -nblk;
2399 for (jjj = j1; jjj < j2+1; jjj++) {
2408 for (jjj = j1; jjj < j2+1; jjj++) {
2427 Int_t apos,iblk,ifr,ilvl,ipiv,ipos,irhs,irow,ist,j,j1=0,j2,j3,jj,jpiv,k,k1,k2,k3,liell,npiv;
2430 const Int_t ifrlvl = 5;
2437 for (irow = 1; irow < n+1; irow++) {
2440 if (iblk > nblk)
break;
2453 if (liell < icntl[ifrlvl+ilvl])
goto hack;
2455 for (jj = j1; jj < j2+1; jj++) {
2463 for (ipiv = 1; ipiv < npiv+1; ipiv++) {
2465 if (jpiv == 1)
continue;
2472 if (liell< ist)
continue;
2475 for (j = ist; j < liell+1; j++) {
2476 w[j] = w[j]+a[k]*w1;
2479 apos = apos+liell-ist+1;
2489 k2 = apos+liell-ipiv;
2490 for (j = ist; j < liell+1; j++) {
2491 w[j] = w[j]+w1*a[k1]+w2*a[k2];
2496 apos = apos+2*(liell-ist+1)+1;
2501 for (jj = j1; jj < j2+1; jj++) {
2517 for (j = j1; j < j2+1; j++) {
2519 rhs[irhs] = rhs[irhs]+a[k]*w1;
2523 apos = apos+j2-j1+1;
2535 for (j = j1; j < j2+1; j++) {
2537 rhs[irhs] = rhs[irhs]+w1*a[k1]+w2*a[k3];
2542 apos = apos+2*(j2-j1+1)+1;
2557 Int_t apos,apos2,i1rhs,i2rhs,iblk,ifr,iipiv,iirhs,ilvl,ipiv,ipos,irhs,ist,
2558 j,j1=0,j2=0,jj,jj1,jj2,jpiv,jpos=0,k,liell,loop,npiv;
2561 const Int_t ifrlvl = 5;
2566 for (loop = 1; loop < n+1; loop++) {
2569 if (iblk < 1)
break;
2581 if (liell < icntl[ifrlvl+ilvl])
goto hack;
2584 for (jj = j1; jj < j2+1; jj++) {
2590 for (iipiv = 1; iipiv < npiv+1; iipiv++) {
2592 if (jpiv == 1)
continue;
2593 ipiv = npiv-iipiv+1;
2594 if (ipiv == 1 || iw[jpos-1] >= 0) {
2596 apos = apos-(liell+1-ipiv);
2598 w1 = w[ipiv]*a[apos];
2601 for (j = ist; j < liell+1; j++) {
2602 w1 = w1+a[jj1]*w[j];
2610 apos2 = apos-(liell+1-ipiv);
2611 apos = apos2-(liell+2-ipiv);
2613 w1 = w[ipiv-1]*a[apos]+w[ipiv]*a[apos+1];
2614 w2 = w[ipiv-1]*a[apos+1]+w[ipiv]*a[apos2];
2618 for (j = ist; j < liell+1; j++) {
2619 w1 = w1+w[j]*a[jj1];
2620 w2 = w2+w[j]*a[jj2];
2631 for (jj = j1; jj < j2+1; jj++) {
2639 if (npiv == 1 || iw[jpos-1] >= 0) {
2641 apos = apos-(j2-jpos+1);
2643 w1 = rhs[iirhs]*a[apos];
2647 for (j = j1; j < j2+1; j++) {
2649 w1 = w1+a[k]*rhs[irhs];
2657 apos2 = apos-(j2-jpos+1);
2658 apos = apos2-(j2-jpos+2);
2659 i1rhs = -iw[jpos-1];
2661 w1 = rhs[i1rhs]*a[apos]+rhs[i2rhs]*a[apos+1];
2662 w2 = rhs[i1rhs]*a[apos+1]+rhs[i2rhs]*a[apos2];
2667 for (j = j1; j < j2+1; j++) {
2669 w1 = w1+rhs[irhs]*a[jj1];
2670 w2 = w2+rhs[irhs]*a[jj2];
2696 fact.
Print(
"fFact");
2704 if (
this != &source) {
2721 fA.
Use(*const_cast<TMatrixDSparse *>(&(source.
fA)));
virtual const Element * GetMatrixArray() const
void SetTreatAsZero(Double_t tol)
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
virtual void SetMatrix(const TMatrixDSparse &a)
Set matrix to be decomposed .
virtual const Int_t * GetRowIndexArray() const
static void InitPivot_sub5(const Int_t n, Int_t *ipe, Int_t *nv, Int_t *ips, Int_t *ne, Int_t *na, Int_t *nd, Int_t &nsteps, const Int_t nemin)
Help routine for pivoting setup.
static void CopyUpperTriang(const TMatrixDSparse &a, Double_t *b)
Static function, copying the non-zero entries in the upper triangle to array b .
const Double_t * GetArray() const
static void InitPivot_sub1(const Int_t n, const Int_t nz, Int_t *irn, Int_t *icn, Int_t *iw, Int_t *ipe, Int_t *iq, Int_t *flag, Int_t &iwfr, Int_t *icntl, Int_t *info)
Help routine for pivoting setup.
const Double_t kThresholdPivotingFactor
Bool_t TestBit(UInt_t f) const
void Print(Option_t *opt="") const
Print class members.
TDecompSparse()
Default constructor.
Decomposition Base class.
static void Factor_sub3(Double_t *a, Int_t *iw, Int_t &j1, Int_t &j2, const Int_t itop, const Int_t ireal, Int_t &ncmpbr, Int_t &ncmpbi)
Help routine for factorization.
TDecompSparse & operator=(const TDecompSparse &source)
Assignment operator.
Short_t Min(Short_t a, Short_t b)
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Array of integers (32 bits per element).
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
virtual Bool_t Decompose()
Decomposition engine .
static Int_t IDiag(Int_t ix, Int_t iy)
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
static void Solve(const Int_t n, TArrayD &Aa, TArrayI &Aiw, TArrayD &Aw, const Int_t maxfrt, TVectorD &b, TArrayI &Aiw1, const Int_t nsteps, Int_t *icntl, Int_t *info)
Main routine for solving Ax=b.
static void Solve_sub1(const Int_t n, Double_t *a, Int_t *iw, Double_t *w, Double_t *rhs, Int_t *iw2, const Int_t nblk, Int_t &latop, Int_t *icntl)
Help routine for solving.
static void Factor(const Int_t n, const Int_t nz, TArrayI &Airn, TArrayI &Aicn, TArrayD &Aa, TArrayI &Aiw, TArrayI &Aikeep, const Int_t nsteps, Int_t &maxfrt, TArrayI &Aiw1, Int_t *icntl, Double_t *cntl, Int_t *info)
Factorization routine, the workhorse for the decomposition step.
void InitParam()
initializing control parameters
const Int_t * GetArray() const
void Set(Int_t n)
Set size of this array to n ints.
virtual void ls(Option_t *option="") const
The ls function lists the contents of a class on stdout.
TVectorT< Element > & Shift(Int_t row_shift)
Element * GetMatrixArray()
static Int_t NonZerosUpperTriang(const TMatrixDSparse &a)
Static function, returning the number of non-zero entries in the upper triangular matrix ...
static void Solve_sub2(const Int_t n, Double_t *a, Int_t *iw, Double_t *w, Double_t *rhs, Int_t *iw2, const Int_t nblk, const Int_t latop, Int_t *icntl)
Help routine for solving.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
void Print(Option_t *opt="") const
Print class members.
Double_t GetThresholdPivoting()
Sparse Symmetric Decomposition class.
void SetThresholdPivoting(Double_t piv)
static void Factor_sub2(const Int_t n, const Int_t nz, Double_t *a, const Int_t la, Int_t *iw, const Int_t liw, Int_t *perm, Int_t *nstk, const Int_t nsteps, Int_t &maxfrt, Int_t *nelim, Int_t *iw2, Int_t *icntl, Double_t *cntl, Int_t *info)
Help routine for factorization.
virtual const Int_t * GetColIndexArray() const
const Double_t kInitPrecision
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)
const Double_t kThresholdPivotingMax
Array of doubles (64 bits per element).
void Print(Option_t *name="") const
Print the matrix as a table of elements.
static void InitPivot_sub2(const Int_t n, Int_t *ipe, Int_t *iw, const Int_t lw, Int_t &iwfr, Int_t *nv, Int_t *nxt, Int_t *lst, Int_t *ipd, Int_t *flag, const Int_t iovflo, Int_t &ncmpa, const Double_t fratio)
Help routine for pivoting setup.
static void InitPivot_sub4(const Int_t n, Int_t *ipe, Int_t *iw, const Int_t lw, Int_t &iwfr, Int_t *ips, Int_t *ipv, Int_t *nv, Int_t *flag, Int_t &ncmpa)
Help routine for pivoting setup.
const Double_t kInitTreatAsZero
static void InitPivot_sub2a(const Int_t n, Int_t *ipe, Int_t *iw, const Int_t lw, Int_t &iwfr, Int_t &ncmpa)
Help routine for pivoting setup.
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
const Double_t kInitThresholdPivoting
static void InitPivot(const Int_t n, const Int_t nz, TArrayI &Airn, TArrayI &Aicn, TArrayI &Aiw, TArrayI &Aikeep, TArrayI &Aiw1, Int_t &nsteps, const Int_t iflag, Int_t *icntl, Double_t *cntl, Int_t *info, Double_t &ops)
Setup Pivoting variables.
static void InitPivot_sub6(const Int_t n, const Int_t nz, Int_t *irn, Int_t *icn, Int_t *perm, Int_t *na, Int_t *ne, Int_t *nd, const Int_t nsteps, Int_t *lstki, Int_t *lstkr, Int_t *iw, Int_t *info, Double_t &ops)
Help routine for pivoting setup.
void Set(Int_t n)
Set size of this array to n doubles.
static char * skip(char **buf, const char *delimiters)
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
static void Factor_sub1(const Int_t n, const Int_t nz, Int_t &nz1, Double_t *a, const Int_t la, Int_t *irn, Int_t *icn, Int_t *iw, const Int_t liw, Int_t *perm, Int_t *iw2, Int_t *icntl, Int_t *info)
Help routine for factorization.
static void InitPivot_sub3(const Int_t n, const Int_t nz, Int_t *irn, Int_t *icn, Int_t *perm, Int_t *iw, Int_t *ipe, Int_t *iq, Int_t *flag, Int_t &iwfr, Int_t *icntl, Int_t *info)
Help routine for pivoting setup.