Logo ROOT   6.16/01
Reference Guide
TCernLib.cxx
Go to the documentation of this file.
1// @(#)root/table:$Id$
2// Author: Valery Fine(fine@bnl.gov) 25/09/99
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////////////////////////////////////////////////////////////////////////////////
13// The set of methods to work with the plain matrix / vector
14// "derived" from https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f110/top.html
15// "derived" from https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html
16//
17// Revision 1.7 2006/05/21 18:05:26 brun
18// Fix more coding conventions violations
19//
20// Revision 1.6 2006/05/20 14:06:09 brun
21// Fix a VERY long list of coding conventions violations
22//
23// Revision 1.5 2003/09/30 09:52:49 brun
24// Add references to the original CERNLIB packages
25//
26// Revision 1.4 2003/05/28 15:17:03 brun
27// From Valeri Fine. A new version of the table package.
28// It fixes a couple of memory leaks:
29// class TTableDescriptorm
30// class TVolumePosition
31// and provides some clean up
32// for the TCL class interface.
33//
34// Revision 1.3 2003/04/03 17:39:39 fine
35// Make merge with ROOT 3.05.03 and add TR package
36//
37// Revision 1.2 2003/02/04 23:35:20 fine
38// Clean up
39//
40// Revision 1.1 2002/04/15 20:23:39 fine
41// New naming schema for RootKErnel classes and a set of classes to back geometry OO
42//
43// Revision 1.2 2001/05/29 19:08:08 brun
44// New version of some STAR classes from Valery.
45//
46// Revision 1.2 2001/05/27 02:38:14 fine
47// New method trsedu to solev Ax=B from Victor
48//
49// Revision 1.1.1.1 2000/11/27 22:57:14 fisyak
50//
51//
52// Revision 1.1.1.1 2000/05/16 17:00:48 rdm
53// Initial import of ROOT into CVS
54//
55////////////////////////////////////////////////////////////////////////////////
56
57#include <assert.h>
58#include "TCernLib.h"
59#include "TMath.h"
60#include "TArrayD.h"
61#include "TError.h"
62
64
65#define TCL_MXMAD(n_,a,b,c,i,j,k) \
66 /* Local variables */ \
67 int l, m, n, ia, ic, ib, ja, jb, iia, iib, ioa, iob; \
68 \
69 /* Parameter adjustments */ \
70 --a; --b; --c; \
71 /* Function Body */ \
72/* MXMAD MXMAD1 MXMAD2 MXMAD3 MXMPY MXMPY1 MXMPY2 MXMPY3 MXMUB MXMUB1 MXMUB2 MXMUB3 */ \
73/* const int iandj1[] = {21, 22, 23, 24, 11, 12, 13, 14, 31, 32, 33, 34 }; */ \
74 const int iandj1[] = {2, 2 , 2 , 2 , 1 , 1 , 1 , 1 , 3 , 3 , 3 , 3 }; \
75 const int iandj2[] = { 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4 }; \
76 int n1 = iandj1[n_]; \
77 int n2 = iandj2[n_]; \
78 if (i == 0 || k == 0) return 0; \
79 \
80 switch (n2) { \
81 case 1: iia = 1; ioa = j; iib = k; iob = 1; break; \
82 case 2: iia = 1; ioa = j; iib = 1; iob = j; break; \
83 case 3: iia = i; ioa = 1; iib = k; iob = 1; break; \
84 case 4: iia = i; ioa = 1; iib = 1; iob = j; break; \
85 default: iia = ioa = iib = iob = 0; assert(iob); \
86 }; \
87 \
88 ia = 1; ic = 1; \
89 for (l = 1; l <= i; ++l) { \
90 ib = 1; \
91 for (m = 1; m <= k; ++m,++ic) { \
92 switch (n1) { \
93 case 1: c[ic] = 0.; break; \
94 case 3: c[ic] = -c[ic]; break; \
95 }; \
96 if (j == 0) continue; \
97 ja = ia; jb = ib; \
98 double cic = c[ic]; \
99 for (n = 1; n <= j; ++n, ja+=iia, jb+=iib) \
100 cic += a[ja] * b[jb]; \
101 c[ic] = cic; \
102 ib += iob; \
103 } \
104 ia += ioa; \
105 }
106
107////////////////////////////////////////////////////////////////////////////////
108
109float *TCL::mxmad_0_(int n_, const float *a, const float *b, float *c, int i, int j, int k)
110{
111 TCL_MXMAD(n_,a,b,c,i,j,k)
112 return c;
113} /* mxmad_ */
114
115////////////////////////////////////////////////////////////////////////////////
116
117double *TCL::mxmad_0_(int n_, const double *a, const double *b, double *c, int i, int j, int k)
118{
119 TCL_MXMAD(n_,a,b,c,i,j,k)
120 return c;
121} /* mxmad_ */
122
123#undef TCL_MXMAD
124
125//___________________________________________________________________________
126//
127// Matrix Multiplication
128//___________________________________________________________________________
129
130#define TCL_MXMLRT( n__, a, b, c, ni,nj) \
131 if (ni <= 0 || nj <= 0) return 0; \
132 double x; \
133 int ia, ib, ic, ja, kc, ii, jj, kj, ki, ia1, ib1, ic1, ja1; \
134 int ipa = 1; int jpa = nj; \
135 if (n__ == 1) { ipa = ni; jpa = 1; } \
136 \
137 --a; --b; --c; \
138 \
139 ic1 = 1; ia1 = 1; \
140 for (ii = 1; ii <= ni; ++ii, ic1+=ni, ia1+=jpa) { \
141 ic = ic1; \
142 for (kc = 1; kc <= ni; ++kc,ic++) c[ic] = 0.; \
143 ib1 = 1; ja1 = 1; \
144 for (jj = 1; jj <= nj; ++jj,++ib1,ja1 += ipa) { \
145 ib = ib1; ia = ia1; \
146 x = 0.; \
147 for (kj = 1;kj <= nj;++kj,ia+=ipa,ib += nj) \
148 x += a[ia] * b[ib]; \
149 ja = ja1; ic = ic1; \
150 for (ki = 1; ki <= ni; ++ki,++ic,ja += jpa) \
151 c[ic] += x * a[ja]; \
152 } \
153 }
154
155////////////////////////////////////////////////////////////////////////////////
156/// Matrix Multiplication
157///
158/// CERN PROGLIB# F110 MXMLRT .VERSION KERNFOR 2.00 720707
159/// ORIG. 01/01/64 RKB
160///
161/// See original documentation of CERNLIB package
162/// [F110](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f110/top.html)
163///
164/// -- ENTRY MXMLRT
165/// -- C = A(I,J) X B(J,J) X A*(J,I)
166/// -- A* TANDS FOR A-TRANSPOSED
167/// mxmlrt (A,B,C,NI,NJ) IS EQUIVALENT TO
168/// CALL MXMPY (A,B,X,NI,NJ,NJ)
169/// CALL MXMPY1 (X,A,C,NI,NJ,NI)
170///
171/// OR CALL MXMPY1 (B,A,Y,NJ,NJ,NI)
172/// CALL MXMPY (A,Y,C,NI,NJ,NI)
173///
174///
175/// -- C = A*(I,J) X B(J,J) X A(J,I)
176///
177/// CALL MXMLTR (A,B,C,NI,NJ) IS EQUIVALENT TO
178/// CALL MXMPY2 (A,B,X,NI,NJ,NJ)
179/// CALL MXMPY (X,A,C,NI,NJ,NI)
180///
181/// OR CALL MXMPY (B,A,Y,NJ,NJ,NI)
182/// CALL MXMPY2 (A,Y,C,NI,NJ,NI)
183
184float *TCL::mxmlrt_0_(int n__, const float *a, const float *b, float *c, int ni,int nj)
185{
186 TCL_MXMLRT( n__, a, b, c, ni,nj)
187 return c;
188} /* mxmlrt_ */
189
190////////////////////////////////////////////////////////////////////////////////
191/// Matrix Multiplication (double precision)
192///
193/// CERN PROGLIB# F110 MXMLRT .VERSION KERNFOR 2.00 720707
194/// ORIG. 01/01/64 RKB
195///
196/// See original documentation of CERNLIB package
197/// [F110](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f110/top.html)
198
199double *TCL::mxmlrt_0_(int n__, const double *a, const double *b, double *c, int ni,int nj)
200{
201 TCL_MXMLRT( n__, a, b, c, ni,nj)
202 return c;
203
204} /* mxmlrt_ */
205
206#undef TCL_MXMLRT
207
208//___________________________________________________________________________
209//
210// Matrix Transposition
211//___________________________________________________________________________
212
213#define TCL_MXTRP(a, b, i, j) \
214 if (i == 0 || j == 0) return 0; \
215 --b; --a; \
216 int ib = 1; \
217 for (int k = 1; k <= j; ++k) \
218 { int ia = k; \
219 for (int l = 1; l <= i; ++l,ia += j,++ib) b[ib] = a[ia]; }
220
221////////////////////////////////////////////////////////////////////////////////
222/// Matrix Transposition
223///
224/// CERN PROGLIB# F110 MXTRP .VERSION KERNFOR 1.0 650809
225/// ORIG. 01/01/64 RKB
226///
227/// See original documentation of CERNLIB package
228/// [F110](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f110/top.html)
229
230float *TCL::mxtrp(const float *a, float *b, int i, int j)
231{
232 TCL_MXTRP(a, b, i, j)
233 return b;
234} /* mxtrp */
235
236////////////////////////////////////////////////////////////////////////////////
237/// Matrix Transposition (double precision)
238///
239/// CERN PROGLIB# F110 MXTRP .VERSION KERNFOR 1.0 650809
240/// ORIG. 01/01/64 RKB
241///
242/// See original documentation of CERNLIB package
243/// [F110](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f110/top.html)
244
245double *TCL::mxtrp(const double *a, double *b, int i, int j)
246{
247 TCL_MXTRP(a, b, i, j)
248 return b;
249
250} /* mxtrp */
251#undef TCL_MXTRP
252
253//___________________________________________________________________________
254//___________________________________________________________________________
255//
256// TRPACK
257//___________________________________________________________________________
258//___________________________________________________________________________
259
260#define TCL_TRAAT(a, s, m, n) \
261 /* Local variables */ \
262 int ipiv, i, j, ipivn, ia, is, iat; \
263 double sum; \
264 --s; --a; \
265 ia = 0; is = 0; \
266 for (i = 1; i <= m; ++i) { \
267 ipiv = ia; \
268 ipivn = ipiv + n; \
269 iat = 0; \
270 for (j = 1; j <= i; ++j) { \
271 ia = ipiv; \
272 sum = 0.; \
273 do { \
274 ++ia; ++iat; \
275 sum += a[ia] * a[iat]; \
276 } while (ia < ipivn); \
277 ++is; \
278 s[is] = sum; \
279 } \
280 } \
281 s++;
282
283////////////////////////////////////////////////////////////////////////////////
284/// Symmetric Multiplication of Rectangular Matrices
285///
286/// CERN PROGLIB# F112 TRAAT .VERSION KERNFOR 4.15 861204
287/// ORIG. 18/12/74 WH
288/// traat.F -- translated by f2c (version 19970219).
289///
290/// See original documentation of CERNLIB package
291/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
292
293float *TCL::traat(const float *a, float *s, int m, int n)
294{
295 TCL_TRAAT(a, s, m, n)
296 return s;
297} /* traat_ */
298
299////////////////////////////////////////////////////////////////////////////////
300/// Symmetric Multiplication of Rectangular Matrices
301///
302/// CERN PROGLIB# F112 TRAAT .VERSION KERNFOR 4.15 861204
303/// ORIG. 18/12/74 WH
304/// traat.F -- translated by f2c (version 19970219).
305///
306/// See original documentation of CERNLIB package
307/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
308
309double *TCL::traat(const double *a, double *s, int m, int n)
310{
311 TCL_TRAAT(a, s, m, n)
312 return s;
313} /* traat_ */
314
315#undef TCL_TRAAT
316
317#define TCL_TRAL(a, u, b, m, n) \
318 int indu, i, j, k, ia, ib, iu; \
319 double sum; \
320 --b; --u; --a; \
321 ib = 1; \
322 for (i = 1; i <= m; ++i) { \
323 indu = 0; \
324 for (j = 1; j <= n; ++j) { \
325 indu += j; \
326 ia = ib; \
327 iu = indu; \
328 sum = 0.; \
329 for (k = j; k <= n; ++k) {\
330 sum += a[ia] * u[iu]; \
331 ++ia; \
332 iu += k; \
333 } \
334 b[ib] = sum; \
335 ++ib; \
336 } \
337 } \
338 b++;
339
340////////////////////////////////////////////////////////////////////////////////
341/// Triangular - Rectangular Multiplication
342///
343/// CERN PROGLIB# F112 TRAL .VERSION KERNFOR 4.15 861204
344/// ORIG. 18/12/74 WH
345/// tral.F -- translated by f2c (version 19970219).
346///
347/// See original documentation of CERNLIB package
348/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
349
350float *TCL::tral(const float *a, const float *u, float *b, int m, int n)
351{
352 TCL_TRAL(a, u, b, m, n)
353 return b;
354} /* tral_ */
355
356////////////////////////////////////////////////////////////////////////////////
357/// Triangular - Rectangular Multiplication
358///
359/// tral.F -- translated by f2c (version 19970219).
360/// CERN PROGLIB# F112 TRAL .VERSION KERNFOR 4.15 861204
361/// ORIG. 18/12/74 WH
362///
363/// See original documentation of CERNLIB package
364/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
365
366double *TCL::tral(const double *a, const double *u, double *b, int m, int n)
367{
368 TCL_TRAL(a, u, b, m, n)
369 return b;
370} /* tral_ */
371
372#undef TCL_TRAL
373
374////////////////////////////////////////////////////////////////////////////////
375
376#define TCL_TRALT(a, u, b, m, n) \
377 int indu, j, k, ia, ib, iu; \
378 double sum; \
379 --b; --u; --a; \
380 ib = m * n; \
381 indu = (n * n + n) / 2; \
382 do { \
383 iu = indu; \
384 for (j = 1; j <= n; ++j) { \
385 ia = ib; \
386 sum = 0.; \
387 for (k = j; k <= n; ++k) {\
388 sum += a[ia] * u[iu]; \
389 --ia; --iu; \
390 } \
391 b[ib] = sum; \
392 --ib; \
393 } \
394 } while (ib > 0); \
395 ++b;
396
397////////////////////////////////////////////////////////////////////////////////
398/// Triangular - Rectangular Multiplication
399///
400/// CERN PROGLIB# F112 TRALT .VERSION KERNFOR 4.15 861204
401/// ORIG. 18/12/74 WH
402/// tralt.F -- translated by f2c (version 19970219).
403///
404/// See original documentation of CERNLIB package
405/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
406
407float *TCL::tralt(const float *a, const float *u, float *b, int m, int n)
408{
409 TCL_TRALT(a, u, b, m, n)
410 return b;
411} /* tralt_ */
412
413////////////////////////////////////////////////////////////////////////////////
414/// Triangular - Rectangular Multiplication
415///
416/// CERN PROGLIB# F112 TRALT .VERSION KERNFOR 4.15 861204
417/// ORIG. 18/12/74 WH
418/// tralt.F -- translated by f2c (version 19970219).
419///
420/// See original documentation of CERNLIB package
421/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
422
423double *TCL::tralt(const double *a, const double *u, double *b, int m, int n)
424{
425 TCL_TRALT(a, u, b, m, n)
426 return b;
427} /* tralt_ */
428
429#undef TCL_TRALT
430
431//____________________________________________________________
432
433#define TCL_TRAS(a, s, b, m, n) \
434 int inds, i__, j, k, ia, ib, is; \
435 double sum; \
436 --b; --s; --a; \
437 ib = 0; inds = 0; i__ = 0; \
438 do { \
439 inds += i__; \
440 ia = 0; \
441 ib = i__ + 1; \
442 for (j = 1; j <= m; ++j) { \
443 is = inds; \
444 sum = 0.; \
445 k = 0; \
446 do { \
447 if (k > i__) is += k; \
448 else ++is; \
449 ++ia; \
450 sum += a[ia] * s[is]; \
451 ++k; \
452 } while (k < n); \
453 b[ib] = sum; \
454 ib += n; \
455 } \
456 ++i__; \
457 } while (i__ < n); \
458 ++b;
459
460////////////////////////////////////////////////////////////////////////////////
461/// Symmetric - Rectangular Multiplication
462///
463/// CERN PROGLIB# F112 TRAS .VERSION KERNFOR 4.15 861204
464/// ORIG. 18/12/74 WH
465/// tras.F -- translated by f2c (version 19970219).
466///
467/// See original documentation of CERNLIB package
468/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
469
470float *TCL::tras(const float *a, const float *s, float *b, int m, int n)
471{
472 TCL_TRAS(a, s, b, m, n)
473 return b;
474} /* tras_ */
475
476////////////////////////////////////////////////////////////////////////////////
477/// Symmetric - Rectangular Multiplication
478///
479/// CERN PROGLIB# F112 TRAS .VERSION KERNFOR 4.15 861204
480/// ORIG. 18/12/74 WH
481/// tras.F -- translated by f2c (version 19970219).
482///
483/// See original documentation of CERNLIB package
484/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
485
486double *TCL::tras(const double *a, const double *s, double *b, int m, int n)
487{
488 TCL_TRAS(a, s, b, m, n)
489 return b;
490} /* tras_ */
491
492#undef TCL_TRAS
493
494////////////////////////////////////////////////////////////////////////////////
495
496#define TCL_TRASAT(a, s, r__, m, n) \
497 int imax, k; \
498 int ia, mn, ir, is, iaa; \
499 double sum; \
500 --r__; --s; --a; \
501 imax = (m * m + m) / 2; \
502 vzero(&r__[1], imax); \
503 mn = m * n; \
504 int ind = 0; \
505 int i__ = 0; \
506 do { \
507 ind += i__; \
508 ia = 0; ir = 0; \
509 do { \
510 is = ind; \
511 sum = 0.; k = 0; \
512 do { \
513 if (k > i__) is += k; \
514 else ++is; \
515 ++ia; \
516 sum += s[is] * a[ia]; \
517 ++k; \
518 } while (k < n); \
519 iaa = i__ + 1; \
520 do { \
521 ++ir; \
522 r__[ir] += sum * a[iaa];\
523 iaa += n; \
524 } while (iaa <= ia); \
525 } while (ia < mn); \
526 ++i__; \
527 } while (i__ < n); \
528 ++r__;
529
530////////////////////////////////////////////////////////////////////////////////
531/// Transformation of Symmetric Matrix
532///
533/// CERN PROGLIB# F112 TRASAT .VERSION KERNFOR 4.15 861204
534/// ORIG. 18/12/74 WH
535/// trasat.F -- translated by f2c (version 19970219).
536///
537/// See original documentation of CERNLIB package
538/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
539
540float *TCL::trasat(const float *a, const float *s, float *r__, int m, int n)
541{
542 TCL_TRASAT(a, s, r__, m, n)
543 return r__;
544} /* trasat_ */
545
546////////////////////////////////////////////////////////////////////////////////
547/// Transformation of Symmetric Matrix
548///
549/// CERN PROGLIB# F112 TRASAT .VERSION KERNFOR 4.15 861204
550/// ORIG. 18/12/74 WH
551/// trasat.F -- translated by f2c (version 19970219).
552///
553/// See original documentation of CERNLIB package
554/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
555
556double *TCL::trasat(const double *a, const double *s, double *r__, int m, int n)
557{
558 TCL_TRASAT(a, s, r__, m, n)
559 return r__;
560} /* trasat_ */
561
562////////////////////////////////////////////////////////////////////////////////
563/// Transformation of Symmetric Matrix
564///
565/// CERN PROGLIB# F112 TRASAT .VERSION KERNFOR 4.15 861204
566/// ORIG. 18/12/74 WH
567/// trasat.F -- translated by f2c (version 19970219).
568///
569/// See original documentation of CERNLIB package
570/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
571
572float *TCL::trasat(const double *a, const float *s, float *r__, int m, int n)
573{
574 TCL_TRASAT(a, s, r__, m, n)
575 return r__;
576} /* trasat_ */
577
578#undef TCL_TRASAT
579
580////////////////////////////////////////////////////////////////////////////////
581/// trata.F -- translated by f2c (version 19970219).
582///
583/// CERN PROGLIB# F112 TRATA .VERSION KERNFOR 4.15 861204
584/// ORIG. 18/12/74 WH
585///
586/// See original documentation of CERNLIB package
587/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
588
589float *TCL::trata(const float *a, float *r__, int m, int n)
590{
591 /* Local variables */
592 int i__, j, ia, mn, ir, iat;
593 double sum;
594
595 /* Parameter adjustments */
596 --r__; --a;
597
598 /* Function Body */
599 mn = m * n;
600 ir = 0;
601
602 for (i__ = 1; i__ <= m; ++i__) {
603 for (j = 1; j <= i__; ++j) {
604 ia = i__;
605 iat = j;
606 sum = 0.;
607 do {
608 sum += a[ia] * a[iat];
609 ia += m;
610 iat += m;
611 } while (ia <= mn);
612 ++ir;
613 r__[ir] = sum;
614 }
615 }
616 ++r__;
617 return r__;
618} /* trata_ */
619
620////////////////////////////////////////////////////////////////////////////////
621/// trats.F -- translated by f2c (version 19970219).
622///
623/// CERN PROGLIB# F112 TRATS .VERSION KERNFOR 4.15 861204
624/// ORIG. 18/12/74 WH
625///
626/// See original documentation of CERNLIB package
627/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
628
629float *TCL::trats(const float *a, const float *s, float *b, int m, int n)
630{
631 /* Local variables */
632 int inds, i__, j, k, ia, ib, is;
633 double sum;
634
635 /* Parameter adjustments */
636 --b; --s; --a;
637
638 /* Function Body */
639 ib = 0; inds = 0; i__ = 0;
640 do {
641 inds += i__;
642 ib = i__ + 1;
643
644 for (j = 1; j <= m; ++j) {
645 ia = j;
646 is = inds;
647 sum = 0.;
648 k = 0;
649
650 do {
651 if (k > i__) is += k;
652 else ++is;
653 sum += a[ia] * s[is];
654 ia += m;
655 ++k;
656 } while (k < n);
657
658 b[ib] = sum;
659 ib += n;
660 }
661 ++i__;
662 } while (i__ < n);
663 ++b;
664 return b;
665} /* trats_ */
666
667
668////////////////////////////////////////////////////////////////////////////////
669/// tratsa.F -- translated by f2c (version 19970219).
670///
671/// CERN PROGLIB# F112 TRATSA .VERSION KERNFOR 4.15 861204
672/// ORIG. 18/12/74 WH
673///
674/// See original documentation of CERNLIB package
675/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
676
677float *TCL::tratsa(const float *a, const float *s, float *r__, int m, int n)
678{
679 /* Local variables */
680 int imax, i__, j, k;
681 int ia, ir, is, iaa, ind;
682 double sum;
683
684 /* Parameter adjustments */
685 --r__; --s; --a;
686
687 /* Function Body */
688 imax = (m * m + m) / 2;
689 vzero(&r__[1], imax);
690 ind = 0;
691 i__ = 0;
692
693 do {
694 ind += i__;
695 ir = 0;
696
697 for (j = 1; j <= m; ++j) {
698 is = ind;
699 ia = j;
700 sum = 0.;
701 k = 0;
702
703 do {
704 if (k > i__) is += k;
705 else ++is;
706 sum += s[is] * a[ia];
707 ia += m;
708 ++k;
709 } while (k < n);
710 iaa = i__ * m;
711
712 for (k = 1; k <= j; ++k) {
713 ++iaa;
714 ++ir;
715 r__[ir] += sum * a[iaa];
716 }
717 }
718 ++i__;
719 } while (i__ < n);
720 ++r__;
721 return r__;
722} /* tratsa_ */
723
724////////////////////////////////////////////////////////////////////////////////
725/// trchlu.F -- translated by f2c (version 19970219).
726///
727/// CERN PROGLIB# F112 TRCHLU .VERSION KERNFOR 4.16 870601
728/// ORIG. 18/12/74 W.HART
729///
730/// See original documentation of CERNLIB package
731/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
732
733float *TCL::trchlu(const float *a, float *b, int n)
734{
735 /* Local variables */
736 int ipiv, kpiv, i__, j;
737 double r__, dc;
738 int id, kd;
739 double sum;
740
741 /* Parameter adjustments */
742 --b; --a;
743
744 /* Function Body */
745 ipiv = 0;
746
747 i__ = 0;
748
749 do {
750 ++i__;
751 ipiv += i__;
752 kpiv = ipiv;
753 r__ = a[ipiv];
754
755 for (j = i__; j <= n; ++j) {
756 sum = 0.;
757 if (i__ == 1) goto L40;
758 if (r__ == 0.) goto L42;
759 id = ipiv - i__ + 1;
760 kd = kpiv - i__ + 1;
761
762 do {
763 sum += b[kd] * b[id];
764 ++kd; ++id;
765 } while (id < ipiv);
766
767L40:
768 sum = a[kpiv] - sum;
769L42:
770 if (j != i__) b[kpiv] = sum * r__;
771 else {
772 dc = TMath::Sqrt(sum);
773 b[kpiv] = dc;
774 if (r__ > 0.) r__ = 1. / dc;
775 }
776 kpiv += j;
777 }
778
779 } while (i__ < n);
780 ++b;
781 return b;
782} /* trchlu_ */
783
784////////////////////////////////////////////////////////////////////////////////
785/// trchul.F -- translated by f2c (version 19970219).
786///
787/// CERN PROGLIB# F112 TRCHUL .VERSION KERNFOR 4.16 870601
788/// ORIG. 18/12/74 WH
789///
790/// See original documentation of CERNLIB package
791/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
792
793float *TCL::trchul(const float *a, float *b, int n)
794{
795 /* Local variables */
796 int ipiv, kpiv, i__;
797 double r__;
798 int nTep;
799 double dc;
800 int id, kd;
801 double sum;
802
803 /* Parameter adjustments */
804 --b; --a;
805
806 /* Function Body */
807 kpiv = (n * n + n) / 2;
808
809 i__ = n;
810 do {
811 ipiv = kpiv;
812 r__ = a[ipiv];
813
814 do {
815 sum = 0.;
816 if (i__ == n) goto L40;
817 if (r__ == 0.) goto L42;
818 id = ipiv;
819 kd = kpiv;
820 nTep = i__;
821
822 do {
823 kd += nTep;
824 id += nTep;
825 ++nTep;
826 sum += b[id] * b[kd];
827 } while (nTep < n);
828
829L40:
830 sum = a[kpiv] - sum;
831L42:
832 if (kpiv < ipiv) b[kpiv] = sum * r__;
833 else {
834 dc = TMath::Sqrt(sum);
835 b[kpiv] = dc;
836 if (r__ > 0.) r__ = 1. / dc;
837 }
838 --kpiv;
839 } while (kpiv > ipiv - i__);
840
841 --i__;
842 } while (i__ > 0);
843
844 ++b;
845 return b;
846} /* trchul_ */
847
848////////////////////////////////////////////////////////////////////////////////
849/// trinv.F -- translated by f2c (version 19970219).
850///
851/// CERN PROGLIB# F112 TRINV .VERSION KERNFOR 4.15 861204
852/// ORIG. 18/12/74 WH
853///
854/// See original documentation of CERNLIB package
855/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
856
857float *TCL::trinv(const float *t, float *s, int n)
858{
859 int lhor, ipiv, lver, j;
860 double sum = 0;
861 double r__ = 0;
862 int mx, ndTep, ind;
863
864 /* Parameter adjustments */
865 --s; --t;
866
867 /* Function Body */
868 mx = (n * n + n) / 2;
869 ipiv = mx;
870
871 int i = n;
872 do {
873 r__ = 0.;
874 if (t[ipiv] > 0.) r__ = 1. / t[ipiv];
875 s[ipiv] = r__;
876 ndTep = n;
877 ind = mx - n + i;
878
879 while (ind != ipiv) {
880 sum = 0.;
881 if (r__ != 0.) {
882 lhor = ipiv;
883 lver = ind;
884 j = i;
885
886 do {
887 lhor += j;
888 ++lver;
889 sum += t[lhor] * s[lver];
890 ++j;
891 } while (lhor < ind);
892 }
893 s[ind] = -sum * r__;
894 --ndTep;
895 ind -= ndTep;
896 }
897
898 ipiv -= i;
899 --i;
900 } while (i > 0);
901
902 ++s;
903 return s;
904} /* trinv_ */
905
906////////////////////////////////////////////////////////////////////////////////
907/// trla.F -- translated by f2c (version 19970219).
908///
909/// CERN PROGLIB# F112 TRLA .VERSION KERNFOR 4.15 861204
910/// ORIG. 18/12/74 WH
911///
912/// See original documentation of CERNLIB package
913/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
914
915float *TCL::trla(const float *u, const float *a, float *b, int m, int n)
916{
917 int ipiv, ia, ib, iu;
918 double sum;
919
920 /* Parameter adjustments */
921 --b; --a; --u;
922
923 /* Function Body */
924 ib = m * n;
925 ipiv = (m * m + m) / 2;
926
927 do {
928 do {
929 ia = ib;
930 iu = ipiv;
931
932 sum = 0.;
933 do {
934 sum += a[ia] * u[iu];
935 --iu;
936 ia -= n;
937 } while (ia > 0);
938
939 b[ib] = sum;
940 --ib;
941 } while (ia > 1 - n);
942
943 ipiv = iu;
944 } while (iu > 0);
945
946 ++b;
947 return b;
948} /* trla_ */
949
950////////////////////////////////////////////////////////////////////////////////
951/// trlta.F -- translated by f2c (version 19970219).
952///
953/// CERN PROGLIB# F112 TRLTA .VERSION KERNFOR 4.15 861204
954/// ORIG. 18/12/74 WH
955///
956/// See original documentation of CERNLIB package
957/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
958
959float *TCL::trlta(const float *u, const float *a, float *b, int m, int n)
960{
961 int ipiv, mxpn, i__, nTep, ia, ib, iu, mx;
962 double sum;
963
964 /* Parameter adjustments */
965 --b; --a; --u;
966
967 /* Function Body */
968 ipiv = 0;
969 mx = m * n;
970 mxpn = mx + n;
971 ib = 0;
972
973 i__ = 0;
974 do {
975 ++i__;
976 ipiv += i__;
977
978 do {
979 iu = ipiv;
980 nTep = i__;
981 ++ib;
982 ia = ib;
983
984 sum = 0.;
985 do {
986 sum += a[ia] * u[iu];
987 ia += n;
988 iu += nTep;
989 ++nTep;
990 } while (ia <= mx);
991
992 b[ib] = sum;
993 } while (ia < mxpn);
994
995 } while (i__ < m);
996
997 ++b;
998 return b;
999} /* trlta_ */
1000
1001////////////////////////////////////////////////////////////////////////////////
1002/// trpck.F -- translated by f2c (version 19970219).
1003///
1004/// CERN PROGLIB# F112 TRPCK .VERSION KERNFOR 2.08 741218
1005/// ORIG. 18/12/74 WH
1006///
1007/// See original documentation of CERNLIB package
1008/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1009
1010float *TCL::trpck(const float *s, float *u, int n)
1011{
1012 int i__, ia, ind, ipiv;
1013
1014 /* Parameter adjustments */
1015 --u; --s;
1016
1017 /* Function Body */
1018 ia = 0;
1019 ind = 0;
1020 ipiv = 0;
1021
1022 for (i__ = 1; i__ <= n; ++i__) {
1023 ipiv += i__;
1024 do {
1025 ++ia;
1026 ++ind;
1027 u[ind] = s[ia];
1028 } while (ind < ipiv);
1029 ia = ia + n - i__;
1030 }
1031
1032 ++u;
1033 return u;
1034} /* trpck_ */
1035
1036////////////////////////////////////////////////////////////////////////////////
1037/// trqsq.F -- translated by f2c (version 19970219).
1038///
1039/// CERN PROGLIB# F112 TRQSQ .VERSION KERNFOR 4.15 861204
1040/// ORIG. 18/12/74 WH
1041///
1042/// See original documentation of CERNLIB package
1043/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1044
1045float *TCL::trqsq(const float *q, const float *s, float *r__, int m)
1046{
1047 int indq, inds, imax, i__, j, k, l;
1048 int iq, ir, is, iqq;
1049 double sum;
1050
1051 /* Parameter adjustments */
1052 --r__; --s; --q;
1053
1054 /* Function Body */
1055 imax = (m * m + m) / 2;
1056 vzero(&r__[1], imax);
1057 inds = 0;
1058 i__ = 0;
1059
1060 do {
1061 inds += i__;
1062 ir = 0;
1063 indq = 0;
1064 j = 0;
1065
1066 do {
1067 indq += j;
1068 is = inds;
1069 iq = indq;
1070 sum = (float)0.;
1071 k = 0;
1072
1073 do {
1074 if (k > i__) is += k;
1075 else ++is;
1076
1077 if (k > j) iq += k;
1078 else ++iq;
1079
1080 sum += s[is] * q[iq];
1081 ++k;
1082 } while (k < m);
1083 iqq = inds;
1084 l = 0;
1085
1086 do {
1087 ++ir;
1088 if (l > i__) iqq += l;
1089 else ++iqq;
1090 r__[ir] += q[iqq] * sum;
1091 ++l;
1092 } while (l <= j);
1093 ++j;
1094 } while (j < m);
1095 ++i__;
1096 } while (i__ < m);
1097
1098 ++r__;
1099 return r__;
1100} /* trqsq_ */
1101
1102////////////////////////////////////////////////////////////////////////////////
1103/// trsa.F -- translated by f2c (version 19970219).
1104///
1105/// CERN PROGLIB# F112 TRSA .VERSION KERNFOR 4.15 861204
1106/// ORIG. 18/12/74 WH
1107///
1108/// See original documentation of CERNLIB package
1109/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1110
1111float *TCL::trsa(const float *s, const float *a, float *b, int m, int n)
1112{
1113 /* Local variables */
1114 int inds, i__, j, k, ia, ib, is;
1115 double sum;
1116
1117 /* Parameter adjustments */
1118 --b; --a; --s;
1119
1120 /* Function Body */
1121 inds = 0;
1122 ib = 0;
1123 i__ = 0;
1124
1125 do {
1126 inds += i__;
1127
1128 for (j = 1; j <= n; ++j) {
1129 ia = j;
1130 is = inds;
1131 sum = 0.;
1132 k = 0;
1133
1134 do {
1135 if (k > i__) is += k;
1136 else ++is;
1137 sum += s[is] * a[ia];
1138 ia += n;
1139 ++k;
1140 } while (k < m);
1141 ++ib;
1142 b[ib] = sum;
1143 }
1144 ++i__;
1145 } while (i__ < m);
1146
1147 ++b;
1148 return b;
1149} /* trsa_ */
1150
1151////////////////////////////////////////////////////////////////////////////////
1152/// trsinv.F -- translated by f2c (version 19970219).
1153///
1154/// CERN PROGLIB# F112 TRSINV .VERSION KERNFOR 2.08 741218
1155/// ORIG. 18/12/74 WH
1156///
1157/// See original documentation of CERNLIB package
1158/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1159
1160float *TCL::trsinv(const float *g, float *gi, int n)
1161{
1162 /* Function Body */
1163 trchlu(g, gi, n);
1164 trinv(gi, gi, n);
1165 return trsmul(gi, gi, n);
1166} /* trsinv_ */
1167
1168////////////////////////////////////////////////////////////////////////////////
1169/// trsmlu.F -- translated by f2c (version 19970219).
1170///
1171/// CERN PROGLIB# F112 TRSMLU .VERSION KERNFOR 4.15 861204
1172/// ORIG. 18/12/74 WH
1173///
1174/// See original documentation of CERNLIB package
1175/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1176
1177float *TCL::trsmlu(const float *u, float *s, int n)
1178{
1179 /* Local variables */
1180 int lhor, lver, i__, k, l, ind;
1181 double sum;
1182
1183 /* Parameter adjustments */
1184 --s; --u;
1185
1186 /* Function Body */
1187 ind = (n * n + n) / 2;
1188
1189 for (i__ = 1; i__ <= n; ++i__) {
1190 lver = ind;
1191
1192 for (k = i__; k <= n; ++k,--ind) {
1193 lhor = ind; sum = 0.;
1194 for (l = k; l <= n; ++l,--lver,--lhor)
1195 sum += u[lver] * u[lhor];
1196 s[ind] = sum;
1197 }
1198 }
1199 ++s;
1200 return s;
1201} /* trsmlu_ */
1202
1203////////////////////////////////////////////////////////////////////////////////
1204/// trsmul.F -- translated by f2c (version 19970219).
1205///
1206/// CERN PROGLIB# F112 TRSMUL .VERSION KERNFOR 4.15 861204
1207/// ORIG. 18/12/74 WH
1208///
1209/// See original documentation of CERNLIB package
1210/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1211
1212float *TCL::trsmul(const float *g, float *gi, int n)
1213{
1214 /* Local variables */
1215 int lhor, lver, lpiv, i__, j, k, ind;
1216 double sum;
1217
1218 /* Parameter adjustments */
1219 --gi; --g;
1220
1221 /* Function Body */
1222 ind = 1;
1223 lpiv = 0;
1224 for (i__ = 1; i__ <= n; ++i__) {
1225 lpiv += i__;
1226 for (j = 1; j <= i__; ++j,++ind) {
1227 lver = lpiv;
1228 lhor = ind;
1229 sum = 0.;
1230 for (k = i__; k <= n; lhor += k,lver += k,++k)
1231 sum += g[lver] * g[lhor];
1232 gi[ind] = sum;
1233 }
1234 }
1235 ++gi;
1236 return gi;
1237} /* trsmul_ */
1238
1239////////////////////////////////////////////////////////////////////////////////
1240/// trupck.F -- translated by f2c (version 19970219).
1241///
1242/// CERN PROGLIB# F112 TRUPCK .VERSION KERNFOR 2.08 741218
1243/// ORIG. 18/12/74 WH
1244///
1245/// See original documentation of CERNLIB package
1246/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1247
1248float *TCL::trupck(const float *u, float *s, int m)
1249{
1250 int i__, im, is, iu, iv, ih, m2;
1251
1252 /* Parameter adjustments */
1253 --s; --u;
1254
1255 /* Function Body */
1256 m2 = m * m;
1257 is = m2;
1258 iu = (m2 + m) / 2;
1259 i__ = m - 1;
1260
1261 do {
1262 im = i__ * m;
1263 do {
1264 s[is] = u[iu];
1265 --is;
1266 --iu;
1267 } while (is > im);
1268 is = is - m + i__;
1269 --i__;
1270 } while (i__ >= 0);
1271
1272 is = 1;
1273 do {
1274 iv = is;
1275 ih = is;
1276 while (1) {
1277 iv += m;
1278 ++ih;
1279 if (iv > m2) break;
1280 s[ih] = s[iv];
1281 }
1282 is = is + m + 1;
1283 } while (is < m2);
1284
1285 ++s;
1286 return s;
1287} /* trupck_ */
1288
1289////////////////////////////////////////////////////////////////////////////////
1290/// trsat.F -- translated by f2c (version 19970219).
1291///
1292/// CERN PROGLIB# F112 TRSAT .VERSION KERNFOR 4.15 861204
1293/// ORIG. 18/12/74 WH
1294///
1295/// See original documentation of CERNLIB package
1296/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1297
1298float *TCL::trsat(const float *s, const float *a, float *b, int m, int n)
1299{
1300
1301 /* Local variables */
1302 int inds, i__, j, k, ia, ib, is;
1303 double sum;
1304
1305 /* Parameter adjustments */
1306 --b; --a; --s;
1307
1308 /* Function Body */
1309 inds = 0;
1310 ib = 0;
1311 i__ = 0;
1312
1313 do {
1314 inds += i__;
1315 ia = 0;
1316
1317 for (j = 1; j <= n; ++j) {
1318 is = inds;
1319 sum = 0.;
1320 k = 0;
1321
1322 do {
1323 if (k > i__) is += k;
1324 else ++is;
1325 ++ia;
1326 sum += s[is] * a[ia];
1327 ++k;
1328 } while (k < m);
1329 ++ib;
1330 b[ib] = sum;
1331 }
1332 ++i__;
1333 } while (i__ < m);
1334
1335 ++b;
1336 return b;
1337} /* trsat_ */
1338
1339// ------ double
1340
1341////////////////////////////////////////////////////////////////////////////////
1342/// trata.F -- translated by f2c (version 19970219).
1343///
1344/// CERN PROGLIB# F112 TRATA .VERSION KERNFOR 4.15 861204
1345/// ORIG. 18/12/74 WH
1346///
1347/// See original documentation of CERNLIB package
1348/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1349
1350double *TCL::trata(const double *a, double *r__, int m, int n)
1351{
1352
1353 /* Local variables */
1354 int i__, j, ia, mn, ir, iat;
1355 double sum;
1356
1357 /* Parameter adjustments */
1358 --r__; --a;
1359
1360 /* Function Body */
1361 mn = m * n;
1362 ir = 0;
1363
1364 for (i__ = 1; i__ <= m; ++i__) {
1365
1366 for (j = 1; j <= i__; ++j) {
1367 ia = i__;
1368 iat = j;
1369
1370 sum = (double)0.;
1371 do {
1372 sum += a[ia] * a[iat];
1373 ia += m;
1374 iat += m;
1375 } while (ia <= mn);
1376 ++ir;
1377 r__[ir] = sum;
1378 }
1379 }
1380
1381 return 0;
1382} /* trata_ */
1383
1384////////////////////////////////////////////////////////////////////////////////
1385/// trats.F -- translated by f2c (version 19970219).
1386///
1387/// CERN PROGLIB# F112 TRATS .VERSION KERNFOR 4.15 861204
1388/// ORIG. 18/12/74 WH
1389///
1390/// See original documentation of CERNLIB package
1391/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1392
1393double *TCL::trats(const double *a, const double *s, double *b, int m, int n)
1394{
1395 /* Local variables */
1396 int inds, i__, j, k, ia, ib, is;
1397 double sum;
1398
1399 /* Parameter adjustments */
1400 --b; --s; --a;
1401
1402 /* Function Body */
1403 ib = 0; inds = 0; i__ = 0;
1404
1405 do {
1406 inds += i__;
1407 ib = i__ + 1;
1408
1409 for (j = 1; j <= m; ++j) {
1410 ia = j;
1411 is = inds;
1412 sum = (double)0.;
1413 k = 0;
1414
1415 do {
1416 if (k > i__) is += k;
1417 else ++is;
1418 sum += a[ia] * s[is];
1419 ia += m;
1420 ++k;
1421 } while (k < n);
1422
1423 b[ib] = sum;
1424 ib += n;
1425 }
1426 ++i__;
1427 } while (i__ < n);
1428
1429 return 0;
1430} /* trats_ */
1431
1432////////////////////////////////////////////////////////////////////////////////
1433/// tratsa.F -- translated by f2c (version 19970219).
1434///
1435/// CERN PROGLIB# F112 TRATSA .VERSION KERNFOR 4.15 861204
1436/// ORIG. 18/12/74 WH
1437///
1438/// See original documentation of CERNLIB package
1439/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1440
1441double *TCL::tratsa(const double *a, const double *s, double *r__, int m, int n)
1442{
1443 /* Local variables */
1444 int imax, i__, j, k;
1445 int ia, ir, is, iaa, ind;
1446 double sum;
1447
1448 /* Parameter adjustments */
1449 --r__; --s; --a;
1450
1451 /* Function Body */
1452 imax = (m * m + m) / 2;
1453 vzero(&r__[1], imax);
1454 ind = 0;
1455 i__ = 0;
1456
1457 do {
1458 ind += i__;
1459 ir = 0;
1460
1461 for (j = 1; j <= m; ++j) {
1462 is = ind;
1463 ia = j;
1464 sum = (double)0.;
1465 k = 0;
1466
1467 do {
1468 if (k > i__) is += k;
1469 else ++is;
1470 sum += s[is] * a[ia];
1471 ia += m;
1472 ++k;
1473 } while (k < n);
1474 iaa = i__ * m;
1475
1476 for (k = 1; k <= j; ++k) {
1477 ++iaa;
1478 ++ir;
1479 r__[ir] += sum * a[iaa];
1480 }
1481 }
1482 ++i__;
1483 } while (i__ < n);
1484
1485 return 0;
1486} /* tratsa_ */
1487
1488////////////////////////////////////////////////////////////////////////////////
1489/// trchlu.F -- translated by f2c (version 19970219).
1490///
1491/// CERN PROGLIB# F112 TRCHLU .VERSION KERNFOR 4.16 870601
1492/// ORIG. 18/12/74 W.HART
1493///
1494/// See original documentation of CERNLIB package
1495/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1496
1497double *TCL::trchlu(const double *a, double *b, int n)
1498{
1499 /* Local variables */
1500 int ipiv, kpiv, i__, j;
1501 double r__, dc;
1502 int id, kd;
1503 double sum;
1504
1505 /* Parameter adjustments */
1506 --b; --a;
1507
1508 /* Function Body */
1509 ipiv = 0;
1510
1511 i__ = 0;
1512
1513 do {
1514 ++i__;
1515 ipiv += i__;
1516 kpiv = ipiv;
1517 r__ = a[ipiv];
1518
1519 for (j = i__; j <= n; ++j) {
1520 sum = 0.;
1521 if (i__ == 1) goto L40;
1522 if (r__ == 0.) goto L42;
1523 id = ipiv - i__ + 1;
1524 kd = kpiv - i__ + 1;
1525
1526 do {
1527 sum += b[kd] * b[id];
1528 ++kd; ++id;
1529 } while (id < ipiv);
1530
1531L40:
1532 sum = a[kpiv] - sum;
1533L42:
1534 if (j != i__) b[kpiv] = sum * r__;
1535 else {
1536 dc = TMath::Sqrt(sum);
1537 b[kpiv] = dc;
1538 if (r__ > 0.) r__ = (double)1. / dc;
1539 }
1540 kpiv += j;
1541 }
1542
1543 } while (i__ < n);
1544
1545 return 0;
1546} /* trchlu_ */
1547
1548////////////////////////////////////////////////////////////////////////////////
1549/// trchul.F -- translated by f2c (version 19970219).
1550///
1551/// CERN PROGLIB# F112 TRCHUL .VERSION KERNFOR 4.16 870601
1552/// ORIG. 18/12/74 WH
1553///
1554/// See original documentation of CERNLIB package
1555/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1556
1557double *TCL::trchul(const double *a, double *b, int n)
1558{
1559 /* Local variables */
1560 int ipiv, kpiv, i__;
1561 double r__;
1562 int nTep;
1563 double dc;
1564 int id, kd;
1565 double sum;
1566
1567 /* Parameter adjustments */
1568 --b; --a;
1569
1570 /* Function Body */
1571 kpiv = (n * n + n) / 2;
1572
1573 i__ = n;
1574 do {
1575 ipiv = kpiv;
1576 r__ = a[ipiv];
1577
1578 do {
1579 sum = 0.;
1580 if (i__ == n) goto L40;
1581 if (r__ == (double)0.) goto L42;
1582 id = ipiv;
1583 kd = kpiv;
1584 nTep = i__;
1585
1586 do {
1587 kd += nTep;
1588 id += nTep;
1589 ++nTep;
1590 sum += b[id] * b[kd];
1591 } while (nTep < n);
1592
1593L40:
1594 sum = a[kpiv] - sum;
1595L42:
1596 if (kpiv < ipiv) b[kpiv] = sum * r__;
1597 else {
1598 dc = TMath::Sqrt(sum);
1599 b[kpiv] = dc;
1600 if (r__ > (double)0.) r__ = (double)1. / dc;
1601 }
1602 --kpiv;
1603 } while (kpiv > ipiv - i__);
1604
1605 --i__;
1606 } while (i__ > 0);
1607
1608 return 0;
1609} /* trchul_ */
1610
1611////////////////////////////////////////////////////////////////////////////////
1612/// trinv.F -- translated by f2c (version 19970219).
1613///
1614/// CERN PROGLIB# F112 TRINV .VERSION KERNFOR 4.15 861204
1615/// ORIG. 18/12/74 WH
1616///
1617/// See original documentation of CERNLIB package
1618/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1619
1620double *TCL::trinv(const double *t, double *s, int n)
1621{
1622 int lhor, ipiv, lver, j;
1623 double r__;
1624 int mx, ndTep, ind;
1625 double sum;
1626
1627 /* Parameter adjustments */
1628 --s; --t;
1629
1630 /* Function Body */
1631 mx = (n * n + n) / 2;
1632 ipiv = mx;
1633
1634 int i = n;
1635 do {
1636 r__ = 0.;
1637 if (t[ipiv] > 0.) r__ = (double)1. / t[ipiv];
1638 s[ipiv] = r__;
1639 ndTep = n;
1640 ind = mx - n + i;
1641
1642 while (ind != ipiv) {
1643 sum = 0.;
1644 if (r__ != 0.) {
1645 lhor = ipiv;
1646 lver = ind;
1647 j = i;
1648
1649 do {
1650 lhor += j;
1651 ++lver;
1652 sum += t[lhor] * s[lver];
1653 ++j;
1654 } while (lhor < ind);
1655 }
1656 s[ind] = -sum * r__;
1657 --ndTep;
1658 ind -= ndTep;
1659 }
1660
1661 ipiv -= i;
1662 --i;
1663 } while (i > 0);
1664
1665 return 0;
1666} /* trinv_ */
1667
1668////////////////////////////////////////////////////////////////////////////////
1669/// trla.F -- translated by f2c (version 19970219).
1670///
1671/// CERN PROGLIB# F112 TRLA .VERSION KERNFOR 4.15 861204
1672/// ORIG. 18/12/74 WH
1673///
1674/// See original documentation of CERNLIB package
1675/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1676
1677double *TCL::trla(const double *u, const double *a, double *b, int m, int n)
1678{
1679 int ipiv, ia, ib, iu;
1680 double sum;
1681
1682 /* Parameter adjustments */
1683 --b; --a; --u;
1684
1685 /* Function Body */
1686 ib = m * n;
1687 ipiv = (m * m + m) / 2;
1688
1689 do {
1690 do {
1691 ia = ib;
1692 iu = ipiv;
1693
1694 sum = 0.;
1695 do {
1696 sum += a[ia] * u[iu];
1697 --iu;
1698 ia -= n;
1699 } while (ia > 0);
1700
1701 b[ib] = sum;
1702 --ib;
1703 } while (ia > 1 - n);
1704
1705 ipiv = iu;
1706 } while (iu > 0);
1707
1708 return 0;
1709} /* trla_ */
1710
1711////////////////////////////////////////////////////////////////////////////////
1712/// trlta.F -- translated by f2c (version 19970219).
1713///
1714/// CERN PROGLIB# F112 TRLTA .VERSION KERNFOR 4.15 861204
1715/// ORIG. 18/12/74 WH
1716///
1717/// See original documentation of CERNLIB package
1718/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1719
1720double *TCL::trlta(const double *u, const double *a, double *b, int m, int n)
1721{
1722 int ipiv, mxpn, i__, nTep, ia, ib, iu, mx;
1723 double sum;
1724
1725 /* Parameter adjustments */
1726 --b; --a; --u;
1727
1728 /* Function Body */
1729 ipiv = 0;
1730 mx = m * n;
1731 mxpn = mx + n;
1732 ib = 0;
1733
1734 i__ = 0;
1735 do {
1736 ++i__;
1737 ipiv += i__;
1738
1739 do {
1740 iu = ipiv;
1741 nTep = i__;
1742 ++ib;
1743 ia = ib;
1744
1745 sum = 0.;
1746 do {
1747 sum += a[ia] * u[iu];
1748 ia += n;
1749 iu += nTep;
1750 ++nTep;
1751 } while (ia <= mx);
1752
1753 b[ib] = sum;
1754 } while (ia < mxpn);
1755
1756 } while (i__ < m);
1757
1758 return 0;
1759} /* trlta_ */
1760
1761////////////////////////////////////////////////////////////////////////////////
1762/// trpck.F -- translated by f2c (version 19970219).
1763///
1764/// CERN PROGLIB# F112 TRPCK .VERSION KERNFOR 2.08 741218
1765/// ORIG. 18/12/74 WH
1766///
1767/// See original documentation of CERNLIB package
1768/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1769
1770double *TCL::trpck(const double *s, double *u, int n)
1771{
1772 int i__, ia, ind, ipiv;
1773
1774 /* Parameter adjustments */
1775 --u; --s;
1776
1777 /* Function Body */
1778 ia = 0;
1779 ind = 0;
1780 ipiv = 0;
1781
1782 for (i__ = 1; i__ <= n; ++i__) {
1783 ipiv += i__;
1784 do {
1785 ++ia;
1786 ++ind;
1787 u[ind] = s[ia];
1788 } while (ind < ipiv);
1789 ia = ia + n - i__;
1790 }
1791
1792 return 0;
1793} /* trpck_ */
1794
1795////////////////////////////////////////////////////////////////////////////////
1796/// trqsq.F -- translated by f2c (version 19970219).
1797///
1798/// CERN PROGLIB# F112 TRQSQ .VERSION KERNFOR 4.15 861204
1799/// ORIG. 18/12/74 WH
1800///
1801/// See original documentation of CERNLIB package
1802/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1803
1804double *TCL::trqsq(const double *q, const double *s, double *r__, int m)
1805{
1806 int indq, inds, imax, i__, j, k, l;
1807 int iq, ir, is, iqq;
1808 double sum;
1809
1810 /* Parameter adjustments */
1811 --r__; --s; --q;
1812
1813 /* Function Body */
1814 imax = (m * m + m) / 2;
1815 vzero(&r__[1], imax);
1816 inds = 0;
1817 i__ = 0;
1818
1819 do {
1820 inds += i__;
1821 ir = 0;
1822 indq = 0;
1823 j = 0;
1824
1825 do {
1826 indq += j;
1827 is = inds;
1828 iq = indq;
1829 sum = 0.;
1830 k = 0;
1831
1832 do {
1833 if (k > i__) is += k;
1834 else ++is;
1835
1836 if (k > j) iq += k;
1837 else ++iq;
1838
1839 sum += s[is] * q[iq];
1840 ++k;
1841 } while (k < m);
1842 iqq = inds;
1843 l = 0;
1844
1845 do {
1846 ++ir;
1847 if (l > i__) iqq += l;
1848 else ++iqq;
1849 r__[ir] += q[iqq] * sum;
1850 ++l;
1851 } while (l <= j);
1852 ++j;
1853 } while (j < m);
1854 ++i__;
1855 } while (i__ < m);
1856
1857 return 0;
1858} /* trqsq_ */
1859
1860////////////////////////////////////////////////////////////////////////////////
1861/// trsa.F -- translated by f2c (version 19970219).
1862///
1863/// CERN PROGLIB# F112 TRSA .VERSION KERNFOR 4.15 861204
1864/// ORIG. 18/12/74 WH
1865///
1866/// See original documentation of CERNLIB package
1867/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1868
1869double *TCL::trsa(const double *s, const double *a, double *b, int m, int n)
1870{
1871 /* Local variables */
1872 int inds, i__, j, k, ia, ib, is;
1873 double sum;
1874
1875 /* Parameter adjustments */
1876 --b; --a; --s;
1877
1878 /* Function Body */
1879 inds = 0;
1880 ib = 0;
1881 i__ = 0;
1882
1883 do {
1884 inds += i__;
1885
1886 for (j = 1; j <= n; ++j) {
1887 ia = j;
1888 is = inds;
1889 sum = 0.;
1890 k = 0;
1891
1892 do {
1893 if (k > i__) is += k;
1894 else ++is;
1895 sum += s[is] * a[ia];
1896 ia += n;
1897 ++k;
1898 } while (k < m);
1899 ++ib;
1900 b[ib] = sum;
1901 }
1902 ++i__;
1903 } while (i__ < m);
1904
1905 return 0;
1906} /* trsa_ */
1907
1908////////////////////////////////////////////////////////////////////////////////
1909/// trsinv.F -- translated by f2c (version 19970219).
1910///
1911/// CERN PROGLIB# F112 TRSINV .VERSION KERNFOR 2.08 741218
1912/// ORIG. 18/12/74 WH
1913///
1914/// See original documentation of CERNLIB package
1915/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1916
1917double *TCL::trsinv(const double *g, double *gi, int n)
1918{
1919
1920 /* Function Body */
1921 trchlu(g, gi, n);
1922 trinv(gi, gi, n);
1923 trsmul(gi, gi, n);
1924
1925 return 0;
1926} /* trsinv_ */
1927
1928////////////////////////////////////////////////////////////////////////////////
1929/// trsmlu.F -- translated by f2c (version 19970219).
1930///
1931/// CERN PROGLIB# F112 TRSMLU .VERSION KERNFOR 4.15 861204
1932/// ORIG. 18/12/74 WH
1933///
1934/// See original documentation of CERNLIB package
1935/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1936
1937double *TCL::trsmlu(const double *u, double *s, int n)
1938{
1939
1940 /* Local variables */
1941 int lhor, lver, i__, k, l, ind;
1942 double sum;
1943
1944 /* Parameter adjustments */
1945 --s; --u;
1946
1947 /* Function Body */
1948 ind = (n * n + n) / 2;
1949
1950 for (i__ = 1; i__ <= n; ++i__) {
1951 lver = ind;
1952
1953 for (k = i__; k <= n; ++k,--ind) {
1954 lhor = ind; sum = 0.;
1955 for (l = k; l <= n; ++l,--lver,--lhor)
1956 sum += u[lver] * u[lhor];
1957 s[ind] = sum;
1958 }
1959 }
1960
1961 return 0;
1962} /* trsmlu_ */
1963
1964////////////////////////////////////////////////////////////////////////////////
1965/// trsmul.F -- translated by f2c (version 19970219).
1966///
1967/// CERN PROGLIB# F112 TRSMUL .VERSION KERNFOR 4.15 861204
1968/// ORIG. 18/12/74 WH
1969///
1970/// See original documentation of CERNLIB package
1971/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
1972
1973double *TCL::trsmul(const double *g, double *gi, int n)
1974{
1975
1976 /* Local variables */
1977 int lhor, lver, lpiv, i__, j, k, ind;
1978 double sum;
1979
1980 /* Parameter adjustments */
1981 --gi; --g;
1982
1983 /* Function Body */
1984 ind = 1;
1985 lpiv = 0;
1986 for (i__ = 1; i__ <= n; ++i__) {
1987 lpiv += i__;
1988 for (j = 1; j <= i__; ++j,++ind) {
1989 lver = lpiv;
1990 lhor = ind;
1991 sum = 0.;
1992 for (k = i__; k <= n;lhor += k,lver += k,++k)
1993 sum += g[lver] * g[lhor];
1994 gi[ind] = sum;
1995 }
1996 }
1997
1998 return 0;
1999} /* trsmul_ */
2000
2001////////////////////////////////////////////////////////////////////////////////
2002/// trupck.F -- translated by f2c (version 19970219).
2003///
2004/// CERN PROGLIB# F112 TRUPCK .VERSION KERNFOR 2.08 741218
2005/// ORIG. 18/12/74 WH
2006///
2007/// See original documentation of CERNLIB package
2008/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
2009
2010double *TCL::trupck(const double *u, double *s, int m)
2011{
2012 int i__, im, is, iu, iv, ih, m2;
2013
2014 /* Parameter adjustments */
2015 --s; --u;
2016
2017 /* Function Body */
2018 m2 = m * m;
2019 is = m2;
2020 iu = (m2 + m) / 2;
2021 i__ = m - 1;
2022
2023 do {
2024 im = i__ * m;
2025 do {
2026 s[is] = u[iu];
2027 --is;
2028 --iu;
2029 } while (is > im);
2030 is = is - m + i__;
2031 --i__;
2032 } while (i__ >= 0);
2033
2034 is = 1;
2035 do {
2036 iv = is;
2037 ih = is;
2038 while (1) {
2039 iv += m;
2040 ++ih;
2041 if (iv > m2) break;
2042 s[ih] = s[iv];
2043 }
2044 is = is + m + 1;
2045 } while (is < m2);
2046
2047 return 0;
2048} /* trupck_ */
2049
2050////////////////////////////////////////////////////////////////////////////////
2051/// trsat.F -- translated by f2c (version 19970219)
2052///
2053/// CERN PROGLIB# F112 TRSAT .VERSION KERNFOR 4.15 861204
2054/// ORIG. 18/12/74 WH
2055///
2056/// See original documentation of CERNLIB package
2057/// [F112](https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/f112/top.html)
2058
2059double *TCL::trsat(const double *s, const double *a, double *b, int m, int n)
2060{
2061 /* Local variables */
2062 int inds, i__, j, k, ia, ib, is;
2063 double sum;
2064
2065 /* Parameter adjustments */
2066 --b; --a; --s;
2067
2068 /* Function Body */
2069 inds = 0;
2070 ib = 0;
2071 i__ = 0;
2072
2073 do {
2074 inds += i__;
2075 ia = 0;
2076
2077 for (j = 1; j <= n; ++j) {
2078 is = inds;
2079 sum = 0.;
2080 k = 0;
2081
2082 do {
2083 if (k > i__) is += k;
2084 else ++is;
2085 ++ia;
2086 sum += s[is] * a[ia];
2087 ++k;
2088 } while (k < m);
2089 ++ib;
2090 b[ib] = sum;
2091 }
2092 ++i__;
2093 } while (i__ < m);
2094
2095 return 0;
2096} /* trsat_ */
2097
2098// ------------ Victor Perevoztchikov's addition
2099
2100////////////////////////////////////////////////////////////////////////////////
2101/// Linear Equations, Matrix Inversion
2102/// trsequ solves the matrix equation
2103///
2104/// SMX*x = B
2105///
2106/// which represents a system of m simultaneous linear equations with n right-hand sides:
2107/// SMX is an unpacked symmetric matrix (all elements) (m x m)
2108/// B is an unpacked matrix of right-hand sides (n x m)
2109///
2110
2111float *TCL::trsequ(float *smx, int m, float *b, int n)
2112{
2113 float *mem = new float[(m*(m+1))/2+m];
2114 float *v = mem;
2115 float *s = v+m;
2116 if (!b) n=0;
2117 TCL::trpck (smx,s ,m);
2118 TCL::trsinv(s ,s, m);
2119
2120 for (int i=0;i<n;i++) {
2121 TCL::trsa (s ,b+i*m, v, m, 1);
2122 TCL::ucopy (v ,b+i*m, m);}
2123 TCL::trupck(s ,smx, m);
2124 delete [] mem;
2125 return b;
2126}
2127
2128////////////////////////////////////////////////////////////////////////////////
2129/// Linear Equations, Matrix Inversion
2130/// trsequ solves the matrix equation
2131///
2132/// SMX*x = B
2133///
2134/// which represents a system of m simultaneous linear equations with n right-hand sides:
2135/// SMX is an unpacked symmetric matrix (all elements) (m x m)
2136/// B is an unpacked matrix of right-hand sides (n x m)
2137///
2138
2139double *TCL::trsequ(double *smx, int m, double *b, int n)
2140{
2141 double *mem = new double[(m*(m+1))/2+m];
2142 double *v = mem;
2143 double *s = v+m;
2144 if (!b) n=0;
2145 TCL::trpck (smx,s ,m);
2146 TCL::trsinv(s ,s, m);
2147
2148 for (int i=0;i<n;i++) {
2149 TCL::trsa (s ,b+i*m, v, m, 1);
2150 TCL::ucopy (v ,b+i*m, m);}
2151 TCL::trupck(s ,smx, m);
2152 delete [] mem;
2153 return b;
2154}
SVector< double, 2 > v
Definition: Dict.h:5
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
#define g(i)
Definition: RSha256.hxx:105
#define ClassImp(name)
Definition: Rtypes.h:363
#define TCL_TRALT(a, u, b, m, n)
Definition: TCernLib.cxx:376
#define TCL_TRASAT(a, s, r__, m, n)
Definition: TCernLib.cxx:496
#define TCL_MXMAD(n_, a, b, c, i, j, k)
Definition: TCernLib.cxx:65
#define TCL_MXMLRT(n__, a, b, c, ni, nj)
Definition: TCernLib.cxx:130
#define TCL_TRAL(a, u, b, m, n)
Definition: TCernLib.cxx:317
#define TCL_MXTRP(a, b, i, j)
Definition: TCernLib.cxx:213
#define TCL_TRAAT(a, s, m, n)
Definition: TCernLib.cxx:260
#define TCL_TRAS(a, s, b, m, n)
Definition: TCernLib.cxx:433
float * q
Definition: THbookFile.cxx:87
int * iq
Definition: THbookFile.cxx:86
Definition: TCernLib.h:36
static float * mxmad_0_(int n, const float *a, const float *b, float *c, int i, int j, int k)
Definition: TCernLib.cxx:109
static float * trpck(const float *s, float *u, int n)
trpck.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:1010
static float * trqsq(const float *q, const float *s, float *r, int m)
trqsq.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:1045
static float * tralt(const float *a, const float *u, float *b, int m, int n)
Triangular - Rectangular Multiplication.
Definition: TCernLib.cxx:407
static float * mxtrp(const float *a, float *b, int i, int j)
Matrix Transposition.
Definition: TCernLib.cxx:230
static float * trata(const float *a, float *r, int m, int n)
trata.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:589
static float * trsmul(const float *g, float *gi, int n)
trsmul.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:1212
static float * tratsa(const float *a, const float *s, float *r, int m, int n)
tratsa.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:677
static int * ucopy(const int *a, int *b, int n)
Definition: TCernLib.h:325
static float * trsa(const float *s, const float *a, float *b, int m, int n)
trsa.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:1111
static float * vzero(float *a, int n2)
Definition: TCernLib.h:485
static float * trchlu(const float *a, float *b, int n)
trchlu.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:733
static float * trsinv(const float *g, float *gi, int n)
trsinv.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:1160
static float * trsmlu(const float *u, float *s, int n)
trsmlu.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:1177
static float * trats(const float *a, const float *s, float *b, int m, int n)
trats.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:629
static float * trupck(const float *u, float *s, int m)
trupck.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:1248
static float * trasat(const float *a, const float *s, float *r, int m, int n)
Transformation of Symmetric Matrix.
Definition: TCernLib.cxx:540
static float * trlta(const float *u, const float *a, float *b, int m, int n)
trlta.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:959
static float * trla(const float *u, const float *a, float *b, int m, int n)
trla.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:915
static float * trinv(const float *t, float *s, int n)
trinv.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:857
static float * tras(const float *a, const float *s, float *b, int m, int n)
Symmetric - Rectangular Multiplication.
Definition: TCernLib.cxx:470
static float * tral(const float *a, const float *u, float *b, int m, int n)
Triangular - Rectangular Multiplication.
Definition: TCernLib.cxx:350
static float * mxmlrt_0_(int n__, const float *a, const float *b, float *c, int ni, int nj)
Matrix Multiplication.
Definition: TCernLib.cxx:184
static float * trsequ(float *smx, int m=3, float *b=0, int n=1)
Linear Equations, Matrix Inversion trsequ solves the matrix equation.
Definition: TCernLib.cxx:2111
static float * traat(const float *a, float *s, int m, int n)
Symmetric Multiplication of Rectangular Matrices.
Definition: TCernLib.cxx:293
static float * trsat(const float *s, const float *a, float *b, int m, int n)
trsat.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:1298
static float * trchul(const float *a, float *b, int n)
trchul.F – translated by f2c (version 19970219).
Definition: TCernLib.cxx:793
const Int_t n
Definition: legend1.C:16
static constexpr double s
static constexpr double m2
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12
static long int sum(long int i)
Definition: Factory.cxx:2258