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