65 #define TCL_MXMAD(n_,a,b,c,i,j,k) \
67 int l, m, n, ia, ic, ib, ja, jb, iia, iib, ioa, iob; \
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; \
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); \
89 for (l = 1; l <= i; ++l) { \
91 for (m = 1; m <= k; ++m,++ic) { \
93 case 1: c[ic] = 0.; break; \
94 case 3: c[ic] = -c[ic]; break; \
96 if (j == 0) continue; \
99 for (n = 1; n <= j; ++n, ja+=iia, jb+=iib) \
100 cic += a[ja] * b[jb]; \
109 float *TCL::mxmad_0_(
int n_,
const float *a,
const float *b,
float *c,
int i,
int j,
int k)
111 TCL_MXMAD(n_,a,b,c,i,j,k)
117 double *
TCL::mxmad_0_(
int n_, const
double *a, const
double *b,
double *c,
int i,
int j,
int k)
119 TCL_MXMAD(n_,a,b,c,i,j,k)
130 #define TCL_MXMLRT( n__, a, b, c, ni,nj) \
131 if (ni <= 0 || nj <= 0) return 0; \
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; } \
140 for (ii = 1; ii <= ni; ++ii, ic1+=ni, ia1+=jpa) { \
142 for (kc = 1; kc <= ni; ++kc,ic++) c[ic] = 0.; \
144 for (jj = 1; jj <= nj; ++jj,++ib1,ja1 += ipa) { \
145 ib = ib1; ia = ia1; \
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]; \
184 float *
TCL::mxmlrt_0_(
int n__,
const float *a,
const float *b,
float *c,
int ni,
int nj)
186 TCL_MXMLRT( n__, a, b, c, ni,nj)
199 double *
TCL::mxmlrt_0_(
int n__,
const double *a,
const double *b,
double *c,
int ni,
int nj)
201 TCL_MXMLRT( n__, a, b, c, ni,nj)
213 #define TCL_MXTRP(a, b, i, j) \
214 if (i == 0 || j == 0) return 0; \
217 for (int k = 1; k <= j; ++k) \
219 for (int l = 1; l <= i; ++l,ia += j,++ib) b[ib] = a[ia]; }
232 TCL_MXTRP(a, b, i, j)
247 TCL_MXTRP(a, b, i, j)
260 #define TCL_TRAAT(a, s, m, n) \
262 int ipiv, i, j, ipivn, ia, is, iat; \
266 for (i = 1; i <= m; ++i) { \
270 for (j = 1; j <= i; ++j) { \
275 sum += a[ia] * a[iat]; \
276 } while (ia < ipivn); \
295 TCL_TRAAT(a, s, m, n)
311 TCL_TRAAT(a, s, m, n)
317 #define TCL_TRAL(a, u, b, m, n) \
318 int indu, i, j, k, ia, ib, iu; \
322 for (i = 1; i <= m; ++i) { \
324 for (j = 1; j <= n; ++j) { \
329 for (k = j; k <= n; ++k) {\
330 sum += a[ia] * u[iu]; \
350 float *
TCL::tral(
const float *a,
const float *u,
float *b,
int m,
int n)
352 TCL_TRAL(a, u, b, m, n)
366 double *
TCL::tral(
const double *a,
const double *u,
double *b,
int m,
int n)
368 TCL_TRAL(a, u, b, m, n)
376 #define TCL_TRALT(a, u, b, m, n) \
377 int indu, j, k, ia, ib, iu; \
381 indu = (n * n + n) / 2; \
384 for (j = 1; j <= n; ++j) { \
387 for (k = j; k <= n; ++k) {\
388 sum += a[ia] * u[iu]; \
407 float *
TCL::tralt(
const float *a,
const float *u,
float *b,
int m,
int n)
409 TCL_TRALT(a, u, b, m, n)
423 double *
TCL::tralt(
const double *a,
const double *u,
double *b,
int m,
int n)
425 TCL_TRALT(a, u, b, m, n)
433 #define TCL_TRAS(a, s, b, m, n) \
434 int inds, i__, j, k, ia, ib, is; \
437 ib = 0; inds = 0; i__ = 0; \
442 for (j = 1; j <= m; ++j) { \
447 if (k > i__) is += k; \
450 sum += a[ia] * s[is]; \
470 float *
TCL::tras(
const float *a,
const float *s,
float *b,
int m,
int n)
472 TCL_TRAS(a, s, b, m, n)
486 double *
TCL::tras(
const double *a,
const double *s,
double *b,
int m,
int n)
488 TCL_TRAS(a, s, b, m, n)
496 #define TCL_TRASAT(a, s, r__, m, n) \
498 int ia, mn, ir, is, iaa; \
501 imax = (m * m + m) / 2; \
502 vzero(&r__[1], imax); \
513 if (k > i__) is += k; \
516 sum += s[is] * a[ia]; \
522 r__[ir] += sum * a[iaa];\
524 } while (iaa <= ia); \
540 float *
TCL::trasat(
const float *a,
const float *s,
float *r__,
int m,
int n)
542 TCL_TRASAT(a, s, r__, m, n)
556 double *
TCL::trasat(
const double *a,
const double *s,
double *r__,
int m,
int n)
558 TCL_TRASAT(a, s, r__, m, n)
572 float *
TCL::trasat(
const double *a,
const float *s,
float *r__,
int m,
int n)
574 TCL_TRASAT(a, s, r__, m, n)
592 int i__, j, ia, mn, ir, iat;
602 for (i__ = 1; i__ <= m; ++i__) {
603 for (j = 1; j <= i__; ++j) {
608 sum += a[ia] * a[iat];
629 float *
TCL::trats(
const float *a,
const float *s,
float *b,
int m,
int n)
632 int inds, i__, j, k, ia, ib, is;
639 ib = 0; inds = 0; i__ = 0;
644 for (j = 1; j <= m; ++j) {
651 if (k > i__) is += k;
653 sum += a[ia] * s[is];
677 float *
TCL::tratsa(
const float *a,
const float *s,
float *r__,
int m,
int n)
681 int ia, ir, is, iaa, ind;
688 imax = (m * m + m) / 2;
689 vzero(&r__[1], imax);
697 for (j = 1; j <= m; ++j) {
704 if (k > i__) is += k;
706 sum += s[is] * a[ia];
712 for (k = 1; k <= j; ++k) {
715 r__[ir] += sum * a[iaa];
736 int ipiv, kpiv, i__, j;
755 for (j = i__; j <= n; ++j) {
757 if (i__ == 1)
goto L40;
758 if (r__ == 0.)
goto L42;
763 sum += b[kd] * b[id];
770 if (j != i__) b[kpiv] = sum * r__;
772 dc = TMath::Sqrt(sum);
774 if (r__ > 0.) r__ = 1. / dc;
807 kpiv = (n * n + n) / 2;
816 if (i__ == n)
goto L40;
817 if (r__ == 0.)
goto L42;
826 sum += b[id] * b[kd];
832 if (kpiv < ipiv) b[kpiv] = sum * r__;
834 dc = TMath::Sqrt(sum);
836 if (r__ > 0.) r__ = 1. / dc;
839 }
while (kpiv > ipiv - i__);
859 int lhor, ipiv, lver, j;
868 mx = (n * n + n) / 2;
874 if (t[ipiv] > 0.) r__ = 1. / t[ipiv];
879 while (ind != ipiv) {
889 sum += t[lhor] * s[lver];
891 }
while (lhor < ind);
915 float *
TCL::trla(
const float *u,
const float *a,
float *b,
int m,
int n)
917 int ipiv, ia, ib, iu;
925 ipiv = (m * m + m) / 2;
934 sum += a[ia] * u[iu];
941 }
while (ia > 1 - n);
959 float *
TCL::trlta(
const float *u,
const float *a,
float *b,
int m,
int n)
961 int ipiv, mxpn, i__, nTep, ia, ib, iu, mx;
986 sum += a[ia] * u[iu];
1012 int i__, ia, ind, ipiv;
1022 for (i__ = 1; i__ <= n; ++i__) {
1028 }
while (ind < ipiv);
1045 float *
TCL::trqsq(
const float *q,
const float *s,
float *r__,
int m)
1047 int indq, inds, imax, i__, j, k, l;
1048 int iq, ir, is, iqq;
1055 imax = (m * m + m) / 2;
1056 vzero(&r__[1], imax);
1074 if (k > i__) is += k;
1080 sum += s[is] * q[iq];
1088 if (l > i__) iqq += l;
1090 r__[ir] += q[iqq] * sum;
1111 float *
TCL::trsa(
const float *s,
const float *a,
float *b,
int m,
int n)
1114 int inds, i__, j, k, ia, ib, is;
1128 for (j = 1; j <= n; ++j) {
1135 if (k > i__) is += k;
1137 sum += s[is] * a[ia];
1165 return trsmul(gi, gi, n);
1180 int lhor, lver, i__, k, l, ind;
1187 ind = (n * n + n) / 2;
1189 for (i__ = 1; i__ <= n; ++i__) {
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];
1215 int lhor, lver, lpiv, i__, j, k, ind;
1224 for (i__ = 1; i__ <= n; ++i__) {
1226 for (j = 1; j <= i__; ++j,++ind) {
1230 for (k = i__; k <= n; lhor += k,lver += k,++k)
1231 sum += g[lver] * g[lhor];
1250 int i__, im, is, iu, iv, ih, m2;
1298 float *
TCL::trsat(
const float *s,
const float *a,
float *b,
int m,
int n)
1302 int inds, i__, j, k, ia, ib, is;
1317 for (j = 1; j <= n; ++j) {
1323 if (k > i__) is += k;
1326 sum += s[is] * a[ia];
1354 int i__, j, ia, mn, ir, iat;
1364 for (i__ = 1; i__ <= m; ++i__) {
1366 for (j = 1; j <= i__; ++j) {
1372 sum += a[ia] * a[iat];
1393 double *
TCL::trats(
const double *a,
const double *s,
double *b,
int m,
int n)
1396 int inds, i__, j, k, ia, ib, is;
1403 ib = 0; inds = 0; i__ = 0;
1409 for (j = 1; j <= m; ++j) {
1416 if (k > i__) is += k;
1418 sum += a[ia] * s[is];
1441 double *
TCL::tratsa(
const double *a,
const double *s,
double *r__,
int m,
int n)
1444 int imax, i__, j, k;
1445 int ia, ir, is, iaa, ind;
1452 imax = (m * m + m) / 2;
1453 vzero(&r__[1], imax);
1461 for (j = 1; j <= m; ++j) {
1468 if (k > i__) is += k;
1470 sum += s[is] * a[ia];
1476 for (k = 1; k <= j; ++k) {
1479 r__[ir] += sum * a[iaa];
1500 int ipiv, kpiv, i__, j;
1519 for (j = i__; j <= n; ++j) {
1521 if (i__ == 1)
goto L40;
1522 if (r__ == 0.)
goto L42;
1523 id = ipiv - i__ + 1;
1524 kd = kpiv - i__ + 1;
1527 sum += b[kd] * b[id];
1529 }
while (
id < ipiv);
1532 sum = a[kpiv] - sum;
1534 if (j != i__) b[kpiv] = sum * r__;
1536 dc = TMath::Sqrt(sum);
1538 if (r__ > 0.) r__ = (double)1. / dc;
1560 int ipiv, kpiv, i__;
1571 kpiv = (n * n + n) / 2;
1580 if (i__ == n)
goto L40;
1581 if (r__ == (
double)0.)
goto L42;
1590 sum += b[id] * b[kd];
1594 sum = a[kpiv] - sum;
1596 if (kpiv < ipiv) b[kpiv] = sum * r__;
1598 dc = TMath::Sqrt(sum);
1600 if (r__ > (
double)0.) r__ = (double)1. / dc;
1603 }
while (kpiv > ipiv - i__);
1622 int lhor, ipiv, lver, j;
1631 mx = (n * n + n) / 2;
1637 if (t[ipiv] > 0.) r__ = (double)1. / t[ipiv];
1642 while (ind != ipiv) {
1652 sum += t[lhor] * s[lver];
1654 }
while (lhor < ind);
1656 s[ind] = -sum * r__;
1677 double *
TCL::trla(
const double *u,
const double *a,
double *b,
int m,
int n)
1679 int ipiv, ia, ib, iu;
1687 ipiv = (m * m + m) / 2;
1696 sum += a[ia] * u[iu];
1703 }
while (ia > 1 - n);
1720 double *
TCL::trlta(
const double *u,
const double *a,
double *b,
int m,
int n)
1722 int ipiv, mxpn, i__, nTep, ia, ib, iu, mx;
1747 sum += a[ia] * u[iu];
1754 }
while (ia < mxpn);
1772 int i__, ia, ind, ipiv;
1782 for (i__ = 1; i__ <= n; ++i__) {
1788 }
while (ind < ipiv);
1804 double *
TCL::trqsq(
const double *q,
const double *s,
double *r__,
int m)
1806 int indq, inds, imax, i__, j, k, l;
1807 int iq, ir, is, iqq;
1814 imax = (m * m + m) / 2;
1815 vzero(&r__[1], imax);
1833 if (k > i__) is += k;
1839 sum += s[is] * q[iq];
1847 if (l > i__) iqq += l;
1849 r__[ir] += q[iqq] * sum;
1869 double *
TCL::trsa(
const double *s,
const double *a,
double *b,
int m,
int n)
1872 int inds, i__, j, k, ia, ib, is;
1886 for (j = 1; j <= n; ++j) {
1893 if (k > i__) is += k;
1895 sum += s[is] * a[ia];
1941 int lhor, lver, i__, k, l, ind;
1948 ind = (n * n + n) / 2;
1950 for (i__ = 1; i__ <= n; ++i__) {
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];
1977 int lhor, lver, lpiv, i__, j, k, ind;
1986 for (i__ = 1; i__ <= n; ++i__) {
1988 for (j = 1; j <= i__; ++j,++ind) {
1992 for (k = i__; k <= n;lhor += k,lver += k,++k)
1993 sum += g[lver] * g[lhor];
2012 int i__, im, is, iu, iv, ih, m2;
2059 double *
TCL::trsat(
const double *s,
const double *a,
double *b,
int m,
int n)
2062 int inds, i__, j, k, ia, ib, is;
2077 for (j = 1; j <= n; ++j) {
2083 if (k > i__) is += k;
2086 sum += s[is] * a[ia];
2113 float *mem =
new float[(m*(m+1))/2+m];
2120 for (
int i=0;i<n;i++) {
2122 TCL::ucopy (v ,b+i*m, m);}
2141 double *mem =
new double[(m*(m+1))/2+m];
2148 for (
int i=0;i<n;i++) {
2150 TCL::ucopy (v ,b+i*m, m);}
static float * traat(const float *a, float *s, int m, int n)
static float * trsmul(const float *g, float *gi, int n)
static float * trata(const float *a, float *r, int m, int n)
static float * mxtrp(const float *a, float *b, int i, int j)
static float * trinv(const float *t, float *s, int n)
static float * trsat(const float *s, const float *a, float *b, int m, int n)
static float * trpck(const float *s, float *u, int n)
static float * trsequ(float *smx, int m=3, float *b=0, int n=1)
static float * trsa(const float *s, const float *a, float *b, int m, int n)
static float * tralt(const float *a, const float *u, float *b, int m, int n)
static float * trlta(const float *u, const float *a, float *b, int m, int n)
static float * trqsq(const float *q, const float *s, float *r, int m)
static float * trats(const float *a, const float *s, float *b, int m, int n)
static float * tral(const float *a, const float *u, float *b, int m, int n)
static float * trchul(const float *a, float *b, int n)
static float * trsinv(const float *g, float *gi, int n)
static float * trasat(const float *a, const float *s, float *r, int m, int n)
static float * trsmlu(const float *u, float *s, int n)
static float * tras(const float *a, const float *s, float *b, int m, int n)
static float * mxmlrt_0_(int n__, const float *a, const float *b, float *c, int ni, int nj)
static float * trupck(const float *u, float *s, int m)
static float * tratsa(const float *a, const float *s, float *r, int m, int n)
static float * trchlu(const float *a, float *b, int n)
static float * trla(const float *u, const float *a, float *b, int m, int n)