StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TCernLib.cxx
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 
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 //
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 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 
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 
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 
183 
184 float *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 
198 
199 double *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 
229 
230 float *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 
244 
245 double *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 
292 
293 float *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 
308 
309 double *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 
349 
350 float *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 
365 
366 double *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 
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 
406 
407 float *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 
422 
423 double *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 
469 
470 float *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 
485 
486 double *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 
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 
539 
540 float *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 
555 
556 double *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 
571 
572 float *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 
588 
589 float *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 
628 
629 float *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 
676 
677 float *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 
732 
733 float *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 
767 L40:
768  sum = a[kpiv] - sum;
769 L42:
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 
792 
793 float *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 
829 L40:
830  sum = a[kpiv] - sum;
831 L42:
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 
856 
857 float *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 
914 
915 float *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 
958 
959 float *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 
1009 
1010 float *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 
1044 
1045 float *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 
1110 
1111 float *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 
1159 
1160 float *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 
1176 
1177 float *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 
1211 
1212 float *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 
1247 
1248 float *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 
1297 
1298 float *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 
1349 
1350 double *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 
1392 
1393 double *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 
1440 
1441 double *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 
1496 
1497 double *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 
1531 L40:
1532  sum = a[kpiv] - sum;
1533 L42:
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 
1556 
1557 double *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 
1593 L40:
1594  sum = a[kpiv] - sum;
1595 L42:
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 
1619 
1620 double *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 
1676 
1677 double *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 
1719 
1720 double *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 
1769 
1770 double *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 
1803 
1804 double *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 
1868 
1869 double *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 
1916 
1917 double *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 
1936 
1937 double *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 
1972 
1973 double *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 
2009 
2010 double *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 
2058 
2059 double *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 
2110 
2111 float *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 
2138 
2139 double *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 }
static float * traat(const float *a, float *s, int m, int n)
Definition: TCernLib.cxx:293
Definition: TCernLib.h:36
static float * trsmul(const float *g, float *gi, int n)
Definition: TCernLib.cxx:1212
static float * trata(const float *a, float *r, int m, int n)
Definition: TCernLib.cxx:589
static float * mxtrp(const float *a, float *b, int i, int j)
Definition: TCernLib.cxx:230
static float * trinv(const float *t, float *s, int n)
Definition: TCernLib.cxx:857
static float * trsat(const float *s, const float *a, float *b, int m, int n)
Definition: TCernLib.cxx:1298
static float * trpck(const float *s, float *u, int n)
Definition: TCernLib.cxx:1010
static float * trsequ(float *smx, int m=3, float *b=0, int n=1)
Definition: TCernLib.cxx:2111
static float * trsa(const float *s, const float *a, float *b, int m, int n)
Definition: TCernLib.cxx:1111
static float * tralt(const float *a, const float *u, float *b, int m, int n)
Definition: TCernLib.cxx:407
static float * trlta(const float *u, const float *a, float *b, int m, int n)
Definition: TCernLib.cxx:959
static float * trqsq(const float *q, const float *s, float *r, int m)
Definition: TCernLib.cxx:1045
static float * trats(const float *a, const float *s, float *b, int m, int n)
Definition: TCernLib.cxx:629
static float * tral(const float *a, const float *u, float *b, int m, int n)
Definition: TCernLib.cxx:350
static float * trchul(const float *a, float *b, int n)
Definition: TCernLib.cxx:793
static float * trsinv(const float *g, float *gi, int n)
Definition: TCernLib.cxx:1160
static float * trasat(const float *a, const float *s, float *r, int m, int n)
Definition: TCernLib.cxx:540
static float * trsmlu(const float *u, float *s, int n)
Definition: TCernLib.cxx:1177
static float * tras(const float *a, const float *s, float *b, int m, int n)
Definition: TCernLib.cxx:470
static float * mxmlrt_0_(int n__, const float *a, const float *b, float *c, int ni, int nj)
Definition: TCernLib.cxx:184
static float * trupck(const float *u, float *s, int m)
Definition: TCernLib.cxx:1248
static float * tratsa(const float *a, const float *s, float *r, int m, int n)
Definition: TCernLib.cxx:677
static float * trchlu(const float *a, float *b, int n)
Definition: TCernLib.cxx:733
static float * trla(const float *u, const float *a, float *b, int m, int n)
Definition: TCernLib.cxx:915