202 Error(
"SetMatrix(const TMatrixDSparse &",
"nRows = %d out of range",fNrows);
205 Error(
"SetMatrix(const TMatrixDSparse &",
"nr_nonzeros = %d out of range",
fNnonZeros);
208 Error(
"SetMatrix(const TMatrixDSparse &",
209 "insufficient space in fIw of %d suggest reset to %d",
fIw.
GetSize(),
this->IError());
212 Error(
"SetMatrix(const TMatrixDSparse &",
213 "detected %d entries out of range in row/col indices; ignored",this->
IError());
238 Error(
"Decompose()",
"Matrix has not been set");
255 Error(
"Decompose()",
"nRows = %d out of range",
fNrows);
263 Info(
"Decompose()",
"insufficient space of fIw: %d",
fIw.
GetSize());
269 Info(
"Decompose()",
"resetting to fIw: %d",
nIw);
283 Info(
"Decompose()",
"resetting to: %d",
nFact);
289 Info(
"Decompose()",
"matrix apparently numerically singular");
290 Info(
"Decompose()",
"detected at stage %d",this->
IError());
291 Info(
"Decompose()",
"accept this factorization and hope for the best..");
297 Info(
"Decompose()",
"change of sign of pivots detected at stage %d",this->
IError());
298 Info(
"Decompose()",
"but who cares ");
303 Error(
"Decompose()",
"value of fNsteps out of range: %d",
fNsteps);
307 Info(
"Decompose()",
"detected %d entries out of range in row/column index",this->
IError());
308 Info(
"Decompose()",
"they are ignored");
314 Info(
"Decompose()",
"rank deficient matrix detected; apparent rank = %d",this->
IError());
328 Error(
"Decompose()",
"did not get a factorization after 10 tries");
344 Error(
"Solve()",
"Matrix is singular");
349 Error(
"Solve()",
"Decomposition failed");
356 Error(
"Solve(TVectorD &",
"vector and matrix incompatible");
383 || refactorizations > 10) {
395 Info(
"Solve",
"Setting ThresholdPivoting parameter to %.4e for future factorizations",
483 for (i = 1; i < 16; i++)
487 ::Info(
"TDecompSparse::InitPivot",
"Start with n = %d nz = %d liw = %d iflag = %d",
n,
nz,
liw,
iflag);
492 printf(
"matrix non-zeros:\n");
493 for (i = 1; i < k+1; i++) {
495 if (i%5 == 0 || i == k)
printf(
"\n");
501 if (
iflag == 1 && k > 0) {
502 for (i = 1; i < k+1; i++) {
504 if (i%10 == 0 || i == k)
printf(
"\n");
513 ::Error(
"TDecompSparse::InitPivot",
"info[1]= %d; value of nz out of range .. = %d",
info[1],
nz);
524 ::Error(
"TDecompSparse::InitPivot",
"info[1]= %d; liw too small, must be increased from %d to at least %d",
info[1],
liw,
info[2]);
527 InitPivot_sub1(
n,
nz,
irn,
icn,
iw,
iw1,
iw1+
n+1,
iw+
l1-1,
iwfr,
icntl,
info);
535 ::Error(
"TDecompSparse::InitPivot",
"info[1]= %d; liw too small, must be increased from %d to at least %d",
info[1],
liw,
info[2]);
538 InitPivot_sub3(
n,
nz,
irn,
icn,
ikeep,
iw,
iw1,
iw1+
n+1,
iw+
l1-1,
iwfr,
icntl,
info);
548 ::Error(
"TDecompSparse::InitPivot",
"info[1]= %d; value of n out of range ... = %d",
info[1],
n);
555 printf(
"nrltot=%d nirtot=%d nrlnec=%d nirnec=%d nrladu=%d niradu=%d ncmpa=%d\n",
562 for (i = 1; i < k+1; i++) {
564 if (k%10 == 0 || i == k)
printf(
"\n");
570 for (i = 1; i < k+1; i++) {
572 if (k%10 == 0 || i == k)
printf(
"\n");
584 Int_t i,
iapos,
iblk,
ipos,
irows,
j1,
j2,
jj,k,
kblk,
kz,
len,
ncols,
nrows,
nz1;
598 printf(
"entering Factor with n=%d nz=%d la=%d liw=%d nsteps=%d u=%10.2e\n",
603 printf(
"matrix non-zeros:\n");
604 for (i = 1; i <
kz+1; i++) {
606 if (i%2 == 0 || i==
kz)
printf(
"\n");
613 for (i = 1; i < k+1; i++) {
615 if (i%10 == 0 || i == k)
printf(
"\n");
621 for (i = 1; i < k+1; i++) {
623 if (i%10 == 0 || i == k)
printf(
"\n");
626 for (i = 1; i < k+1; i++) {
628 if (i%10 == 0 || i == k)
printf(
"\n");
640 }
else if (
la <
nz+
n) {
646 Factor_sub1(
n,
nz,
nz1,
a,
la,
irn,
icn,
iw,
liw,
ikeep,
iw1,
icntl,
info);
647 if (
info[1] != -3 &&
info[1] != -4) {
648 Factor_sub2(
n,
nz1,
a,
la,
iw,
liw,
ikeep,
ikeep+2*(
n+1),
nsteps,
maxfrt,
ikeep+
n+1,
iw1,
icntl,
cntl,
info);
650 ::Warning(
"TDecompSparse::Factor",
"info[1]= %d; matrix is singular. rank=%d",
info[1],
info[2]);
657 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; value of n out of range ... =%d",
info[1],
n);
661 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; value of nz out of range ... =%d",
info[1],
nz);
665 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; liw too small, must be increased from %d to at least %d",
info[1],
liw,
info[2]);
669 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; la too small, must be increased from %d to at least %d",
info[1],
la,
info[2]);
673 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; zero pivot at stage %d zero pivot at stage",
info[1],
info[2]);
677 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; change in sign of pivot encountered when factoring allegedly definite matrix",
info[1]);
681 ::Error(
"TDecompSparse::Factor",
"info[1]= %d; nsteps is out of range",
info[1]);
689 ::Info(
"TDecompSparse::Factor",
"leaving Factor with maxfrt=%d info[1]=%d nrlbdu=%d nirbdu=%d ncmpbr=%d ncmpbi=%d ntwo=%d ierror=%d",
692 if (
info[1] < 0)
return;
695 if (
kblk == 0)
return;
713 printf(
" column indices =\n");
719 printf(
" real entries .. each row starts on a new line\n");
741 Int_t i,
iapos,
iblk,
ipos,
irows,
j1,
j2,
jj,k,
kblk,
latop,
len,
nblk,
ncols,
nrows;
775 printf(
"column indices =\n");
780 printf(
"real entries .. each row starts on a new line\n");
799 for (i = 1; i < k+1; i++) {
801 if (i%5 == 0 || i == k)
printf(
"\n");
809 for (i = 1; i <
n+1; i++)
818 printf(
"leaving Solve with:\n");
821 for (i = 1; i < k+1; i++) {
823 if (i%5 == 0 || i == k)
printf(
"\n");
839 Int_t i,
id,
j,
jn,k,
k1,
k2,
l,last,
lr,
n1,
ndup;
842 for (i = 1; i <
n+1; i++)
847 for (k = 1; k <
nz+1; k++) {
856 ::Warning(
"TDecompSparse::InitPivot_sub1",
"info[1]= %d; %d th non-zero (in row=%d and column=%d) ignored",
info[1],k,i,
j);
875 for (i = 1; i <
n1+1; i++) {
877 if (
ipe[i] == 0)
ipe[i] = -1;
887 for (k =
k1; k < last+1; k++)
893 for (k = 1; k <
nz+1; k++) {
895 if (
j <= 0)
continue;
898 for (
id = 1;
id <
nz+1;
id++) {
924 for (i = 1; i <
n+1; i++) {
931 for (k =
k1; k <
k2+1; k++) {
953 for (i = 1; i <
n+1; i++) {
955 if (
k1 == 1)
continue;
960 for (k =
k1; k <
k2+1; k++) {
961 if (
iw[k] == 0)
continue;
979 Int_t i,
id,
idl,
idn,
ie,
ip,
is,
jp,
jp1,
jp2,
js,k,
k1,
k2,
ke,
kp,
kp0,
kp1,
980 kp2,
ks,
l,
len,limit,
ln,
ls,
lwfr,
md,me,
ml,ms,
nel,
nflg,
np,
983 for (i = 1; i <
n+1; i++) {
1002 if (ns > 0)
lst[ns] =
is;
1014 for (
ml = 1;
ml <
n+1;
ml++) {
1016 for (
id =
md;
id <
n+1;
id++) {
1026 if (ns > 0)
lst[ns] = -
id;
1041 if (
ipe[
ke] != -root)
continue;
1043 if (
flag[
ke] <= 0)
continue;
1058 if (
ipe[
is] == -root) {
1061 if (
flag[
is] <= 0)
continue;
1091 if (ns > 0)
lst[ns] =
ls;
1111 for (k =
k1; k <
k2+1; k++) {
1113 if (
is == root)
continue;
1115 for (i = 1; i <
n+1; i++) {
1131 if (
ipe[
ke] != -root)
continue;
1134 if (
flag[
ke] == -1)
continue;
1181 if (
ipe[
ks] == -root) {
1202 for (
l = 1;
l <
n+1;
l++) {
1208 if (
iw[
kp1] != me) {
1236 if (ns > 0)
lst[ns] =
is;
1243 }
else if (
doit == 2) {
1255 }
else if (
doit == 3) {
1257 if (ns > 0)
lst[ns] =
is;
1265 for (k =
k1; k <
k2+1; k++) {
1267 if (
nv[
is] == 0)
continue;
1282 for (
is = 1;
is <
n+1;
is++) {
1310 for (i = 1; i <
n+1; i++) {
1312 if (
k1 <= 0)
continue;
1319 for (
ir = 1;
ir <
n+1;
ir++) {
1322 for (k =
lwfr; k <
lw+1; k++) {
1336 for (k =
k1; k <
k2+1; k++) {
1352 Int_t i,
id,
in,
j,
jdummy,k,
k1,
k2,
l,
lbig,
len;
1356 for (i = 1; i <
n+1; i++)
1360 for (k = 1; k <
nz+1; k++) {
1370 ::Warning(
"TDecompSparse::InitPivot_sub3",
"info[1]= %d; %d 'th non-zero (in row %d and column %d) ignored",
info[1],k,i,
j);
1386 for (i = 1; i <
n+1; i++) {
1394 for (k = 1; k <
nz+1; k++) {
1396 if (i <= 0)
continue;
1399 for (
id = 1;
id <
nz+1;
id++) {
1413 if (i <= 0)
continue;
1420 for (i = 1; i <
n+1; i++) {
1436 for (i = 1; i <
n+1; i++) {
1439 if (
iq[i] == 0)
ipe[i] = 0;
1443 for (i = 1; i <
n+1; i++) {
1451 for (k =
k1; k <
k2+1; k++) {
1453 if (
flag[
j] == i)
continue;
1474 Int_t i,
ie,
ip,
j,
je,
jp,
jp1,
jp2,
js,
kdummy,
ln,
lwfr,me,
minjs,
ml,ms;
1476 for (i = 1; i <
n+1; i++) {
1484 for (
ml = 1;
ml <
n+1;
ml++) {
1500 if (
flag[
js] == me)
continue;
1553 for (i = 1; i <
n+1; i++) {
1557 for (i = 1; i <
n+1; i++) {
1558 if (
nv[i] > 0)
continue;
1566 for (i = 1; i <
n+1; i++) {
1567 if (
nv[i] <= 0)
continue;
1582 for (k = 1; k <
n+1; k++) {
1590 for (
l = 1;
l <
n+1;
l++) {
1591 if (
ips[i] >= 0)
break;
1640 Int_t i,
inew,
iold,
iorg,
irow,
istki,
istkr,
itop,
itree,
jold,
jorg,k,
lstk,
nassr,
nelim,
nfr,
nstk,
1644 if (
nz != 0 &&
irn[1] ==
iw[1]) {
1655 for (i = 1; i <
n+1; i++)
1659 for (i = 1; i <
nz+1; i++) {
1701 for (k = 1; k <
nstk+1; k++) {
1744 Int_t i,
ia,
ich,
ii,
iiw,
inew,
iold,
ipos,
j1,
j2,
jj,
jnew,
jold,
jpos,k;
1759 for (k = 1; k <
nz+1; k++) {
1781 ::Warning(
"TDecompSparse::Factor_sub1",
"info[1]= %d; %d 'th non-zero (in row %d and column %d) ignored",
1791 for (i = 1; i <
n+1; i++) {
1797 for (i = 1; i <
n+1; i++) {
1816 for (k = 1; k <
nz+1; k++) {
1818 if (
iold <= 0)
continue;
1833 if (
iold == 0)
break;
1842 for (
ii = 1;
ii <
n+1;
ii++) {
1870 for (i = 1; i <
nz1+1; i++) {
1891 Int_t j2,
jcol,
jdummy,
jfirst,
jj,
jj1,
jjj,
jlast,
jmax,
jmxmip,
jnew;
1892 Int_t jnext,
jpiv,
jpos,k,
k1,
k2,
kdummy,
kk,kmax,
krow,
laell,
lapos2;
1917 for (i = 1; i <
n+1; i++)
1950 if (
iw2[
j] > 0)
continue;
1994 if (
j1 >
liw)
break;
2079 for (k =
j1; k <
j2+1; k++) {
2094 if (
jpiv == 1)
continue;
2113 ::Warning(
"TDecompSparse::Factor_sub2",
"info[1]= %d; pivot %d has different sign from the previous one",
2145 for (k = 1; k <
lt+1; k++) {
2169 for (k = 1; k <
jmxmip+1; k++) {
2179 for (k = 1; k <
ipmnp+1; k++) {
2241 if (
pivsiz == 1)
continue;
2323 for (k = 1; k <
nfront+1; k++) {
2341 for (k = 1; k <
liell+1; k++) {
2352 for (k = 1; k <
laell+1; k++) {
2358 for ( k =
kk; k < kmax+1; k++)
2426 Int_t apos,
iblk,
ifr,
ilvl,
ipiv,
ipos,
irhs,
irow,
ist,
j,
j1=0,
j2,
j3,
jj,
jpiv,k,
k1,
k2,
k3,
liell,
npiv;
2464 if (
jpiv == 1)
continue;
2556 Int_t apos,
apos2,
i1rhs,
i2rhs,
iblk,
ifr,
iipiv,
iirhs,
ilvl,
ipiv,
ipos,
irhs,
ist,
2557 j,
j1=0,
j2=0,
jj,
jj1,
jj2,
jpiv,
jpos=0,k,
liell,loop,
npiv;
2565 for (loop = 1; loop <
n+1; loop++) {
2568 if (
iblk < 1)
break;
2591 if (
jpiv == 1)
continue;
2695 fact.Print(
"fFact");
int Int_t
Signed integer 4 bytes (int)
double Double_t
Double 8 bytes.
const char Option_t
Option string (const char)
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
const Double_t kThresholdPivotingFactor
const Double_t kThresholdPivotingMax
const Double_t kInitThresholdPivoting
const Double_t kInitTreatAsZero
const Double_t kInitPrecision
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize id
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t UChar_t len
Array of doubles (64 bits per element).
void Set(Int_t n) override
Set size of this array to n doubles.
const Double_t * GetArray() const
Array of integers (32 bits per element).
void Set(Int_t n) override
Set size of this array to n ints.
const Int_t * GetArray() const
Decomposition Base class.
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
void Print(Option_t *opt="") const override
Print class members.
Sparse Symmetric Decomposition class.
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.
TDecompSparse()
Default constructor.
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.
static Int_t IDiag(Int_t ix, Int_t iy)
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_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 Int_t NonZerosUpperTriang(const TMatrixDSparse &a)
Static function, returning the number of non-zero entries in the upper triangular matrix .
void Print(Option_t *opt="") const override
Print class members.
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_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.
static void CopyUpperTriang(const TMatrixDSparse &a, Double_t *b)
Static function, copying the non-zero entries in the upper triangle to array b .
Double_t GetThresholdPivoting()
void SetTreatAsZero(Double_t tol)
TDecompSparse & operator=(const TDecompSparse &source)
Assignment operator.
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.
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.
virtual void SetMatrix(const TMatrixDSparse &a)
Set matrix to be decomposed .
void SetThresholdPivoting(Double_t piv)
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.
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.
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.
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_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.
void InitParam()
initializing control parameters
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 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.
Bool_t Decompose() override
Decomposition engine .
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)
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
virtual void ls(Option_t *option="") const
The ls function lists the contents of a class on stdout.
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Int_t Nint(T x)
Round to nearest integer. Rounds half integers to the nearest even integer.
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.