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