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)
117 double *
TCL::mxmad_0_(
int n_,
const double *
a,
const double *
b,
double *c,
int i,
int j,
int 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]; \ 199 double *
TCL::mxmlrt_0_(
int n__,
const double *
a,
const double *
b,
double *c,
int ni,
int 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]; } 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); \ 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]; \ 366 double *
TCL::tral(
const double *
a,
const double *u,
double *
b,
int m,
int 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]; \ 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]; \ 486 double *
TCL::tras(
const double *
a,
const double *s,
double *
b,
int m,
int 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); \ 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];
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];
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__;
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__;
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);
917 int ipiv, ia, ib, iu;
925 ipiv = (m * m +
m) / 2;
934 sum += a[ia] * u[iu];
941 }
while (ia > 1 - 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);
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;
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;
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];
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];
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__;
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__;
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__;
1679 int ipiv, ia, ib, iu;
1687 ipiv = (m * m +
m) / 2;
1696 sum += a[ia] * u[iu];
1703 }
while (ia > 1 - 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);
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;
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;
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++) {
2141 double *mem =
new double[(m*(m+1))/2+
m];
2148 for (
int i=0;i<
n;i++) {
static float * traat(const float *a, float *s, int m, int n)
Symmetric Multiplication of Rectangular Matrices.
static long int sum(long int i)
static float * trsmul(const float *g, float *gi, int n)
trsmul.F – translated by f2c (version 19970219).
#define TCL_TRAAT(a, s, m, n)
#define TCL_TRASAT(a, s, r__, m, n)
static float * trata(const float *a, float *r, int m, int n)
trata.F – translated by f2c (version 19970219).
#define TCL_MXTRP(a, b, i, j)
static float * mxtrp(const float *a, float *b, int i, int j)
Matrix Transposition.
static float * trinv(const float *t, float *s, int n)
trinv.F – translated by f2c (version 19970219).
static float * trsat(const float *s, const float *a, float *b, int m, int n)
trsat.F – translated by f2c (version 19970219).
static float * mxmad_0_(int n, const float *a, const float *b, float *c, int i, int j, int k)
static float * trpck(const float *s, float *u, int n)
trpck.F – translated by f2c (version 19970219).
static float * trsequ(float *smx, int m=3, float *b=0, int n=1)
Linear Equations, Matrix Inversion trsequ solves the matrix equation.
static float * trsa(const float *s, const float *a, float *b, int m, int n)
trsa.F – translated by f2c (version 19970219).
static float * tralt(const float *a, const float *u, float *b, int m, int n)
Triangular - Rectangular Multiplication.
#define TCL_TRAS(a, s, b, m, n)
static float * trlta(const float *u, const float *a, float *b, int m, int n)
trlta.F – translated by f2c (version 19970219).
static float * trqsq(const float *q, const float *s, float *r, int m)
trqsq.F – translated by f2c (version 19970219).
static float * trats(const float *a, const float *s, float *b, int m, int n)
trats.F – translated by f2c (version 19970219).
static int * ucopy(const int *a, int *b, int n)
static float * tral(const float *a, const float *u, float *b, int m, int n)
Triangular - Rectangular Multiplication.
static float * trchul(const float *a, float *b, int n)
trchul.F – translated by f2c (version 19970219).
static float * trsinv(const float *g, float *gi, int n)
trsinv.F – translated by f2c (version 19970219).
#define TCL_TRALT(a, u, b, m, n)
static float * trasat(const float *a, const float *s, float *r, int m, int n)
Transformation of Symmetric Matrix.
static float * trsmlu(const float *u, float *s, int n)
trsmlu.F – translated by f2c (version 19970219).
static float * tras(const float *a, const float *s, float *b, int m, int n)
Symmetric - Rectangular Multiplication.
static float * mxmlrt_0_(int n__, const float *a, const float *b, float *c, int ni, int nj)
Matrix Multiplication.
static float * vzero(float *a, int n2)
#define TCL_TRAL(a, u, b, m, n)
#define TCL_MXMLRT(n__, a, b, c, ni, nj)
static float * trupck(const float *u, float *s, int m)
trupck.F – translated by f2c (version 19970219).
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
#define TCL_MXMAD(n_, a, b, c, i, j, k)
static float * tratsa(const float *a, const float *s, float *r, int m, int n)
tratsa.F – translated by f2c (version 19970219).
Double_t Sqrt(Double_t x)
static float * trchlu(const float *a, float *b, int n)
trchlu.F – translated by f2c (version 19970219).
static float * trla(const float *u, const float *a, float *b, int m, int n)
trla.F – translated by f2c (version 19970219).