Logo ROOT   6.12/07
Reference Guide
TDecompSparse.cxx
Go to the documentation of this file.
1 // @(#)root/matrix:$Id$
2 // Authors: Fons Rademakers, Eddy Offermann Apr 2004
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #include "TDecompSparse.h"
13 #include "TMath.h"
14 
16 
17 /** \class TDecompSparse
18  \ingroup Matrix
19 
20  Sparse Symmetric Decomposition class
21 
22  Solve a sparse symmetric system of linear equations using a method
23  based on Gaussian elimination as discussed in Duff and Reid,
24  ACM Trans. Math. Software 9 (1983), 302-325.
25 */
26 
27 ////////////////////////////////////////////////////////////////////////////////
28 /// Default constructor
29 
31 {
32  fVerbose = 0;
33  InitParam();
34  memset(fInfo,0,21*sizeof(Int_t));
35 }
36 
37 ////////////////////////////////////////////////////////////////////////////////
38 /// Constructor for a matrix with nrows and unspecified number of columns .
39 /// nr_nonZeros is the total number of non-zero entries in the matrix .
40 
41 TDecompSparse::TDecompSparse(Int_t nRows,Int_t nr_nonZeros,Int_t verbose)
42 {
43  fVerbose = verbose;
44  InitParam();
45 
46  fNrows = nRows;
47  fNnonZeros = nr_nonZeros;
48 
49  fRowFact.Set(nr_nonZeros+1);
50  fColFact.Set(nr_nonZeros+1);
51  fW .Set(fNrows+1);
52  fIkeep .Set(3*(fNrows+1));
53  fIw .Set((Int_t)(1.3 * (2*fNnonZeros+3*fNrows+1)+1));
54  fIw1 .Set(2*(fNrows+1));
55 
56  memset(fInfo,0,21*sizeof(Int_t));
57 
58  // These parameters can only be set after sparsity/pivoting pattern is known
59  fNsteps = 0;
60  fMaxfrt = 0;
61 }
62 
63 ////////////////////////////////////////////////////////////////////////////////
64 /// Constructor for a matrix with row range, [row_lwb..row_upb] and unspecified column
65 /// range . nr_nonZeros is the total number of non-zero entries in the matrix .
66 
67 TDecompSparse::TDecompSparse(Int_t row_lwb,Int_t row_upb,Int_t nr_nonZeros,Int_t verbose)
68 {
69  fVerbose = verbose;
70  InitParam();
71 
72  fRowLwb = row_lwb;
73  fColLwb = row_lwb;
74  fNrows = row_upb-row_lwb+1;
75  fNnonZeros = nr_nonZeros;
76 
77  fRowFact.Set(nr_nonZeros+1);
78  fColFact.Set(nr_nonZeros+1);
79  fW .Set(fNrows+1);
80  fIkeep .Set(3*(fNrows+1));
81  fIw .Set((Int_t)(1.3 * (2*fNnonZeros+3*fNrows+1)+1));
82  fIw1 .Set(2*(fNrows+1));
83 
84  memset(fInfo,0,21*sizeof(Int_t));
85 
86  // These parameters can only be set after sparsity/pivoting pattern is known
87  fNsteps = 0;
88  fMaxfrt = 0;
89 }
90 
91 ////////////////////////////////////////////////////////////////////////////////
92 /// Constructor for matrix A .
93 
95 {
96  fVerbose = verbose;
97 
98  InitParam();
99  SetMatrix(a);
100 
101  memset(fInfo,0,21*sizeof(Int_t));
102 }
103 
104 ////////////////////////////////////////////////////////////////////////////////
105 /// Copy constructor
106 
108 {
109  *this = another;
110 }
111 
112 ////////////////////////////////////////////////////////////////////////////////
113 /// Static function, returning the number of non-zero entries in the upper triangular matrix .
114 
116 {
117  const Int_t rowLwb = a.GetRowLwb();
118  const Int_t colLwb = a.GetColLwb();
119  const Int_t nrows = a.GetNrows();;
120  const Int_t *pRowIndex = a.GetRowIndexArray();
121  const Int_t *pColIndex = a.GetColIndexArray();
122 
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++;
129  }
130  }
131 
132  return nr_nonzeros;
133 }
134 
135 ////////////////////////////////////////////////////////////////////////////////
136 /// Static function, copying the non-zero entries in the upper triangle to
137 /// array b . User should allocate enough memory for array b .
138 
140 {
141  const Int_t rowLwb = a.GetRowLwb();
142  const Int_t colLwb = a.GetColLwb();
143  const Int_t nrows = a.GetNrows();;
144  const Int_t *pRowIndex = a.GetRowIndexArray();
145  const Int_t *pColIndex = a.GetColIndexArray();
146  const Double_t *pData = a.GetMatrixArray();
147 
148  Int_t nr = 0;
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];
154  }
155  }
156 }
157 
158 ////////////////////////////////////////////////////////////////////////////////
159 /// Set matrix to be decomposed .
160 
162 {
163  ResetStatus();
164 
165  fA.Use(*const_cast<TMatrixDSparse *>(&a));
166  fRowLwb = fA.GetRowLwb();
167  fColLwb = fA.GetColLwb();
168  fNrows = fA.GetNrows();
170 
173 
174  const Int_t *rowIndex = a.GetRowIndexArray();
175  const Int_t *colIndex = a.GetColIndexArray();
176 
177  Int_t nr = 0;
178  for (Int_t irow = 0; irow < fNrows; irow++ ) {
179  const Int_t rown = irow+fRowLwb;
180  for (Int_t index = rowIndex[irow]; index < rowIndex[irow+1]; index++ ) {
181  const Int_t coln = colIndex[index]+fColLwb;
182  if (coln >= rown) {
183  fRowFact[nr+1] = irow+1;
184  fColFact[nr+1] = colIndex[index]+1;
185  nr++;
186  }
187  }
188  }
189 
190  fW .Set(fNrows+1);
191  fIkeep.Set(3*(fNrows+1));
192  fIw .Set((Int_t)(1.3 * (2*fNnonZeros+3*fNrows+1)+1));
193  fIw1 .Set(2*(fNrows+1));
194 
195  // Determine pivot sequence, set iflag = 0 in order to make InitPivot choose the order.
196  Int_t iflag = 0;
197  Double_t ops;
199  fIcntl,fCntl,fInfo,ops);
200 
201  switch ( this->ErrorFlag() ) {
202  case -1 :
203  Error("SetMatrix(const TMatrixDSparse &","nRows = %d out of range",fNrows);
204  return;
205  case -2 :
206  Error("SetMatrix(const TMatrixDSparse &","nr_nonzeros = %d out of range",fNnonZeros);
207  return;
208  case -3 :
209  Error("SetMatrix(const TMatrixDSparse &",
210  "insufficient space in fIw of %d suggest reset to %d",fIw.GetSize(),this->IError());
211  return;
212  case 1 :
213  Error("SetMatrix(const TMatrixDSparse &",
214  "detected %d entries out of rage in row/col indices; ignored",this->IError());
215  return;
216  }
217 
218  // set fIw and fIw1 in prep for calls to Factor and Solve
219 
220 // fIw .Set((Int_t) 1.2*this->MinRealWorkspace()+1);
221  fIw .Set((Int_t) 3*this->MinRealWorkspace()+1);
222  fIw1 .Set(fNrows+1);
223  fIw2 .Set(fNsteps+1);
224 // fFact.Set((Int_t) 1.2*this->MinRealWorkspace()+1);
225  fFact.Set((Int_t) 3*this->MinRealWorkspace()+1);
226 
228 }
229 
230 ////////////////////////////////////////////////////////////////////////////////
231 /// Decomposition engine .
232 /// If the decomposition succeeds, bit kDecomposed is set .
233 
235 {
236  if (TestBit(kDecomposed)) return kTRUE;
237 
238  if ( !TestBit(kMatrixSet) ) {
239  Error("Decompose()","Matrix has not been set");
240  return kFALSE;
241  }
242 
243  Int_t done = 0; Int_t tries = 0;
244  do {
245  fFact[0] = 0.;
247 
250 
251  switch ( this->ErrorFlag() ) {
252  case 0 :
253  done = 1;
254  break;
255  case -1 :
256  Error("Decompose()","nRows = %d out of range",fNrows);
257  return kFALSE;
258  case -2 :
259  Error("Decompose()","nr_nonzeros = %d out of range",fNnonZeros);
260  return kFALSE;
261  case -3 :
262  {
263  if (fVerbose)
264  Info("Decompose()","insufficient space of fIw: %d",fIw.GetSize());
265  const Int_t nIw_old = fIw.GetSize();
266  const Int_t nIw = (this->IError() > fIPessimism*nIw_old) ? this->IError() :
267  (Int_t)(fIPessimism*nIw_old);
268  fIw.Set(nIw);
269  if (fVerbose)
270  Info("Decompose()","resetting to fIw: %d",nIw);
271  fIPessimism *= 1.1;
272  break;
273  }
274  case -4 :
275  {
276  if (fVerbose)
277  Info("Decompose()","insufficient factorization space: %d",fFact.GetSize());
278  const Int_t nFact_old = fFact.GetSize();
279  const Int_t nFact = (this->IError() > fRPessimism*nFact_old) ? this->IError() :
280  (Int_t) (fRPessimism*nFact_old);
281  fFact.Set(nFact); fFact.Reset(0.0);
283  if (fVerbose)
284  Info("Decompose()","reseting to: %d",nFact);
285  fRPessimism *= 1.1;
286  break;
287  }
288  case -5 :
289  if (fVerbose) {
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..");
293  }
294  done = 1;
295  break;
296  case -6 :
297  if (fVerbose) {
298  Info("Decompose()","change of sign of pivots detected at stage %d",this->IError());
299  Info("Decompose()","but who cares ");
300  }
301  done = 1;
302  break;
303  case -7 :
304  Error("Decompose()","value of fNsteps out of range: %d",fNsteps);
305  return kFALSE;
306  case 1 :
307  if (fVerbose) {
308  Info("Decompose()","detected %d entries out of range in row/column index",this->IError());
309  Info("Decompose()","they are ignored");
310  }
311  done = 1;
312  break;
313  case 3 :
314  if (fVerbose)
315  Info("Decompose()","rank deficient matrix detected; apparent rank = %d",this->IError());
316  done = 1;
317  break;
318  default:
319  break;
320  }
321 
322  tries++;
323  } while (!done && tries < 10);
324 
325  Int_t ok;
326  if ( !done && tries >= 10) {
327  ok = kFALSE;
328  if (fVerbose)
329  Error("Decompose()","did not get a factorization after 10 tries");
330  } else {
331  ok = kTRUE;
333  }
334 
335  return ok;
336 }
337 
338 ////////////////////////////////////////////////////////////////////////////////
339 /// Solve Ax=b . Solution returned in b.
340 
342 {
343  R__ASSERT(b.IsValid());
344  if (TestBit(kSingular)) {
345  Error("Solve()","Matrix is singular");
346  return kFALSE;
347  }
348  if ( !TestBit(kDecomposed) ) {
349  if (!Decompose()) {
350  Error("Solve()","Decomposition failed");
351  return kFALSE;
352  }
353  }
354 
355  if (fNrows != b.GetNrows() || fRowLwb != b.GetLwb())
356  {
357  Error("Solve(TVectorD &","vector and matrix incompatible");
358  return kFALSE;
359  }
360  b.Shift(-fRowLwb); // make sure rowlwb = 0
361 
362  // save bs and store residuals
363  TVectorD resid = b;
364  TVectorD bSave = b;
365 
366  Double_t bnorm = b.NormInf();
367  Double_t rnorm = 0.0;
368 
369  Int_t done = 0;
370  Int_t refactorizations = 0;
371 
372  while (!done && refactorizations < 10) {
373 
375 
376  // compute residuals
377  resid = fA*b-resid;
378  rnorm = resid.NormInf();
379 
380  if (rnorm < fPrecision*(1.+bnorm)) {
381  // residuals are small enough, use this solution
382  done = 1;
383  } else if (this->GetThresholdPivoting() >= kThresholdPivotingMax
384  || refactorizations > 10) {
385  // ThresholdPivoting parameter is already too high; give up and
386  // use this solution, whatever it is (the algorithm may bomb in
387  // an outer loop).
388  done = 1;
389  } else {
390  // refactor with a higher Threshold Pivoting parameter
391  Double_t tp = this->GetThresholdPivoting();
394  this->SetThresholdPivoting(tp);
395  if (fVerbose)
396  Info("Solve","Setting ThresholdPivoting parameter to %.4e for future factorizations",
397  this->GetThresholdPivoting());
398 
399  SetMatrix(fA);
400  refactorizations++;
401  resid = bSave;
402  b = bSave;
403  }
404  }
405 
406  b.Shift(fRowLwb);
407  return kTRUE;
408 }
409 
410 ////////////////////////////////////////////////////////////////////////////////
411 /// initializing control parameters
412 
414 {
416  fIPessimism = 1.2;
417  fRPessimism = 1.2;
418 
419  const Int_t ifrlvl = 5;
420 
422  fIcntl[4] = 2139062143;
423  fIcntl[5] = 1;
424  fIcntl[ifrlvl+1] = 32639;
425  fIcntl[ifrlvl+2] = 32639;
426  fIcntl[ifrlvl+3] = 32639;
427  fIcntl[ifrlvl+4] = 32639;
428  fIcntl[ifrlvl+5] = 14;
429  fIcntl[ifrlvl+6] = 9;
430  fIcntl[ifrlvl+7] = 8;
431  fIcntl[ifrlvl+8] = 8;
432  fIcntl[ifrlvl+9] = 9;
433  fIcntl[ifrlvl+10] = 10;
434  fIcntl[ifrlvl+11] = 32639;
435  fIcntl[ifrlvl+12] = 32639;
436  fIcntl[ifrlvl+13] = 32639;
437  fIcntl[ifrlvl+14] = 32689;
438  fIcntl[ifrlvl+15] = 24;
439  fIcntl[ifrlvl+16] = 11;
440  fIcntl[ifrlvl+17] = 9;
441  fIcntl[ifrlvl+18] = 8;
442  fIcntl[ifrlvl+19] = 9;
443  fIcntl[ifrlvl+20] = 10;
444  fIcntl[26] = 0;
445  fIcntl[27] = 0;
446  fIcntl[28] = 0;
447  fIcntl[29] = 0;
448  fIcntl[30] = 0;
449  fCntl[1] = 0.10;
450  fCntl[2] = 1.00;
451  fCntl[3] = 0.00;
452  fCntl[4] = 0.0;
453  fCntl[5] = 0.0;
454 
455  // set initial value of "Treat As Zero" parameter
457 
458  // set initial value of Threshold parameter
460 
461  fNsteps = 0;
462  fMaxfrt = 0;
463  fNrows = 0;
464  fNnonZeros = 0;
465 }
466 
467 ////////////////////////////////////////////////////////////////////////////////
468 /// Setup Pivoting variables
469 
470 void TDecompSparse::InitPivot(const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,
471  TArrayI &Aiw,TArrayI &Aikeep,TArrayI &Aiw1,Int_t &nsteps,
472  const Int_t iflag,Int_t *icntl,Double_t *cntl,Int_t *info,
473  Double_t &ops)
474 {
475  Int_t i,iwfr,k,l1,l2,lliw;
476 
477  Int_t *irn = Airn.GetArray();
478  Int_t *icn = Aicn.GetArray();
479  Int_t *iw = Aiw.GetArray();
480  Int_t *ikeep = Aikeep.GetArray();
481  Int_t *iw1 = Aiw1.GetArray();
482  const Int_t liw = Aiw.GetSize()-1;
483 
484  for (i = 1; i < 16; i++)
485  info[i] = 0;
486 
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);
489  nsteps = 0;
490  k = TMath::Min(8,nz);
491  if (icntl[3] > 1) k = nz;
492  if (k > 0) {
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");
497  }
498  }
499 
500  k = TMath::Min(10,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");
506  }
507  }
508  }
509 
510  if (n >= 1 && n <= icntl[4]) {
511  if (nz < 0) {
512  info[1] = -2;
513  if (icntl[1] > 0)
514  ::Error("TDecompSparse::InitPivot","info[1]= %d; value of nz out of range .. = %d",info[1],nz);
515  return;
516  }
517  lliw = liw-2*n;
518  l1 = lliw+1;
519  l2 = l1+n;
520  if (iflag != 1) {
521  if (liw < 2*nz+3*n+1) {
522  info[1] = -3;
523  info[2] = 2*nz+3*n+1;
524  if (icntl[1] > 0)
525  ::Error("TDecompSparse::InitPivot","info[1]= %d; liw too small, must be increased from %d to at least %d",info[1],liw,info[2]);
526  return;
527  }
528  InitPivot_sub1(n,nz,irn,icn,iw,iw1,iw1+n+1,iw+l1-1,iwfr,icntl,info);
529  InitPivot_sub2(n,iw1,iw,lliw,iwfr,iw+l1-1,iw+l2-1,ikeep+n+1,
530  ikeep+2*(n+1),ikeep,icntl[4],info[11],cntl[2]);
531  } else {
532  if (liw < nz+3*n+1) {
533  info[1] = -3;
534  info[2] = nz+3*n+1;
535  if (icntl[1] > 0)
536  ::Error("TDecompSparse::InitPivot","info[1]= %d; liw too small, must be increased from %d to at least %d",info[1],liw,info[2]);
537  return;
538  }
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]);
541  }
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);
546  } else {
547  info[1] = -1;
548  if (icntl[1] > 0)
549  ::Error("TDecompSparse::InitPivot","info[1]= %d; value of n out of range ... = %d",info[1],n);
550  return;
551  }
552 
553  if (icntl[3] <= 0 || icntl[2] <= 0) return;
554 
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]);
558 
559  k = TMath::Min(9,n);
560  if (icntl[3] > 1) k = n;
561  if (k > 0) {
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");
566  }
567  }
568  k = TMath::Min(k,nsteps);
569  if (k > 0) {
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");
574  }
575  }
576 }
577 
578 ////////////////////////////////////////////////////////////////////////////////
579 /// Factorization routine, the workhorse for the decomposition step
580 
581 void TDecompSparse::Factor(const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,TArrayD &Aa,
582  TArrayI &Aiw,TArrayI &Aikeep,const Int_t nsteps,Int_t &maxfrt,
583  TArrayI &Aiw1,Int_t *icntl,Double_t *cntl,Int_t *info)
584 {
585  Int_t i,iapos,iblk,ipos,irows,j1,j2,jj,k,kblk,kz,len,ncols,nrows,nz1;
586 
587  Int_t *irn = Airn.GetArray();
588  Int_t *icn = Aicn.GetArray();
589  Int_t *iw = Aiw.GetArray();
590  Int_t *ikeep = Aikeep.GetArray();
591  Int_t *iw1 = Aiw1.GetArray();
592  Double_t *a = Aa.GetArray();
593 
594  const Int_t la = Aa.GetSize()-1;
595  const Int_t liw = Aiw.GetSize()-1;
596 
597  info[1] = 0;
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]);
601  kz = TMath::Min(6,nz);
602  if (icntl[3] > 1) kz = nz;
603  if (nz > 0) {
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");
608  }
609  }
610  k = TMath::Min(9,n);
611  if (icntl[3] > 1) k = n;
612  if (k > 0) {
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");
617  }
618  }
619  k = TMath::Min(k,nsteps);
620  if (k > 0) {
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");
625  }
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");
630  }
631  }
632  }
633 
634  if (n < 1 || n > icntl[4])
635  info[1] = -1;
636  else if (nz < 0)
637  info[1] = -2;
638  else if (liw < nz) {
639  info[1] = -3;
640  info[2] = nz;
641  } else if (la < nz+n) {
642  info[1] = -4;
643  info[2] = nz+n;
644  } else if (nsteps < 1 || nsteps > n)
645  info[1] = -7;
646  else {
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]);
652  }
653  }
654 
655  if (icntl[1] > 0) {
656  switch(info[1]) {
657  case -1:
658  ::Error("TDecompSparse::Factor","info[1]= %d; value of n out of range ... =%d",info[1],n);
659  break;
660 
661  case -2:
662  ::Error("TDecompSparse::Factor","info[1]= %d; value of nz out of range ... =%d",info[1],nz);
663  break;
664 
665  case -3:
666  ::Error("TDecompSparse::Factor","info[1]= %d; liw too small, must be increased from %d to at least %d",info[1],liw,info[2]);
667  break;
668 
669  case -4:
670  ::Error("TDecompSparse::Factor","info[1]= %d; la too small, must be increased from %d to at least %d",info[1],la,info[2]);
671  break;
672 
673  case -5:
674  ::Error("TDecompSparse::Factor","info[1]= %d; zero pivot at stage %d zero pivot at stage",info[1],info[2]);
675  break;
676 
677  case -6:
678  ::Error("TDecompSparse::Factor","info[1]= %d; change in sign of pivot encountered when factoring allegedly definite matrix",info[1]);
679  break;
680 
681  case -7:
682  ::Error("TDecompSparse::Factor","info[1]= %d; nsteps is out of range",info[1]);
683  break;
684  }
685  }
686 
687  if (icntl[3] <= 0 || icntl[2] <= 0 || info[1] < 0)
688  return;
689 
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]);
692 
693  if (info[1] < 0) return;
694 
695  kblk = TMath::Abs(iw[1]+0);
696  if (kblk == 0) return;
697  if (icntl[3] == 1) kblk = 1;
698  ipos = 2;
699  iapos = 1;
700 
701  for (iblk = 1; iblk < kblk+1; iblk++) {
702  ncols = iw[ipos];
703  nrows = iw[ipos+1];
704  j1 = ipos+2;
705  if (ncols <= 0) {
706  ncols = -ncols;
707  nrows = 1;
708  j1 = j1-1;
709  }
710  ::Info("TDecompSparse::Factor","block pivot =%d nrows =%d ncols =%d",iblk,nrows,ncols);
711  j2 = j1+ncols-1;
712  ipos = j2+1;
713 
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");
718  }
719 
720  printf(" real entries .. each row starts on a new line\n");
721  len = ncols;
722  for (irows = 1; irows < nrows+1; irows++) {
723  j1 = iapos;
724  j2 = iapos+len-1;
725  for (jj = j1; jj < j2+1; jj++) {
726  printf("%13.4e ",a[jj]);
727  if (jj%5 == 0 || jj == j2) printf("\n");
728  }
729  len = len-1;
730  iapos = j2+1;
731  }
732  }
733 }
734 
735 ////////////////////////////////////////////////////////////////////////////////
736 /// Main routine for solving Ax=b
737 
739  TArrayD &Aw,const Int_t maxfrt,TVectorD &b,TArrayI &Aiw1,
740  const Int_t nsteps,Int_t *icntl,Int_t *info)
741 {
742  Int_t i,iapos,iblk,ipos,irows,j1,j2,jj,k,kblk,latop,len,nblk,ncols,nrows;
743 
744  Double_t *a = Aa.GetArray();
745  Double_t *w = Aw.GetArray();
746  Int_t *iw = Aiw.GetArray();
747  Int_t *iw1 = Aiw1.GetArray();
748  Double_t *rhs = new Double_t[n+1];
749  rhs[0] = 0.;
750  memcpy(rhs+1,b.GetMatrixArray(),n*sizeof(Double_t));
751  const Int_t la = Aa.GetSize()-1;
752  const Int_t liw = Aiw.GetSize()-1;
753 
754  info[1] = 0;
755  k = 0;
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);
758 
759  kblk = TMath::Abs(iw[1]+0);
760  if (kblk != 0) {
761  if (icntl[3] == 1) kblk = 1;
762  ipos = 2;
763  iapos = 1;
764  for (iblk = 1; iblk < kblk+1; iblk++) {
765  ncols = iw[ipos];
766  nrows = iw[ipos+1];
767  j1 = ipos+2;
768  if (ncols <= 0) {
769  ncols = -ncols;
770  nrows = 1;
771  j1 = j1-1;
772  }
773  printf("block pivot=%d nrows=%d ncols=%d\n",iblk,nrows,ncols);
774  j2 = j1+ncols-1;
775  ipos = j2+1;
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");
780  }
781  printf("real entries .. each row starts on a new line\n");
782  len = ncols;
783  for (irows = 1; irows < nrows+1; irows++) {
784  j1 = iapos;
785  j2 = iapos+len-1;
786  for (jj = j1; jj < j2+1; jj++) {
787  printf("%13.3e ",a[jj]);
788  if (jj%5 == 0 || jj == j2) printf("\n");
789  }
790  len = len-1;
791  iapos = j2+1;
792  }
793  }
794  }
795 
796  k = TMath::Min(10,n);
797  if (icntl[3] > 1) k = n;
798  if (n > 0) {
799  printf("rhs =\n");
800  for (i = 1; i < k+1; i++) {
801  printf("%13.3e ",rhs[i]);
802  if (i%5 == 0 || i == k) printf("\n");
803  }
804  }
805  }
806 
807  nblk = 0;
808  if (iw[1] == 0) {
809  nblk = 0;
810  for (i = 1; i < n+1; i++)
811  rhs[i] = 0.0;
812  } else {
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);
816  }
817 
818  if (icntl[3] > 0 && icntl[2] > 0) {
819  printf("leaving Solve with:\n");
820  if (n > 0) {
821  printf("rhs =\n");
822  for (i = 1; i < k+1; i++) {
823  printf("%13.3e ",rhs[i]);
824  if (i%5 == 0 || i == k) printf("\n");
825  }
826  }
827  }
828 
829  memcpy(b.GetMatrixArray(),rhs+1,n*sizeof(Double_t));
830  delete [] rhs;
831 }
832 
833 ////////////////////////////////////////////////////////////////////////////////
834 /// Help routine for pivoting setup
835 
836 void TDecompSparse::InitPivot_sub1(const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,
837  Int_t *iw,Int_t *ipe,Int_t *iq,Int_t *flag,
838  Int_t &iwfr,Int_t *icntl,Int_t *info)
839 {
840  Int_t i,id,j,jn,k,k1,k2,l,last,lr,n1,ndup;
841 
842  info[2] = 0;
843  for (i = 1; i < n+1; i++)
844  ipe[i] = 0;
845  lr = nz;
846 
847  if (nz != 0) {
848  for (k = 1; k < nz+1; k++) {
849  i = irn[k];
850  j = icn[k];
851 
852  Bool_t outRange = (i < 1 || i > n || j < 1 || j > n);
853  if (outRange) {
854  info[2] = info[2]+1;
855  info[1] = 1;
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);
858  }
859 
860  if (outRange || i == j) {
861  i = 0;
862  j = 0;
863  } else {
864  ipe[i] = ipe[i]+1;
865  ipe[j] = ipe[j]+1;
866  }
867  iw[k] = j;
868  lr = lr+1;
869  iw[lr] = i;
870  }
871  }
872 
873  iq[1] = 1;
874  n1 = n-1;
875  if (n1 > 0) {
876  for (i = 1; i < n1+1; i++) {
877  flag[i] = 0;
878  if (ipe[i] == 0) ipe[i] = -1;
879  iq[i+1] = ipe[i]+iq[i]+1;
880  ipe[i] = iq[i];
881  }
882  }
883 
884  last = ipe[n]+iq[n];
885  flag[n] = 0;
886  if (lr < last) {
887  k1 = lr+1;
888  for (k = k1; k < last+1; k++)
889  iw[k] = 0;
890  }
891  ipe[n] = iq[n];
892  iwfr = last+1;
893  if (nz != 0) {
894  for (k = 1; k < nz+1; k++) {
895  j = iw[k];
896  if (j <= 0) continue;
897  l = k;
898  iw[k] = 0;
899  for (id = 1; id < nz+1; id++) {
900  if (l <= nz)
901  l = l+nz;
902  else
903  l = l-nz;
904  i = iw[l];
905  iw[l] = 0;
906  if (i >= j) {
907  l = iq[j]+1;
908  iq[j] = l;
909  jn = iw[l];
910  iw[l] = -i;
911  } else {
912  l = iq[i]+1;
913  iq[i] = l;
914  jn = iw[l];
915  iw[l] = -j;
916  }
917  j = jn;
918  if (j <= 0) break;
919  }
920  }
921  }
922 
923  ndup = 0;
924 
925  for (i = 1; i < n+1; i++) {
926  k1 = ipe[i]+1;
927  k2 = iq[i];
928  if (k1 > k2) {
929  ipe[i] = 0;
930  iq[i] = 0;
931  } else {
932  for (k = k1; k < k2+1; k++) {
933  j = -iw[k];
934  if (j <= 0) break;
935  l = iq[j]+1;
936  iq[j] = l;
937  iw[l] = i;
938  iw[k] = j;
939  if (flag[j] == i) {
940  ndup = ndup + 1;
941  iw[l] = 0;
942  iw[k] = 0;
943  }
944  flag[j] = i;
945  }
946 
947  iq[i] = iq[i]-ipe[i];
948  if (ndup == 0) iw[k1-1] = iq[i];
949  }
950  }
951 
952  if (ndup != 0) {
953  iwfr = 1;
954  for (i = 1; i < n+1; i++) {
955  k1 = ipe[i]+1;
956  if (k1 == 1) continue;
957  k2 = iq[i]+ipe[i];
958  l = iwfr;
959  ipe[i] = iwfr;
960  iwfr = iwfr+1;
961  for (k = k1; k < k2+1; k++) {
962  if (iw[k] == 0) continue;
963  iw[iwfr] = iw[k];
964  iwfr = iwfr+1;
965  }
966  iw[l] = iwfr-l-1;
967  }
968  }
969 
970 }
971 
972 ////////////////////////////////////////////////////////////////////////////////
973 /// Help routine for pivoting setup
974 
975 void TDecompSparse::InitPivot_sub2(const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,
976  Int_t &iwfr,Int_t *nv,Int_t *nxt,Int_t *lst,Int_t *ipd,
977  Int_t *flag,const Int_t iovflo,Int_t &ncmpa,
978  const Double_t fratio)
979 {
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;
983 
984  for (i = 1; i < n+1; i++) {
985  ipd[i] = 0;
986  nv[i] = 1;
987  flag[i] = iovflo;
988  }
989 
990  js = 0;
991  ms = 0;
992  ncmpa = 0;
993  md = 1;
994  nflg = iovflo;
995  nel = 0;
996  root = n+1;
997  nvroot = 0;
998  for (is = 1; is < n+1; is++) {
999  k = ipe[is];
1000  if (k > 0) {
1001  id = iw[k]+1;
1002  ns = ipd[id];
1003  if (ns > 0) lst[ns] = is;
1004  nxt[is] = ns;
1005  ipd[id] = is;
1006  lst[is] = -id;
1007  } else {
1008  nel = nel+1;
1009  flag[is] = -1;
1010  nxt[is] = 0;
1011  lst[is] = 0;
1012  }
1013  }
1014 
1015  for (ml = 1; ml < n+1; ml++) {
1016  if (nel+nvroot+1 >= n) break;
1017  for (id = md; id < n+1; id++) {
1018  ms = ipd[id];
1019  if (ms > 0) break;
1020  }
1021 
1022  md = id;
1023  nvpiv = nv[ms];
1024  ns = nxt[ms];
1025  nxt[ms] = 0;
1026  lst[ms] = 0;
1027  if (ns > 0) lst[ns] = -id;
1028  ipd[id] = ns;
1029  me = ms;
1030  nel = nel+nvpiv;
1031  idn = 0;
1032  kp = ipe[me];
1033  flag[ms] = -1;
1034  ip = iwfr;
1035  len = iw[kp];
1036  jp = 0;
1037  for (kp1 = 1; kp1 < len+1; kp1++) {
1038  kp = kp+1;
1039  ke = iw[kp];
1040  if (flag[ke] > -2) {
1041  if (flag[ke] <= 0) {
1042  if (ipe[ke] != -root) continue;
1043  ke = root;
1044  if (flag[ke] <= 0) continue;
1045  }
1046  jp = kp-1;
1047  ln = len-kp1+1;
1048  ie = ms;
1049  } else {
1050  ie = ke;
1051  jp = ipe[ie];
1052  ln = iw[jp];
1053  }
1054 
1055  for (jp1 = 1; jp1 < ln+1; jp1++) {
1056  jp = jp+1;
1057  is = iw[jp];
1058  if (flag[is] <= 0) {
1059  if (ipe[is] == -root) {
1060  is = root;
1061  iw[jp] = root;
1062  if (flag[is] <= 0) continue;
1063  } else
1064  continue;
1065  }
1066  flag[is] = 0;
1067  if (iwfr >= lw) {
1068  ipe[ms] = kp;
1069  iw[kp] = len-kp1;
1070  ipe[ie] = jp;
1071  iw[jp] = ln-jp1;
1072  InitPivot_sub2a(n,ipe,iw,ip-1,lwfr,ncmpa);
1073  jp2 = iwfr-1;
1074  iwfr = lwfr;
1075  if (ip <= jp2) {
1076  for (jp = ip; jp < jp2+1; jp++) {
1077  iw[iwfr] = iw[jp];
1078  iwfr = iwfr+1;
1079  }
1080  }
1081  ip = lwfr;
1082  jp = ipe[ie];
1083  kp = ipe[me];
1084  }
1085  iw[iwfr] = is;
1086  idn = idn+nv[is];
1087  iwfr = iwfr+1;
1088  ls = lst[is];
1089  lst[is] = 0;
1090  ns = nxt[is];
1091  nxt[is] = 0;
1092  if (ns > 0) lst[ns] = ls;
1093  if (ls < 0) {
1094  ls = -ls;
1095  ipd[ls] = ns;
1096  } else if (ls > 0)
1097  nxt[ls] = ns;
1098  }
1099 
1100  if (ie == ms)
1101  break;
1102  ipe[ie] = -me;
1103  flag[ie] = -1;
1104  }
1105 
1106  nv[ms] = idn+nvpiv;
1107  if (iwfr != ip) {
1108  k1 = ip;
1109  k2 = iwfr-1;
1110  limit = TMath::Nint(fratio*(n-nel));
1111 
1112  for (k = k1; k < k2+1; k++) {
1113  is = iw[k];
1114  if (is == root) continue;
1115  if (nflg <= 2) {
1116  for (i = 1; i < n+1; i++) {
1117  if (flag[i] > 0) flag[i] = iovflo;
1118  if (flag[i] <= -2) flag[i] = -iovflo;
1119  }
1120  nflg = iovflo;
1121  }
1122  nflg = nflg-1;
1123  id = idn;
1124  kp1 = ipe[is]+1;
1125  np = kp1;
1126  kp2 = iw[kp1-1]+kp1-1;
1127 
1128  Int_t skip = 0;
1129  for (kp = kp1; kp < kp2+1; kp++) {
1130  ke = iw[kp];
1131  if (flag[ke] == -1) {
1132  if (ipe[ke] != -root) continue;
1133  ke = root;
1134  iw[kp] = root;
1135  if (flag[ke] == -1) continue;
1136  }
1137  if (flag[ke] >= 0) {
1138  skip = 1;
1139  break;
1140  }
1141  jp1 = ipe[ke]+1;
1142  jp2 = iw[jp1-1]+jp1-1;
1143  idl = id;
1144  for (jp = jp1; jp < jp2+1; jp++) {
1145  js = iw[jp];
1146  if (flag[js] <= nflg) continue;
1147  id = id+nv[js];
1148  flag[js] = nflg;
1149  }
1150  if (id <= idl) {
1151  Int_t skip2 = 0;
1152  for (jp = jp1; jp < jp2+1; jp++) {
1153  js = iw[jp];
1154  if (flag[js] != 0) {
1155  skip2 = 1;
1156  break;
1157  }
1158  }
1159  if (skip2) {
1160  iw[np] = ke;
1161  flag[ke] = -nflg;
1162  np = np+1;
1163  } else {
1164  ipe[ke] = -me;
1165  flag[ke] = -1;
1166  }
1167  } else {
1168  iw[np] = ke;
1169  flag[ke] = -nflg;
1170  np = np+1;
1171  }
1172  }
1173 
1174  if (!skip)
1175  np0 = np;
1176  else {
1177  np0 = np;
1178  kp0 = kp;
1179  for (kp = kp0; kp < kp2+1; kp++) {
1180  ks = iw[kp];
1181  if (flag[ks] <= nflg) {
1182  if (ipe[ks] == -root) {
1183  ks = root;
1184  iw[kp] = root;
1185  if (flag[ks] <= nflg) continue;
1186  } else
1187  continue;
1188  }
1189  id = id+nv[ks];
1190  flag[ks] = nflg;
1191  iw[np] = ks;
1192  np = np+1;
1193  }
1194  }
1195 
1196  Int_t doit = 2;
1197  if (id < limit) {
1198  iw[np] = iw[np0];
1199  iw[np0] = iw[kp1];
1200  iw[kp1] = me;
1201  iw[kp1-1] = np-kp1+1;
1202  js = ipd[id];
1203  for (l = 1; l < n+1; l++) {
1204  if (js <= 0) {
1205  doit = 3;
1206  break;
1207  }
1208  kp1 = ipe[js]+1;
1209  if (iw[kp1] != me) {
1210  doit = 3;
1211  break;
1212  }
1213  kp2 = kp1-1+iw[kp1-1];
1214  Int_t stayInLoop = 0;
1215  for (kp = kp1; kp < kp2+1; kp++) {
1216  ie = iw[kp];
1217  if (TMath::Abs(flag[ie]+0) > nflg) {
1218  stayInLoop = 1;
1219  break;
1220  }
1221  }
1222  if (!stayInLoop) {
1223  doit = 1;
1224  break;
1225  }
1226  js = nxt[js];
1227  }
1228  }
1229 
1230  if (doit == 1) {
1231  ipe[js] = -is;
1232  nv[is] = nv[is]+nv[js];
1233  nv[js] = 0;
1234  flag[js] = -1;
1235  ns = nxt[js];
1236  ls = lst[js];
1237  if (ns > 0) lst[ns] = is;
1238  if (ls > 0) nxt[ls] = is;
1239  lst[is] = ls;
1240  nxt[is] = ns;
1241  lst[js] = 0;
1242  nxt[js] = 0;
1243  if (ipd[id] == js) ipd[id] = is;
1244  } else if (doit == 2) {
1245  if (nvroot == 0) {
1246  root = is;
1247  ipe[is] = 0;
1248  } else {
1249  iw[k] = root;
1250  ipe[is] = -root;
1251  nv[root] = nv[root]+nv[is];
1252  nv[is] = 0;
1253  flag[is] = -1;
1254  }
1255  nvroot = nv[root];
1256  } else if (doit == 3) {
1257  ns = ipd[id];
1258  if (ns > 0) lst[ns] = is;
1259  nxt[is] = ns;
1260  ipd[id] = is;
1261  lst[is] = -id;
1262  md = TMath::Min(md,id);
1263  }
1264  }
1265 
1266  for (k = k1; k < k2+1; k++) {
1267  is = iw[k];
1268  if (nv[is] == 0) continue;
1269  flag[is] = nflg;
1270  iw[ip] = is;
1271  ip = ip+1;
1272  }
1273  iwfr = k1;
1274  flag[me] = -nflg;
1275  iw[ip] = iw[k1];
1276  iw[k1] = ip-k1;
1277  ipe[me] = k1;
1278  iwfr = ip+1;
1279  } else
1280  ipe[me] = 0;
1281  }
1282 
1283  for (is = 1; is < n+1; is++) {
1284  if (nxt[is] != 0 || lst[is] != 0) {
1285  if (nvroot == 0) {
1286  root = is;
1287  ipe[is] = 0;
1288  } else {
1289  ipe[is] = -root;
1290  }
1291  nvroot = nvroot+nv[is];
1292  nv[is] = 0;
1293  }
1294  }
1295 
1296  for (ie = 1; ie < n+1; ie++)
1297  if (ipe[ie] > 0) ipe[ie] = -root;
1298 
1299  if (nvroot> 0) nv[root] = nvroot;
1300 }
1301 
1302 ////////////////////////////////////////////////////////////////////////////////
1303 /// Help routine for pivoting setup
1304 
1306  Int_t &iwfr,Int_t &ncmpa)
1307 {
1308  Int_t i,ir,k,k1,k2,lwfr;
1309 
1310  ncmpa = ncmpa+1;
1311  for (i = 1; i < n+1; i++) {
1312  k1 = ipe[i];
1313  if (k1 <= 0) continue;
1314  ipe[i] = iw[k1];
1315  iw[k1] = -i;
1316  }
1317 
1318  iwfr = 1;
1319  lwfr = iwfr;
1320  for (ir = 1; ir < n+1; ir++) {
1321  if (lwfr > lw) break;
1322  Int_t skip = 1;
1323  for (k = lwfr; k < lw+1; k++) {
1324  if (iw[k] < 0) {
1325  skip = 0;
1326  break;
1327  }
1328  }
1329  if (skip) break;
1330  i = -iw[k];
1331  iw[iwfr] = ipe[i];
1332  ipe[i] = iwfr;
1333  k1 = k+1;
1334  k2 = k+iw[iwfr];
1335  iwfr = iwfr+1;
1336  if (k1 <= k2) {
1337  for (k = k1; k < k2+1; k++) {
1338  iw[iwfr] = iw[k];
1339  iwfr = iwfr+1;
1340  }
1341  }
1342  lwfr = k2+1;
1343  }
1344 }
1345 
1346 ////////////////////////////////////////////////////////////////////////////////
1347 /// Help routine for pivoting setup
1348 
1349 void TDecompSparse::InitPivot_sub3(const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,
1350  Int_t *perm,Int_t *iw,Int_t *ipe,Int_t *iq,
1351  Int_t *flag,Int_t &iwfr,Int_t *icntl,Int_t *info)
1352 {
1353  Int_t i,id,in,j,jdummy,k,k1,k2,l,lbig,len;
1354 
1355  info[1] = 0;
1356  info[2] = 0;
1357  for (i = 1; i < n+1; i++)
1358  iq[i] = 0;
1359 
1360  if (nz != 0) {
1361  for (k = 1; k < nz+1; k++) {
1362  i = irn[k];
1363  j = icn[k];
1364  iw[k] = -i;
1365 
1366  Bool_t outRange = (i < 1 || i > n || j < 1 || j > n);
1367  if (outRange) {
1368  info[2] = info[2]+1;
1369  info[1] = 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);
1372  }
1373 
1374  if (outRange || i==j) {
1375  iw[k] = 0;
1376  } else {
1377  if (perm[j] <= perm[i])
1378  iq[j] = iq[j]+1;
1379  else
1380  iq[i] = iq[i]+1;
1381  }
1382  }
1383  }
1384 
1385  iwfr = 1;
1386  lbig = 0;
1387  for (i = 1; i < n+1; i++) {
1388  l = iq[i];
1389  lbig = TMath::Max(l,lbig);
1390  iwfr = iwfr+l;
1391  ipe[i] = iwfr-1;
1392  }
1393 
1394  if (nz != 0) {
1395  for (k = 1; k < nz+1; k++) {
1396  i = -iw[k];
1397  if (i <= 0) continue;
1398  l = k;
1399  iw[k] = 0;
1400  for (id = 1; id < nz+1; id++) {
1401  j = icn[l];
1402  if (perm[i] >= perm[j]) {
1403  l = ipe[j];
1404  ipe[j] = l-1;
1405  in = iw[l];
1406  iw[l] = i;
1407  } else {
1408  l = ipe[i];
1409  ipe[i] = l-1;
1410  in = iw[l];
1411  iw[l] = j;
1412  }
1413  i = -in;
1414  if (i <= 0) continue;
1415  }
1416  }
1417 
1418  k = iwfr-1;
1419  l = k+n;
1420  iwfr = l+1;
1421  for (i = 1; i < n+1; i++) {
1422  flag[i] = 0;
1423  j = n+1-i;
1424  len = iq[j];
1425  if (len > 0) {
1426  for (jdummy = 1; jdummy < len+1; jdummy++) {
1427  iw[l] = iw[k];
1428  k = k-1;
1429  l = l-1;
1430  }
1431  }
1432  ipe[j] = l;
1433  l = l-1;
1434  }
1435 
1436  if (lbig < icntl[4]) {
1437  for (i = 1; i < n+1; i++) {
1438  k = ipe[i];
1439  iw[k] = iq[i];
1440  if (iq[i] == 0) ipe[i] = 0;
1441  }
1442  } else {
1443  iwfr = 1;
1444  for (i = 1; i < n+1; i++) {
1445  k1 = ipe[i]+1;
1446  k2 = ipe[i]+iq[i];
1447  if (k1 > k2) {
1448  ipe[i] = 0;
1449  } else {
1450  ipe[i] = iwfr;
1451  iwfr = iwfr+1;
1452  for (k = k1; k < k2+1; k++) {
1453  j = iw[k];
1454  if (flag[j] == i) continue;
1455  iw[iwfr] = j;
1456  iwfr = iwfr+1;
1457  flag[j] = i;
1458  }
1459  k = ipe[i];
1460  iw[k] = iwfr-k-1;
1461  }
1462  }
1463  }
1464  }
1465 
1466 }
1467 
1468 ////////////////////////////////////////////////////////////////////////////////
1469 /// Help routine for pivoting setup
1470 
1471 void TDecompSparse::InitPivot_sub4(const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,
1472  Int_t &iwfr,Int_t *ips,Int_t *ipv,Int_t *nv,Int_t *flag,
1473  Int_t &ncmpa)
1474 {
1475  Int_t i,ie,ip,j,je,jp,jp1,jp2,js,kdummy,ln,lwfr,me,minjs,ml,ms;
1476 
1477  for (i = 1; i < n+1; i++) {
1478  flag[i] = 0;
1479  nv[i] = 0;
1480  j = ips[i];
1481  ipv[j] = i;
1482  }
1483 
1484  ncmpa = 0;
1485  for (ml = 1; ml < n+1; ml++) {
1486  ms = ipv[ml];
1487  me = ms;
1488  flag[ms] = me;
1489  ip = iwfr;
1490  minjs = n;
1491  ie = me;
1492 
1493  for (kdummy = 1; kdummy < n+1; kdummy++) {
1494  jp = ipe[ie];
1495  ln = 0;
1496  if (jp > 0) {
1497  ln = iw[jp];
1498  for (jp1 = 1; jp1 < ln+1; jp1++) {
1499  jp = jp+1;
1500  js = iw[jp];
1501  if (flag[js] == me) continue;
1502  flag[js] = me;
1503  if (iwfr >= lw) {
1504  ipe[ie] = jp;
1505  iw[jp] = ln-jp1;
1506  InitPivot_sub2a(n,ipe,iw,ip-1,lwfr,ncmpa);
1507  jp2 = iwfr-1;
1508  iwfr = lwfr;
1509  if (ip <= jp2) {
1510  for (jp = ip; jp < jp2+1; jp++) {
1511  iw[iwfr] = iw[jp];
1512  iwfr = iwfr+1;
1513  }
1514  }
1515  ip = lwfr;
1516  jp = ipe[ie];
1517  }
1518  iw[iwfr] = js;
1519  minjs = TMath::Min(minjs,ips[js]+0);
1520  iwfr = iwfr+1;
1521  }
1522  }
1523  ipe[ie] = -me;
1524  je = nv[ie];
1525  nv[ie] = ln+1;
1526  ie = je;
1527  if (ie == 0) break;
1528  }
1529 
1530  if (iwfr <= ip) {
1531  ipe[me] = 0;
1532  nv[me] = 1;
1533  } else {
1534  minjs = ipv[minjs];
1535  nv[me] = nv[minjs];
1536  nv[minjs] = me;
1537  iw[iwfr] = iw[ip];
1538  iw[ip] = iwfr-ip;
1539  ipe[me] = ip;
1540  iwfr = iwfr+1;
1541  }
1542  }
1543 }
1544 
1545 ////////////////////////////////////////////////////////////////////////////////
1546 /// Help routine for pivoting setup
1547 
1549  Int_t *na,Int_t *nd,Int_t &nsteps,const Int_t nemin)
1550 {
1551  Int_t i,ib,iff,il,is,ison,k,l,nr;
1552 
1553  il = 0;
1554  for (i = 1; i < n+1; i++) {
1555  ips[i] = 0;
1556  ne[i] = 0;
1557  }
1558  for (i = 1; i < n+1; i++) {
1559  if (nv[i] > 0) continue;
1560  iff = -ipe[i];
1561  is = -ips[iff];
1562  if (is > 0) ipe[i] = is;
1563  ips[iff] = -i;
1564  }
1565 
1566  nr = n+1;
1567  for (i = 1; i < n+1; i++) {
1568  if (nv[i] <= 0) continue;
1569  iff = -ipe[i];
1570  if (iff != 0) {
1571  is = -ips[iff];
1572  if (is > 0)
1573  ipe[i] = is;
1574  ips[iff] = -i;
1575  } else {
1576  nr = nr-1;
1577  ne[nr] = i;
1578  }
1579  }
1580 
1581  is = 1;
1582  i = 0;
1583  for (k = 1; k < n+1; k++) {
1584  if (i <= 0) {
1585  i = ne[nr];
1586  ne[nr] = 0;
1587  nr = nr+1;
1588  il = n;
1589  na[n] = 0;
1590  }
1591  for (l = 1; l < n+1; l++) {
1592  if (ips[i] >= 0) break;
1593  ison = -ips[i];
1594  ips[i] = 0;
1595  i = ison;
1596  il = il-1;
1597  na[il] = 0;
1598  }
1599 
1600  ips[i] = k;
1601  ne[is] = ne[is]+1;
1602  if (nv[i] > 0) {
1603  if (il < n) na[il+1] = na[il+1]+1;
1604  na[is] = na[il];
1605  nd[is] = nv[i];
1606 
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);
1609 
1610  if (doit) {
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];
1614  ne[is] = 0;
1615  } else {
1616  is = is+1;
1617  }
1618  }
1619 
1620  ib = ipe[i];
1621  if (ib >= 0) {
1622  if (ib > 0)
1623  na[il] = 0;
1624  i = ib;
1625  } else {
1626  i = -ib;
1627  il = il+1;
1628  }
1629  }
1630 
1631  nsteps = is-1;
1632 }
1633 
1634 ////////////////////////////////////////////////////////////////////////////////
1635 /// Help routine for pivoting setup
1636 
1637 void TDecompSparse::InitPivot_sub6(const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,
1638  Int_t *perm,Int_t *na,Int_t *ne,Int_t *nd,const Int_t nsteps,
1639  Int_t *lstki,Int_t *lstkr,Int_t *iw,Int_t *info,Double_t &ops)
1640 {
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;
1643  Double_t delim;
1644 
1645  if (nz != 0 && irn[1] == iw[1]) {
1646  irn[1] = iw[1]-1;
1647  nz2 = 0;
1648  for (iold = 1; iold < n+1; iold++) {
1649  inew = perm[iold];
1650  lstki[inew] = lstkr[iold]+1;
1651  nz2 = nz2+lstkr[iold];
1652  }
1653  nz1 = nz2/2+n;
1654  nz2 = nz2+n;
1655  } else {
1656  for (i = 1; i < n+1; i++)
1657  lstki[i] = 1;
1658  nz1 = n;
1659  if (nz != 0) {
1660  for (i = 1; i < nz+1; i++) {
1661  iold = irn[i];
1662  jold = icn[i];
1663  if (iold < 1 || iold > n) continue;
1664  if (jold < 1 || jold > n) continue;
1665  if (iold == jold) continue;
1666  nz1 = nz1+1;
1667  irow = TMath::Min(perm[iold]+0,perm[jold]+0);
1668  lstki[irow] = lstki[irow]+1;
1669  }
1670  }
1671  nz2 = nz1;
1672  }
1673 
1674  ops = 0.0;
1675  istki = 0;
1676  istkr = 0;
1677  nrladu = 0;
1678  niradu = 1;
1679  nirtot = nz1;
1680  nrltot = nz1;
1681  nirnec = nz2;
1682  nrlnec = nz2;
1683  numorg = 0;
1684  itop = 0;
1685  for (itree = 1; itree < nsteps+1; itree++) {
1686  nelim = ne[itree];
1687  delim = Double_t(nelim);
1688  nfr = nd[itree];
1689  nstk = na[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++) {
1697  jorg = numorg+iorg;
1698  nz2 = nz2-lstki[jorg];
1699  }
1700  numorg = numorg+nelim;
1701  if (nstk > 0) {
1702  for (k = 1; k < nstk+1; k++) {
1703  lstk = lstkr[itop];
1704  istkr = istkr-lstk;
1705  lstk = lstki[itop];
1706  istki = istki-lstk;
1707  itop = itop-1;
1708  }
1709  }
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;
1715  itop = itop+1;
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);
1722  }
1723 
1724  nrlnec = TMath::Max(nrlnec,n+TMath::Max(nz,nz1));
1725  nrltot = TMath::Max(nrltot,n+TMath::Max(nz,nz1));
1726  nrlnec = TMath::Min(nrlnec,nrltot);
1727  nirnec = TMath::Max(nz,nirnec);
1728  nirtot = TMath::Max(nz,nirtot);
1729  nirnec = TMath::Min(nirnec,nirtot);
1730  info[3] = nrltot;
1731  info[4] = nirtot;
1732  info[5] = nrlnec;
1733  info[6] = nirnec;
1734  info[7] = nrladu;
1735  info[8] = niradu;
1736 }
1737 
1738 ////////////////////////////////////////////////////////////////////////////////
1739 /// Help routine for factorization
1740 
1742  const Int_t la,Int_t *irn,Int_t *icn,Int_t *iw,const Int_t liw,
1743  Int_t *perm,Int_t *iw2,Int_t *icntl,Int_t *info)
1744 {
1745  Int_t i,ia,ich,ii,iiw,inew,iold,ipos,j1,j2,jj,jnew,jold,jpos,k;
1746  Double_t anext,anow;
1747 
1748  const Double_t zero = 0.0;
1749  info[1] = 0;
1750  ia = la;
1751  for (iold = 1; iold < n+1; iold++) {
1752  iw2[iold] = 1;
1753  a[ia] = zero;
1754  ia = ia-1;
1755  }
1756 
1757  info[2] = 0;
1758  nz1 = n;
1759  if (nz != 0) {
1760  for (k = 1; k < nz+1; k++) {
1761  iold = irn[k];
1762  jold = icn[k];
1763  Bool_t outRange = (iold < 1 || iold > n || jold < 1 || jold > n);
1764 
1765  inew = perm[iold];
1766  jnew = perm[jold];
1767 
1768  if (!outRange && inew == jnew) {
1769  ia = la-n+iold;
1770  a[ia] = a[ia]+a[k];
1771  iw[k] = 0;
1772  } else {
1773  if (!outRange) {
1774  inew = TMath::Min(inew,jnew);
1775  iw2[inew] = iw2[inew]+1;
1776  iw[k] = -iold;
1777  nz1 = nz1+1;
1778  } else {
1779  info[1] = 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]);
1784  iw[k] = 0;
1785  }
1786  }
1787  }
1788  }
1789 
1790  if (nz >= nz1 || nz1 == n) {
1791  k = 1;
1792  for (i = 1; i < n+1; i++) {
1793  k = k+iw2[i];
1794  iw2[i] = k;
1795  }
1796  } else {
1797  k = 1;
1798  for (i = 1; i < n+1; i++) {
1799  k = k+iw2[i]-1;
1800  iw2[i] = k;
1801  }
1802  }
1803 
1804  if (nz1 > liw) {
1805  info[1] = -3;
1806  info[2] = nz1;
1807  return;
1808  }
1809 
1810  if (nz1+n > la) {
1811  info[1] = -4;
1812  info[2] = nz1+n;
1813  return;
1814  }
1815 
1816  if (nz1 != n) {
1817  for (k = 1; k < nz+1; k++) {
1818  iold = -iw[k];
1819  if (iold <= 0) continue;
1820  jold = icn[k];
1821  anow = a[k];
1822  iw[k] = 0;
1823  for (ich = 1; ich < nz+1; ich++) {
1824  inew = perm[iold];
1825  jnew = perm[jold];
1826  inew = TMath::Min(inew,jnew);
1827  if (inew == perm[jold]) jold = iold;
1828  jpos = iw2[inew]-1;
1829  iold = -iw[jpos];
1830  anext = a[jpos];
1831  a[jpos] = anow;
1832  iw[jpos] = jold;
1833  iw2[inew] = jpos;
1834  if (iold == 0) break;
1835  anow = anext;
1836  jold = icn[jpos];
1837  }
1838  }
1839 
1840  if (nz < nz1) {
1841  ipos = nz1;
1842  jpos = nz1-n;
1843  for (ii = 1; ii < n+1; ii++) {
1844  i = n-ii+1;
1845  j1 = iw2[i];
1846  j2 = jpos;
1847  if (j1 <= jpos) {
1848  for (jj = j1; jj < j2+1; jj++) {
1849  iw[ipos] = iw[jpos];
1850  a[ipos] = a[jpos];
1851  ipos = ipos-1;
1852  jpos = jpos-1;
1853  }
1854  }
1855  iw2[i] = ipos+1;
1856  ipos = ipos-1;
1857  }
1858  }
1859  }
1860 
1861  for (iold = 1; iold < n+1; iold++) {
1862  inew = perm[iold];
1863  jpos = iw2[inew]-1;
1864  ia = la-n+iold;
1865  a[jpos] = a[ia];
1866  iw[jpos] = -iold;
1867  }
1868  ipos = nz1;
1869  ia = la;
1870  iiw = liw;
1871  for (i = 1; i < nz1+1; i++) {
1872  a[ia] = a[ipos];
1873  iw[iiw] = iw[ipos];
1874  ipos = ipos-1;
1875  ia = ia-1;
1876  iiw = iiw-1;
1877  }
1878 }
1879 
1880 ////////////////////////////////////////////////////////////////////////////////
1881 /// Help routine for factorization
1882 
1883 void TDecompSparse::Factor_sub2(const Int_t n,const Int_t nz,Double_t *a,const Int_t la,
1884  Int_t *iw,const Int_t liw,Int_t *perm,Int_t *nstk,
1885  const Int_t nsteps,Int_t &maxfrt,Int_t *nelim,Int_t *iw2,
1886  Int_t *icntl,Double_t *cntl,Int_t *info)
1887 {
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;
1897 
1898  const Double_t zero = 0.0;
1899  const Double_t half = 0.5;
1900  const Double_t one = 1.0;
1901 
1902  detpiv = 0.0;
1903  apos3 = 0;
1904  isnpiv = 0;
1905  jmax = 0;
1906  jpos = 0;
1907 
1908  nblk = 0;
1909  ntwo = 0;
1910  neig = 0;
1911  ncmpbi = 0;
1912  ncmpbr = 0;
1913  maxfrt = 0;
1914  nrlbdu = 0;
1915  nirbdu = 0;
1916  uu = TMath::Min(cntl[1],half);
1917  uu = TMath::Max(uu,-half);
1918  for (i = 1; i < n+1; i++)
1919  iw2[i] = 0;
1920  iwpos = 2;
1921  posfac = 1;
1922  istk = liw-nz+1;
1923  istk2 = istk-1;
1924  astk = la-nz+1;
1925  astk2 = astk-1;
1926  iinput = istk;
1927  ainput = astk;
1928  azero = 0;
1929  ntotpv = 0;
1930  numass = 0;
1931 
1932  for (iass = 1; iass < nsteps+1; iass++) {
1933  nass = nelim[iass];
1934  newel = iwpos+1;
1935  jfirst = n+1;
1936  nfront = 0;
1937  numstk = nstk[iass];
1938  ltopst = 1;
1939  lnass = 0;
1940  if (numstk != 0) {
1941  j2 = istk-1;
1942  lnass = nass;
1943  ltopst = ((iw[istk]+1)*iw[istk])/2;
1944  for (iell = 1; iell < numstk+1; iell++) {
1945  jnext = jfirst;
1946  jlast = n+1;
1947  j1 = j2+2;
1948  j2 = j1-1+iw[j1-1];
1949  for (jj = j1; jj < j2+1; jj++) {
1950  j = iw[jj];
1951  if (iw2[j] > 0) continue;
1952  jnew = perm[j];
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;
1957  jlast = jnext;
1958  jnext = iw2[jlast];
1959  }
1960  if (jlast == n+1)
1961  jfirst = j;
1962  else
1963  iw2[jlast] = j;
1964  iw2[j] = jnext;
1965  jlast = j;
1966  nfront = nfront+1;
1967  }
1968  }
1969  lnass = nass-lnass;
1970  }
1971 
1972  numorg = nelim[iass];
1973  j1 = iinput;
1974  for (iorg = 1; iorg < numorg+1; iorg++) {
1975  j = -iw[j1];
1976  for (idummy = 1; idummy < liw+1; idummy++) {
1977  jnew = perm[j];
1978  if (iw2[j] <= 0) {
1979  jlast = n+1;
1980  jnext = jfirst;
1981  for (jdummy = 1; jdummy < n+1; jdummy++) {
1982  if (jnext == n+1) break;
1983  if (perm[jnext] > jnew) break;
1984  jlast = jnext;
1985  jnext = iw2[jlast];
1986  }
1987  if (jlast == n+1)
1988  jfirst = j;
1989  else
1990  iw2[jlast] = j;
1991  iw2[j] = jnext;
1992  nfront = nfront+1;
1993  }
1994  j1 = j1+1;
1995  if (j1 > liw) break;
1996  j = iw[j1];
1997  if (j < 0) break;
1998  }
1999  }
2000 
2001  if (newel+nfront >= istk)
2002  Factor_sub3(a,iw,istk,istk2,iinput,2,ncmpbr,ncmpbi);
2003  if (newel+nfront >= istk) {
2004  info[1] = -3;
2005  info[2] = liw+1+newel+nfront-istk;
2006  goto finish;
2007  }
2008 
2009  j = jfirst;
2010  for (ifr = 1; ifr < nfront+1; ifr++) {
2011  newel = newel+1;
2012  iw[newel] = j;
2013  jnext = iw2[j];
2014  iw2[j] = newel-(iwpos+1);
2015  j = jnext;
2016  }
2017 
2018  maxfrt = TMath::Max(maxfrt,nfront);
2019  iw[iwpos] = nfront;
2020  laell = ((nfront+1)*nfront)/2;
2021  apos2 = posfac+laell-1;
2022  if (numstk != 0) lnass = lnass*(2*nfront-lnass+1)/2;
2023 
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) {
2027  info[1] = -4;
2028  info[2] = la+TMath::Max(posfac+lnass,apos2-ltopst+2)-astk;
2029  goto finish;
2030  }
2031  }
2032 
2033  if (apos2 > azero) {
2034  apos = azero+1;
2035  lapos2 = TMath::Min(apos2,astk-1);
2036  if (lapos2 >= apos) {
2037  for (k= apos; k< lapos2+1; k++)
2038  a[k] = zero;
2039  }
2040  azero = apos2;
2041  }
2042 
2043  if (numstk != 0) {
2044  for (iell = 1; iell < numstk+1; iell++) {
2045  j1 = istk+1;
2046  j2 = istk+iw[istk];
2047  for (jj = j1; jj < j2+1; jj++) {
2048  irow = iw[jj];
2049  irow = iw2[irow];
2050  apos = posfac+IDiag(nfront,irow);
2051  for (jjj = jj; jjj < j2+1; jjj++) {
2052  j = iw[jjj];
2053  apos2 = apos+iw2[j]-irow;
2054  a[apos2] = a[apos2]+a[astk];
2055  a[astk] = zero;
2056  astk = astk+1;
2057  }
2058  }
2059  istk = j2+1;
2060  }
2061  }
2062 
2063  for (iorg = 1; iorg < numorg+1; iorg++) {
2064  j = -iw[iinput];
2065  irow = iw2[j];
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];
2070  ainput = ainput+1;
2071  iinput = iinput+1;
2072  if (iinput > liw) break;
2073  j = iw[iinput];
2074  if (j < 0) break;
2075  }
2076  }
2077  numass = numass+numorg;
2078  j1 = iwpos+2;
2079  j2 = iwpos+nfront+1;
2080  for (k = j1; k < j2+1; k++) {
2081  j = iw[k];
2082  iw2[j] = 0;
2083  }
2084 
2085  lnpiv = -1;
2086  npiv = 0;
2087  for (kdummy = 1; kdummy < nass+1; kdummy++) {
2088  if (npiv == nass) break;
2089  if (npiv == lnpiv) break;
2090  lnpiv = npiv;
2091  npivp1 = npiv+1;
2092  jpiv = 1;
2093  for (ipiv = npivp1; ipiv < nass+1; ipiv++) {
2094  jpiv = jpiv-1;
2095  if (jpiv == 1) continue;
2096  apos = posfac+IDiag(nfront-npiv,ipiv-npiv);
2097 
2098  if (uu <= zero) {
2099  if (TMath::Abs(a[apos]) <= cntl[3]) {
2100  info[1] = -5;
2101  info[2] = ntotpv+1;
2102  goto finish;
2103  }
2104  if (ntotpv <= 0) {
2105  if (a[apos] > zero) isnpiv = 1;
2106  if (a[apos] < zero) isnpiv = -1;
2107  }
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;
2111  info[1] = 2;
2112  i = ntotpv+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",
2115  info[1],i);
2116  isnpiv = -isnpiv;
2117  }
2118  if ((a[apos] > zero && isnpiv == 1) || (a[apos] < zero && isnpiv == -1) || (uu == zero)) goto hack;
2119  info[1] = -6;
2120  info[2] = ntotpv+1;
2121  goto finish;
2122  }
2123 
2124  amax = zero;
2125  tmax = amax;
2126  j1 = apos+1;
2127  j2 = apos+nass-ipiv;
2128  if (j2 >= j1) {
2129  for (jj = j1; jj < j2+1; jj++) {
2130  if (TMath::Abs(a[jj]) <= amax) continue;
2131  jmax = ipiv+jj-j1+1;
2132  amax = TMath::Abs(a[jj]);
2133  }
2134  }
2135  j1 = j2+1;
2136  j2 = apos+nfront-ipiv;
2137  if (j2 >= j1) {
2138  for (jj = j1; jj < j2+1; jj++)
2139  tmax = TMath::Max(TMath::Abs(a[jj]),tmax);
2140  }
2141  rmax = TMath::Max(tmax,amax);
2142  apos1 = apos;
2143  kk = nfront-ipiv;
2144  lt = ipiv-(npiv+1);
2145  if (lt != 0) {
2146  for (k = 1; k < lt+1; k++) {
2147  kk = kk+1;
2148  apos1 = apos1-kk;
2149  rmax = TMath::Max(rmax,TMath::Abs(a[apos1]));
2150  }
2151  }
2152  if (TMath::Abs(a[apos]) <= TMath::Max(cntl[3],uu*rmax)) {
2153  if (TMath::Abs(amax) <= cntl[3]) continue;
2154  apos2 = posfac+IDiag(nfront-npiv,jmax-npiv);
2155  detpiv = a[apos]*a[apos2]-amax*amax;
2156  thresh = TMath::Abs(detpiv);
2157  thresh = thresh/(uu*TMath::Max(TMath::Abs(a[apos])+amax,TMath::Abs(a[apos2])+amax));
2158  if (thresh <= rmax) continue;
2159  rmax = zero;
2160  j1 = apos2+1;
2161  j2 = apos2+nfront-jmax;
2162  if (j2 >= j1) {
2163  for (jj = j1; jj < j2+1; jj++)
2164  rmax = TMath::Max(rmax,TMath::Abs(a[jj]));
2165  }
2166  kk = nfront-jmax+1;
2167  apos3 = apos2;
2168  jmxmip = jmax-ipiv-1;
2169  if (jmxmip != 0) {
2170  for (k = 1; k < jmxmip+1; k++) {
2171  apos2 = apos2-kk;
2172  kk = kk+1;
2173  rmax = TMath::Max(rmax,TMath::Abs(a[apos2]));
2174  }
2175  }
2176  ipmnp = ipiv-npiv-1;
2177  if (ipmnp != 0) {
2178  apos2 = apos2-kk;
2179  kk = kk+1;
2180  for (k = 1; k < ipmnp+1; k++) {
2181  apos2 = apos2-kk;
2182  kk = kk+1;
2183  rmax = TMath::Max(rmax,TMath::Abs(a[apos2]));
2184  }
2185  }
2186  if (thresh <= rmax) continue;
2187  pivsiz = 2;
2188  } else {
2189  pivsiz = 1;
2190  }
2191 
2192  irow = ipiv-npiv;
2193  for (krow = 1; krow < pivsiz+1; krow++) {
2194  if (irow != 1) {
2195  j1 = posfac+irow;
2196  j2 = posfac+nfront-(npiv+1);
2197  if (j2 >= j1) {
2198  apos2 = apos+1;
2199  for (jj = j1; jj < j2+1; jj++) {
2200  swop = a[apos2];
2201  a[apos2] = a[jj];
2202  a[jj] = swop;
2203  apos2 = apos2+1;
2204  }
2205  }
2206  j1 = posfac+1;
2207  j2 = posfac+irow-2;
2208  apos2 = apos;
2209  kk = nfront-(irow+npiv);
2210  if (j2 >= j1) {
2211  for (jjj = j1; jjj < j2+1; jjj++) {
2212  jj = j2-jjj+j1;
2213  kk = kk+1;
2214  apos2 = apos2-kk;
2215  swop = a[apos2];
2216  a[apos2] = a[jj];
2217  a[jj] = swop;
2218  }
2219  }
2220  if (npiv != 0) {
2221  apos1 = posfac;
2222  kk = kk+1;
2223  apos2 = apos2-kk;
2224  for (jj = 1; jj < npiv+1; jj++) {
2225  kk = kk+1;
2226  apos1 = apos1-kk;
2227  apos2 = apos2-kk;
2228  swop = a[apos2];
2229  a[apos2] = a[apos1];
2230  a[apos1] = swop;
2231  }
2232  }
2233  swop = a[apos];
2234  a[apos] = a[posfac];
2235  a[posfac] = swop;
2236  ipos = iwpos+npiv+2;
2237  iexch = iwpos+irow+npiv+1;
2238  iswop = iw[ipos];
2239  iw[ipos] = iw[iexch];
2240  iw[iexch] = iswop;
2241  }
2242  if (pivsiz == 1) continue;
2243  if (krow != 2) {
2244  irow = jmax-(npiv+1);
2245  jpos = posfac;
2246  posfac = posfac+nfront-npiv;
2247  npiv = npiv+1;
2248  apos = apos3;
2249  } else {
2250  npiv = npiv-1;
2251  posfac = jpos;
2252  }
2253  }
2254 
2255  if (pivsiz != 2) {
2256 hack:
2257  a[posfac] = one/a[posfac];
2258  if (a[posfac] < zero) neig = neig+1;
2259  j1 = posfac+1;
2260  j2 = posfac+nfront-(npiv+1);
2261  if (j2 >= j1) {
2262  ibeg = j2+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];
2269  }
2270  ibeg = iend+1;
2271  a[jj] = amult;
2272  }
2273  }
2274  npiv = npiv+1;
2275  ntotpv = ntotpv+1;
2276  jpiv = 1;
2277  posfac = posfac+nfront-npiv+1;
2278  } else {
2279  ipos = iwpos+npiv+2;
2280  ntwo = ntwo+1;
2281  iw[ipos] = -iw[ipos];
2282  pospv1 = posfac;
2283  pospv2 = posfac+nfront-npiv;
2284  swop = a[pospv2];
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;
2290  j1 = pospv1+2;
2291  j2 = pospv1+nfront-(npiv+1);
2292  if (j2 >= j1) {
2293  jj1 = pospv2;
2294  ibeg = pospv2+nfront-(npiv+1);
2295  for (jj = j1; jj < j2+1; jj++) {
2296  jj1 = jj1+1;
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++) {
2301  k1 = jj+irow-ibeg;
2302  k2 = jj1+irow-ibeg;
2303  a[irow] = a[irow]+amult1*a[k1]+amult2*a[k2];
2304  }
2305  ibeg = iend+1;
2306  a[jj] = amult1;
2307  a[jj1] = amult2;
2308  }
2309  }
2310  npiv = npiv+2;
2311  ntotpv = ntotpv+2;
2312  jpiv = 2;
2313  posfac = pospv2+nfront-npiv+1;
2314  }
2315  }
2316  }
2317 
2318  if (npiv != 0) nblk = nblk+1;
2319  ioldps = iwpos;
2320  iwpos = iwpos+nfront+2;
2321  if (npiv != 0) {
2322  if (npiv <= 1) {
2323  iw[ioldps] = -iw[ioldps];
2324  for (k = 1; k < nfront+1; k++) {
2325  j1 = ioldps+k;
2326  iw[j1] = iw[j1+1];
2327  }
2328  iwpos = iwpos-1;
2329  } else {
2330  iw[ioldps+1] = npiv;
2331  }
2332  }
2333  liell = nfront-npiv;
2334 
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;
2339  iw[istk] = liell;
2340  j1 = istk;
2341  kk = iwpos-liell-1;
2342  for (k = 1; k < liell+1; k++) {
2343  j1 = j1+1;
2344  kk = kk+1;
2345  iw[j1] = iw[kk];
2346  }
2347  laell = ((liell+1)*liell)/2;
2348  kk = posfac+laell;
2349  if (kk == astk) {
2350  astk = astk-laell;
2351  } else {
2352  kmax = kk-1;
2353  for (k = 1; k < laell+1; k++) {
2354  kk = kk-1;
2355  astk = astk-1;
2356  a[astk] = a[kk];
2357  }
2358  kmax = TMath::Min(kmax,astk-1);
2359  for ( k = kk; k < kmax+1; k++)
2360  a[k] = zero;
2361  }
2362  azero = TMath::Min(azero,astk-1);
2363  }
2364  if (npiv == 0) iwpos = ioldps;
2365  }
2366 
2367  iw[1] = nblk;
2368  if (ntwo > 0) iw[1] = -nblk;
2369  nrlbdu = posfac-1;
2370  nirbdu = iwpos-1;
2371 
2372  if (ntotpv != n) {
2373  info[1] = 3;
2374  info[2] = ntotpv;
2375  }
2376 
2377 finish:
2378  info[9] = nrlbdu;
2379  info[10] = nirbdu;
2380  info[12] = ncmpbr;
2381  info[13] = ncmpbi;
2382  info[14] = ntwo;
2383  info[15] = neig;
2384 }
2385 
2386 ////////////////////////////////////////////////////////////////////////////////
2387 /// Help routine for factorization
2388 
2390  const Int_t ireal,Int_t &ncmpbr,Int_t &ncmpbi)
2391 {
2392  Int_t ipos,jj,jjj;
2393 
2394  ipos = itop-1;
2395  if (j2 != ipos) {
2396  if (ireal != 2) {
2397  ncmpbr = ncmpbr+1;
2398  if (j1 <= j2) {
2399  for (jjj = j1; jjj < j2+1; jjj++) {
2400  jj = j2-jjj+j1;
2401  a[ipos] = a[jj];
2402  ipos = ipos-1;
2403  }
2404  }
2405  } else {
2406  ncmpbi = ncmpbi+1;
2407  if (j1 <= j2) {
2408  for (jjj = j1; jjj < j2+1; jjj++) {
2409  jj = j2-jjj+j1;
2410  iw[ipos] = iw[jj];
2411  ipos = ipos-1;
2412  }
2413  }
2414  }
2415  j2 = itop-1;
2416  j1 = ipos+1;
2417  }
2418 }
2419 
2420 ////////////////////////////////////////////////////////////////////////////////
2421 /// Help routine for solving
2422 
2424  Double_t *rhs,Int_t *iw2,const Int_t nblk,Int_t &latop,
2425  Int_t *icntl)
2426 {
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;
2428  Double_t w1,w2;
2429 
2430  const Int_t ifrlvl = 5;
2431 
2432  apos = 1;
2433  ipos = 1;
2434  j2 = 0;
2435  iblk = 0;
2436  npiv = 0;
2437  for (irow = 1; irow < n+1; irow++) {
2438  if (npiv <= 0) {
2439  iblk = iblk+1;
2440  if (iblk > nblk) break;
2441  ipos = j2+1;
2442  iw2[iblk] = ipos;
2443  liell = -iw[ipos];
2444  npiv = 1;
2445  if (liell <= 0) {
2446  liell = -liell;
2447  ipos = ipos+1;
2448  npiv = iw[ipos];
2449  }
2450  j1 = ipos+1;
2451  j2 = ipos+liell;
2452  ilvl = TMath::Min(npiv,10);
2453  if (liell < icntl[ifrlvl+ilvl]) goto hack;
2454  ifr = 0;
2455  for (jj = j1; jj < j2+1; jj++) {
2456  j = TMath::Abs(iw[jj]+0);
2457  ifr = ifr+1;
2458  w[ifr] = rhs[j];
2459  }
2460  jpiv = 1;
2461  j3 = j1;
2462 
2463  for (ipiv = 1; ipiv < npiv+1; ipiv++) {
2464  jpiv = jpiv-1;
2465  if (jpiv == 1) continue;
2466 
2467  if (iw[j3] >= 0) {
2468  jpiv = 1;
2469  j3 = j3+1;
2470  apos = apos+1;
2471  ist = ipiv+1;
2472  if (liell< ist) continue;
2473  w1 = w[ipiv];
2474  k = apos;
2475  for (j = ist; j < liell+1; j++) {
2476  w[j] = w[j]+a[k]*w1;
2477  k = k+1;
2478  }
2479  apos = apos+liell-ist+1;
2480  } else {
2481  jpiv = 2;
2482  j3 = j3+2;
2483  apos = apos+2;
2484  ist = ipiv+2;
2485  if (liell >= ist) {
2486  w1 = w[ipiv];
2487  w2 = w[ipiv+1];
2488  k1 = apos;
2489  k2 = apos+liell-ipiv;
2490  for (j = ist; j < liell+1; j++) {
2491  w[j] = w[j]+w1*a[k1]+w2*a[k2];
2492  k1 = k1+1;
2493  k2 = k2+1;
2494  }
2495  }
2496  apos = apos+2*(liell-ist+1)+1;
2497  }
2498  }
2499 
2500  ifr = 0;
2501  for (jj = j1; jj < j2+1; jj++) {
2502  j = TMath::Abs(iw[jj]+0);
2503  ifr = ifr+1;
2504  rhs[j] = w[ifr];
2505  }
2506  npiv = 0;
2507  } else {
2508 hack:
2509  if (iw[j1] >= 0) {
2510  npiv = npiv-1;
2511  apos = apos+1;
2512  j1 = j1+1;
2513  if (j1 <= j2) {
2514  irhs = iw[j1-1];
2515  w1 = rhs[irhs];
2516  k = apos;
2517  for (j = j1; j < j2+1; j++) {
2518  irhs = TMath::Abs(iw[j]+0);
2519  rhs[irhs] = rhs[irhs]+a[k]*w1;
2520  k = k+1;
2521  }
2522  }
2523  apos = apos+j2-j1+1;
2524  } else {
2525  npiv = npiv-2;
2526  j1 = j1+2;
2527  apos = apos+2;
2528  if (j1 <= j2) {
2529  irhs = -iw[j1-2];
2530  w1 = rhs[irhs];
2531  irhs = iw[j1-1];
2532  w2 = rhs[irhs];
2533  k1 = apos;
2534  k3 = apos+j2-j1+2;
2535  for (j = j1; j < j2+1; j++) {
2536  irhs = TMath::Abs(iw[j]+0);
2537  rhs[irhs] = rhs[irhs]+w1*a[k1]+w2*a[k3];
2538  k1 = k1+1;
2539  k3 = k3+1;
2540  }
2541  }
2542  apos = apos+2*(j2-j1+1)+1;
2543  }
2544  }
2545  }
2546 
2547  latop = apos-1;
2548 }
2549 
2550 ////////////////////////////////////////////////////////////////////////////////
2551 /// Help routine for solving
2552 
2554  Double_t *rhs,Int_t *iw2,const Int_t nblk,
2555  const Int_t latop,Int_t *icntl)
2556 {
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;
2559  Double_t w1,w2;
2560 
2561  const Int_t ifrlvl = 5;
2562 
2563  apos = latop+1;
2564  npiv = 0;
2565  iblk = nblk+1;
2566  for (loop = 1; loop < n+1; loop++) {
2567  if (npiv <= 0) {
2568  iblk = iblk-1;
2569  if (iblk < 1) break;
2570  ipos = iw2[iblk];
2571  liell = -iw[ipos];
2572  npiv = 1;
2573  if (liell <= 0) {
2574  liell = -liell;
2575  ipos = ipos+1;
2576  npiv = iw[ipos];
2577  }
2578  jpos = ipos+npiv;
2579  j2 = ipos+liell;
2580  ilvl = TMath::Min(10,npiv)+10;
2581  if (liell < icntl[ifrlvl+ilvl]) goto hack;
2582  j1 = ipos+1;
2583  ifr = 0;
2584  for (jj = j1; jj < j2+1; jj++) {
2585  j = TMath::Abs(iw[jj]+0);
2586  ifr = ifr+1;
2587  w[ifr] = rhs[j];
2588  }
2589  jpiv = 1;
2590  for (iipiv = 1; iipiv < npiv+1; iipiv++) {
2591  jpiv = jpiv-1;
2592  if (jpiv == 1) continue;
2593  ipiv = npiv-iipiv+1;
2594  if (ipiv == 1 || iw[jpos-1] >= 0) {
2595  jpiv = 1;
2596  apos = apos-(liell+1-ipiv);
2597  ist = ipiv+1;
2598  w1 = w[ipiv]*a[apos];
2599  if (liell >= ist) {
2600  jj1 = apos+1;
2601  for (j = ist; j < liell+1; j++) {
2602  w1 = w1+a[jj1]*w[j];
2603  jj1 = jj1+1;
2604  }
2605  }
2606  w[ipiv] = w1;
2607  jpos = jpos-1;
2608  } else {
2609  jpiv = 2;
2610  apos2 = apos-(liell+1-ipiv);
2611  apos = apos2-(liell+2-ipiv);
2612  ist = ipiv+1;
2613  w1 = w[ipiv-1]*a[apos]+w[ipiv]*a[apos+1];
2614  w2 = w[ipiv-1]*a[apos+1]+w[ipiv]*a[apos2];
2615  if (liell >= ist) {
2616  jj1 = apos+2;
2617  jj2 = apos2+1;
2618  for (j = ist; j < liell+1; j++) {
2619  w1 = w1+w[j]*a[jj1];
2620  w2 = w2+w[j]*a[jj2];
2621  jj1 = jj1+1;
2622  jj2 = jj2+1;
2623  }
2624  }
2625  w[ipiv-1] = w1;
2626  w[ipiv] = w2;
2627  jpos = jpos-2;
2628  }
2629  }
2630  ifr = 0;
2631  for (jj = j1; jj < j2+1; jj++) {
2632  j = TMath::Abs(iw[jj]+0);
2633  ifr = ifr+1;
2634  rhs[j] = w[ifr];
2635  }
2636  npiv = 0;
2637  } else {
2638 hack:
2639  if (npiv == 1 || iw[jpos-1] >= 0) {
2640  npiv = npiv-1;
2641  apos = apos-(j2-jpos+1);
2642  iirhs = iw[jpos];
2643  w1 = rhs[iirhs]*a[apos];
2644  j1 = jpos+1;
2645  if (j1 <= j2) {
2646  k = apos+1;
2647  for (j = j1; j < j2+1; j++) {
2648  irhs = TMath::Abs(iw[j]+0);
2649  w1 = w1+a[k]*rhs[irhs];
2650  k = k+1;
2651  }
2652  }
2653  rhs[iirhs] = w1;
2654  jpos = jpos-1;
2655  } else {
2656  npiv = npiv-2;
2657  apos2 = apos-(j2-jpos+1);
2658  apos = apos2-(j2-jpos+2);
2659  i1rhs = -iw[jpos-1];
2660  i2rhs = iw[jpos];
2661  w1 = rhs[i1rhs]*a[apos]+rhs[i2rhs]*a[apos+1];
2662  w2 = rhs[i1rhs]*a[apos+1]+rhs[i2rhs]*a[apos2];
2663  j1 = jpos+1;
2664  if (j1 <= j2) {
2665  jj1 = apos+2;
2666  jj2 = apos2+1;
2667  for (j = j1; j < j2+1; j++) {
2668  irhs = TMath::Abs(iw[j]+0);
2669  w1 = w1+rhs[irhs]*a[jj1];
2670  w2 = w2+rhs[irhs]*a[jj2];
2671  jj1 = jj1+1;
2672  jj2 = jj2+1;
2673  }
2674  }
2675  rhs[i1rhs] = w1;
2676  rhs[i2rhs] = w2;
2677  jpos = jpos-2;
2678  }
2679  }
2680  }
2681 }
2682 
2683 ////////////////////////////////////////////////////////////////////////////////
2684 /// Print class members
2685 
2687 {
2688  TDecompBase::Print(opt);
2689 
2690  printf("fPrecision = %.3f\n",fPrecision);
2691  printf("fIPessimism = %.3f\n",fIPessimism);
2692  printf("fRPessimism = %.3f\n",fRPessimism);
2693 
2694  TMatrixDSparse fact(0,fNrows-1,0,fNrows-1,fNnonZeros,
2696  fact.Print("fFact");
2697 }
2698 
2699 ////////////////////////////////////////////////////////////////////////////////
2700 /// Assignment operator
2701 
2703 {
2704  if (this != &source) {
2705  TDecompBase::operator=(source);
2706  memcpy(fIcntl,source.fIcntl,31*sizeof(Int_t));
2707  memcpy(fCntl,source.fCntl,6*sizeof(Double_t));
2708  memcpy(fInfo,source.fInfo,21*sizeof(Int_t));
2709  fVerbose = source.fVerbose;
2710  fPrecision = source.fPrecision;
2711  fIkeep = source.fIkeep;
2712  fIw = source.fIw;
2713  fIw1 = source.fIw1;
2714  fIw2 = source.fIw2;
2715  fNsteps = source.fNsteps;
2716  fMaxfrt = source.fMaxfrt;
2717  fW = source.fW;
2718  fIPessimism = source.fIPessimism;
2719  fRPessimism = source.fRPessimism;
2720  if (fA.IsValid())
2721  fA.Use(*const_cast<TMatrixDSparse *>(&(source.fA)));
2722  fNrows = source.fNrows;
2723  fNnonZeros = source.fNnonZeros;
2724  fFact = source.fFact;
2725  fRowFact = source.fRowFact;
2726  fColFact = source.fColFact;
2727  }
2728  return *this;
2729 }
virtual const Element * GetMatrixArray() const
void SetTreatAsZero(Double_t tol)
Double_t fCntl[6]
Definition: TDecompSparse.h:50
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:854
TMatrixDSparse fA
Definition: TDecompSparse.h:69
TArrayI fColFact
Definition: TDecompSparse.h:75
Double_t fRPessimism
Definition: TDecompSparse.h:66
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:120
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.
Double_t fIPessimism
Definition: TDecompSparse.h:65
Int_t fRowLwb
Definition: TDecompBase.h:40
const char Option_t
Definition: RtypesCore.h:62
void Reset()
Definition: TArrayD.h:47
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
Definition: TArrayD.h:43
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
Definition: TDecompSparse.h:41
Int_t GetLwb() const
Definition: TVectorT.h:73
TArrayI fIkeep
Definition: TDecompSparse.h:57
void Print(Option_t *opt="") const
Print class members.
TDecompSparse()
Default constructor.
Decomposition Base class.
Definition: TDecompBase.h:33
#define R__ASSERT(e)
Definition: TError.h:96
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
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.
Int_t GetNrows() const
Definition: TVectorT.h:75
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition: TVectorT.cxx:601
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
Array of integers (32 bits per element).
Definition: TArrayI.h:27
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual Bool_t Decompose()
Decomposition engine .
Int_t GetColLwb() const
Definition: TMatrixTBase.h:123
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
Definition: TArrayI.h:43
Int_t fColLwb
Definition: TDecompBase.h:41
XFontStruct * id
Definition: TGX11.cxx:108
void Set(Int_t n)
Set size of this array to n ints.
Definition: TArrayI.cxx:105
TArrayI fRowFact
Definition: TDecompSparse.h:74
virtual void ls(Option_t *option="") const
The ls function lists the contents of a class on stdout.
Definition: TObject.cxx:492
TVectorT< Element > & Shift(Int_t row_shift)
Definition: TVectorT.h:87
Element * GetMatrixArray()
Definition: TVectorT.h:78
static Int_t NonZerosUpperTriang(const TMatrixDSparse &a)
Static function, returning the number of non-zero entries in the upper triangular matrix ...
Int_t fInfo[21]
Definition: TDecompSparse.h:51
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.
Int_t GetSize() const
Definition: TArray.h:47
auto * a
Definition: textangle.C:12
Bool_t IsValid() const
Definition: TVectorT.h:83
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
void Print(Option_t *opt="") const
Print class members.
Double_t GetThresholdPivoting()
Sparse Symmetric Decomposition class.
Definition: TDecompSparse.h:43
void ResetStatus()
Definition: TDecompBase.h:43
void SetThresholdPivoting(Double_t piv)
static constexpr double ms
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.
Int_t fIcntl[31]
Definition: TDecompSparse.h:49
virtual const Int_t * GetColIndexArray() const
Int_t GetNrows() const
Definition: TMatrixTBase.h:122
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t fPrecision
Definition: TDecompSparse.h:53
const Double_t kInitPrecision
Definition: TDecompSparse.h:29
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
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
Definition: TDecompSparse.h:36
Int_t MinRealWorkspace()
Array of doubles (64 bits per element).
Definition: TArrayD.h:27
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
Definition: TDecompSparse.h:27
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.
auto * l
Definition: textangle.C:4
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
Int_t ErrorFlag()
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
Definition: TRolke.cxx:630
const Double_t kInitThresholdPivoting
Definition: TDecompSparse.h:28
int * iq
Definition: THbookFile.cxx:86
void SetVerbose(Int_t v)
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.
const Bool_t kTRUE
Definition: RtypesCore.h:87
Int_t Nint(T x)
Definition: TMath.h:606
void Set(Int_t n)
Set size of this array to n doubles.
Definition: TArrayD.cxx:106
static char * skip(char **buf, const char *delimiters)
Definition: civetweb.c:2039
static constexpr double ns
const Int_t n
Definition: legend1.C:16
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
Bool_t IsValid() const
Definition: TMatrixTBase.h:145
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.