Logo ROOT   6.16/01
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
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
67TDecompSparse::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
470void 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
581void 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
836void 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
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
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
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
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
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) {
2256hack:
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
2377finish:
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 {
2508hack:
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 {
2638hack:
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
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}
#define b(i)
Definition: RSha256.hxx:100
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:363
const Double_t kThresholdPivotingFactor
Definition: TDecompSparse.h:41
const Double_t kThresholdPivotingMax
Definition: TDecompSparse.h:36
const Double_t kInitThresholdPivoting
Definition: TDecompSparse.h:28
const Double_t kInitTreatAsZero
Definition: TDecompSparse.h:27
const Double_t kInitPrecision
Definition: TDecompSparse.h:29
#define R__ASSERT(e)
Definition: TError.h:96
int * iq
Definition: THbookFile.cxx:86
Array of doubles (64 bits per element).
Definition: TArrayD.h:27
void Set(Int_t n)
Set size of this array to n doubles.
Definition: TArrayD.cxx:106
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)
Set size of this array to n ints.
Definition: TArrayI.cxx:105
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
void Print(Option_t *opt="") const
Print class members.
Int_t fRowLwb
Definition: TDecompBase.h:40
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
Int_t fColLwb
Definition: TDecompBase.h:41
Sparse Symmetric Decomposition class.
Definition: TDecompSparse.h:44
Double_t fRPessimism
Definition: TDecompSparse.h:66
virtual Bool_t Decompose()
Decomposition engine .
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 .
TArrayI fIkeep
Definition: TDecompSparse.h:57
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)
void Print(Option_t *opt="") const
Print class members.
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 .
TArrayI fRowFact
Definition: TDecompSparse.h:74
Double_t fPrecision
Definition: TDecompSparse.h:53
void SetThresholdPivoting(Double_t piv)
Int_t fIcntl[31]
Definition: TDecompSparse.h:49
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.
TArrayI fColFact
Definition: TDecompSparse.h:75
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
Definition: TDecompSparse.h:65
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
Definition: TDecompSparse.h:69
Int_t fInfo[21]
Definition: TDecompSparse.h:51
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.
Int_t ErrorFlag()
Double_t fCntl[6]
Definition: TDecompSparse.h:50
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
void Print(Option_t *name="") const
Print the matrix as a table of elements.
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:122
Int_t GetColLwb() const
Definition: TMatrixTBase.h:125
Bool_t IsValid() const
Definition: TMatrixTBase.h:147
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:172
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
virtual void ls(Option_t *option="") const
The ls function lists the contents of a class on stdout.
Definition: TObject.cxx:492
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:854
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition: TVectorT.cxx:601
const Int_t n
Definition: legend1.C:16
static constexpr double ms
static constexpr double ns
Int_t Nint(T x)
Round to nearest integer. Rounds half integers to the nearest even integer.
Definition: TMath.h:701
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12