// @(#)root/star:$Name: $:$Id: TCL.cxx,v 1.1 2002/05/27 16:26:59 rdm Exp $
// Author: Valery Fine(fine@bnl.gov) 25/09/99
//
// The set of methods to work with the plain matrix / vector
// "derived" from http://wwwinfo.cern.ch/asdoc/shortwrupsdir/f110/top.html
// "derived" from http://wwwinfo.cern.ch/asdoc/shortwrupsdir/f112/top.html
//
// $Id: TCL.cxx,v 1.1 2002/05/27 16:26:59 rdm Exp $
// $Log: TCL.cxx,v $
// Revision 1.1 2002/05/27 16:26:59 rdm
// rename star to table.
//
// Revision 1.2 2001/05/29 19:08:08 brun
// New version of some STAR classes from Valery.
//
// Revision 1.2 2001/05/27 02:38:14 fine
// New method trsedu to solev Ax=B from Victor
//
// Revision 1.1.1.1 2000/11/27 22:57:14 fisyak
//
//
// Revision 1.1.1.1 2000/05/16 17:00:48 rdm
// Initial import of ROOT into CVS
//
////////////////////////////////////////////////////////////////////////////////////////////////////////
#include <assert.h>
#include "TCL.h"
#include "TMath.h"
#include <TArrayD.h>
#include <TError.h>
ClassImp(TCL)
#define TCL_MXMAD(n_,a,b,c,i,j,k) \
/* Local variables */ \
int l, m, n, ia, ic, ib, ja, jb, iia, iib, ioa, iob; \
\
/* Parameter adjuTments */ \
--a; --b; --c; \
/* Function Body */ \
/* MXMAD MXMAD1 MXMAD2 MXMAD3 MXMPY MXMPY1 MXMPY2 MXMPY3 MXMUB MXMUB1 MXMUB2 MXMUB3 */ \
/* const int iandj1[] = {21, 22, 23, 24, 11, 12, 13, 14, 31, 32, 33, 34 }; */ \
const int iandj1[] = {2, 2 , 2 , 2 , 1 , 1 , 1 , 1 , 3 , 3 , 3 , 3 }; \
const int iandj2[] = { 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4 }; \
int n1 = iandj1[n_]; \
int n2 = iandj2[n_]; \
if (i == 0 || k == 0) return 0; \
\
switch (n2) { \
case 1: iia = 1; ioa = j; iib = k; iob = 1; break; \
case 2: iia = 1; ioa = j; iib = 1; iob = j; break; \
case 3: iia = i; ioa = 1; iib = k; iob = 1; break; \
case 4: iia = i; ioa = 1; iib = 1; iob = j; break; \
default: iia = ioa = iib = iob = 0; assert(iob); \
}; \
\
ia = 1; ic = 1; \
for (l = 1; l <= i; ++l) { \
ib = 1; \
for (m = 1; m <= k; ++m,++ic) { \
switch (n1) { \
case 1: c[ic] = 0.; break; \
case 3: c[ic] = -c[ic]; break; \
}; \
if (j == 0) continue; \
ja = ia; jb = ib; \
double cic = c[ic]; \
for (n = 1; n <= j; ++n, ja+=iia, jb+=iib) \
cic += a[ja] * b[jb]; \
c[ic] = cic; \
ib += iob; \
} \
ia += ioa; \
}
//___________________________________________________________________________
float *TCL::mxmad_0_(int n_, const float *a, const float *b, float *c, int i, int j, int k)
{
TCL_MXMAD(n_,a,b,c,i,j,k)
return c;
} /* mxmad_ */
//___________________________________________________________________________
double *TCL::mxmad_0_(int n_, const double *a, const double *b, double *c, int i, int j, int k)
{
TCL_MXMAD(n_,a,b,c,i,j,k)
return c;
} /* mxmad_ */
#undef TCL_MXMAD
//___________________________________________________________________________
//
// Matrix Multiplication
//___________________________________________________________________________
#define TCL_MXMLRT( n__, a, b, c, ni,nj) \
if (ni <= 0 || nj <= 0) return 0; \
double x; \
int ia, ib, ic, ja, kc, ii, jj, kj, ki, ia1, ib1, ic1, ja1; \
int ipa = 1; int jpa = nj; \
if (n__ == 1) { ipa = ni; jpa = 1; } \
\
--a; --b; --c; \
\
ic1 = 1; ia1 = 1; \
for (ii = 1; ii <= ni; ++ii, ic1+=ni, ia1+=jpa) { \
ic = ic1; \
for (kc = 1; kc <= ni; ++kc,ic++) c[ic] = 0.; \
ib1 = 1; ja1 = 1; \
for (jj = 1; jj <= nj; ++jj,++ib1,ja1 += ipa) { \
ib = ib1; ia = ia1; \
x = 0.; \
for (kj = 1;kj <= nj;++kj,ia+=ipa,ib += nj) \
x += a[ia] * b[ib]; \
ja = ja1; ic = ic1; \
for (ki = 1; ki <= ni; ++ki,++ic,ja += jpa) \
c[ic] += x * a[ja]; \
} \
}
//___________________________________________________________________________
float *TCL::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. 01/01/64 RKB
// -- ENTRY MXMLRT */
// -- C = A(I,J) X B(J,J) X A*(J,I) */
// -- A* TANDS FOR A-TRANSPOSED */
//
//
// mxmlrt (A,B,C,NI,NJ) IS EQUIVALENT TO */
// CALL MXMPY (A,B,X,NI,NJ,NJ) */
// CALL MXMPY1 (X,A,C,NI,NJ,NI) */
/* OR CALL MXMPY1 (B,A,Y,NJ,NJ,NI) */
/* CALL MXMPY (A,Y,C,NI,NJ,NI) */
// -- C = A*(I,J) X B(J,J) X A(J,I)
//
//
// CALL MXMLTR (A,B,C,NI,NJ) IS EQUIVALENT TO
// CALL MXMPY2 (A,B,X,NI,NJ,NJ)
// CALL MXMPY (X,A,C,NI,NJ,NI)
// OR CALL MXMPY (B,A,Y,NJ,NJ,NI)
// CALL MXMPY2 (A,Y,C,NI,NJ,NI)
TCL_MXMLRT( n__, a, b, c, ni,nj)
return c;
} /* mxmlrt_ */
//___________________________________________________________________________
double *TCL::mxmlrt_0_(int n__, const double *a, const double *b, double *c, int ni,int nj)
{
// Matrix Multiplication (double precision)
//
//
//
//
TCL_MXMLRT( n__, a, b, c, ni,nj)
return c;
} /* mxmlrt_ */
#undef TCL_MXMLRT
//___________________________________________________________________________
//
// Matrix Transposition
//___________________________________________________________________________
#define TCL_MXTRP(a, b, i, j) \
if (i == 0 || j == 0) return 0; \
--b; --a; \
int ib = 1; \
for (int k = 1; k <= j; ++k) \
{ int ia = k; \
for (int l = 1; l <= i; ++l,ia += j,++ib) b[ib] = a[ia]; }
//___________________________________________________________________________
float *TCL::mxtrp(const float *a, float *b, int i, int j)
{
//
// Matrix Transposition
//
//
// CERN PROGLIB# F110 MXTRP .VERSION KERNFOR 1.0 650809
// ORIG. 01/01/64 RKB
TCL_MXTRP(a, b, i, j)
return b;
} /* mxtrp */
//___________________________________________________________________________
double *TCL::mxtrp(const double *a, double *b, int i, int j)
{
// Matrix Transposition (double precision)
//
//
// CERN PROGLIB# F110 MXTRP .VERSION KERNFOR 1.0 650809
// ORIG. 01/01/64 RKB
TCL_MXTRP(a, b, i, j)
return b;
} /* mxtrp */
#undef TCL_MXTRP
//___________________________________________________________________________
//___________________________________________________________________________
//
// TRPACK
//___________________________________________________________________________
//___________________________________________________________________________
#define TCL_TRAAT(a, s, m, n) \
/* Local variables */ \
int ipiv, i, j, ipivn, ia, is, iat; \
double sum; \
--s; --a; \
ia = 0; is = 0; \
for (i = 1; i <= m; ++i) { \
ipiv = ia; \
ipivn = ipiv + n; \
iat = 0; \
for (j = 1; j <= i; ++j) { \
ia = ipiv; \
sum = 0.; \
do { \
++ia; ++iat; \
sum += a[ia] * a[iat]; \
} while (ia < ipivn); \
++is; \
s[is] = sum; \
} \
} \
s++;
//____________________________________________________________
float *TCL::traat(float *a, float *s, int m, int n)
{
//
// Symmetric Multiplication of Rectangular Matrices
// CERN PROGLIB# F112 TRAAT .VERSION KERNFOR 4.15 861204
// ORIG. 18/12/74 WH */
// traat.F -- translated by f2c (version 19970219).
//
TCL_TRAAT(a, s, m, n)
return s;
} /* traat_ */
//____________________________________________________________
double *TCL::traat(double *a, double *s, int m, int n)
{
// Symmetric Multiplication of Rectangular Matrices
// CERN PROGLIB# F112 TRAAT .VERSION KERNFOR 4.15 861204
// ORIG. 18/12/74 WH */
// traat.F -- translated by f2c (version 19970219).
//
TCL_TRAAT(a, s, m, n)
return s;
} /* traat_ */
#undef TCL_TRAAT
#define TCL_TRAL(a, u, b, m, n) \
int indu, i, j, k, ia, ib, iu; \
double sum; \
--b; --u; --a; \
ib = 1; \
for (i = 1; i <= m; ++i) { \
indu = 0; \
for (j = 1; j <= n; ++j) { \
indu += j; \
ia = ib; \
iu = indu; \
sum = 0.; \
for (k = j; k <= n; ++k) { \
sum += a[ia] * u[iu]; \
++ia; \
iu += k; \
} \
b[ib] = sum; \
++ib; \
} \
} \
b++;
//____________________________________________________________
float *TCL::tral(float *a, float *u, float *b, int m, int n)
{
// Triangular - Rectangular Multiplication
// CERN PROGLIB# F112 TRAL .VERSION KERNFOR 4.15 861204
// ORIG. 18/12/74 WH
// tral.F -- translated by f2c (version 19970219).
TCL_TRAL(a, u, b, m, n)
return b;
} /* tral_ */
//____________________________________________________________
double *TCL::tral(double *a, double *u, double *b, int m, int n)
{
// Triangular - Rectangular Multiplication
// tral.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRAL .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
TCL_TRAL(a, u, b, m, n)
return b;
} /* tral_ */
#undef TCL_TRAL
//____________________________________________________________
#define TCL_TRALT(a, u, b, m, n) \
int indu, j, k, ia, ib, iu; \
double sum; \
--b; --u; --a; \
ib = m * n; \
indu = (n * n + n) / 2; \
do { \
iu = indu; \
for (j = 1; j <= n; ++j) { \
ia = ib; \
sum = 0.; \
for (k = j; k <= n; ++k) {\
sum += a[ia] * u[iu]; \
--ia; --iu; \
} \
b[ib] = sum; \
--ib; \
} \
} while (ib > 0); \
++b;
//____________________________________________________________
float *TCL::tralt(float *a, float *u, float *b, int m, int n)
{
// Triangular - Rectangular Multiplication
// CERN PROGLIB# F112 TRALT .VERSION KERNFOR 4.15 861204
// ORIG. 18/12/74 WH
// tralt.F -- translated by f2c (version 19970219).
TCL_TRALT(a, u, b, m, n)
return b;
} /* tralt_ */
//____________________________________________________________
double *TCL::tralt(double *a, double *u, double *b, int m, int n)
{
// Triangular - Rectangular Multiplication
// CERN PROGLIB# F112 TRALT .VERSION KERNFOR 4.15 861204
// ORIG. 18/12/74 WH
// tralt.F -- translated by f2c (version 19970219).
TCL_TRALT(a, u, b, m, n)
return b;
} /* tralt_ */
#undef TCL_TRALT
//____________________________________________________________
#define TCL_TRAS(a, s, b, m, n) \
int inds, i__, j, k, ia, ib, is; \
double sum; \
--b; --s; --a; \
ib = 0; inds = 0; i__ = 0; \
do { \
inds += i__; \
ia = 0; \
ib = i__ + 1; \
for (j = 1; j <= m; ++j) { \
is = inds; \
sum = 0.; \
k = 0; \
do { \
if (k > i__) is += k; \
else ++is; \
++ia; \
sum += a[ia] * s[is]; \
++k; \
} while (k < n); \
b[ib] = sum; \
ib += n; \
} \
++i__; \
} while (i__ < n); \
++b;
//____________________________________________________________
float *TCL::tras(float *a, float *s, float *b, int m, int n)
{
// Symmetric - Rectangular Multiplication
// CERN PROGLIB# F112 TRAS .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
// tras.F -- translated by f2c (version 19970219).
TCL_TRAS(a, s, b, m, n)
return b;
} /* tras_ */
//____________________________________________________________
double *TCL::tras(double *a, double *s, double *b, int m, int n)
{
// Symmetric - Rectangular Multiplication
// CERN PROGLIB# F112 TRAS .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
// tras.F -- translated by f2c (version 19970219).
TCL_TRAS(a, s, b, m, n)
return b;
} /* tras_ */
#undef TCL_TRAS
//____________________________________________________________
#define TCL_TRASAT(a, s, r__, m, n) \
int imax, k; \
int ia, mn, ir, is, iaa; \
double sum; \
--r__; --s; --a; \
imax = (m * m + m) / 2; \
vzero(&r__[1], imax); \
mn = m * n; \
int ind = 0; \
int i__ = 0; \
do { \
ind += i__; \
ia = 0; ir = 0; \
do { \
is = ind; \
sum = 0.; k = 0; \
do { \
if (k > i__) is += k; \
else ++is; \
++ia; \
sum += s[is] * a[ia]; \
++k; \
} while (k < n); \
iaa = i__ + 1; \
do { \
++ir; \
r__[ir] += sum * a[iaa]; \
iaa += n; \
} while (iaa <= ia); \
} while (ia < mn); \
++i__; \
} while (i__ < n); \
++r__;
//____________________________________________________________
float *TCL::trasat(float *a, float *s, float *r__, int m, int n)
{
// Transformation of Symmetric Matrix
// CERN PROGLIB# F112 TRASAT .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
// trasat.F -- translated by f2c (version 19970219).
TCL_TRASAT(a, s, r__, m, n)
return r__;
} /* trasat_ */
//____________________________________________________________
double *TCL::trasat(double *a, double *s, double *r__, int m, int n)
{
// Transformation of Symmetric Matrix
// CERN PROGLIB# F112 TRASAT .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
// trasat.F -- translated by f2c (version 19970219).
TCL_TRASAT(a, s, r__, m, n)
return r__;
} /* trasat_ */
//____________________________________________________________
float *TCL::trasat(double *a, float *s, float *r__, int m, int n)
{
// Transformation of Symmetric Matrix
// CERN PROGLIB# F112 TRASAT .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
// trasat.F -- translated by f2c (version 19970219).
TCL_TRASAT(a, s, r__, m, n)
return r__;
} /* trasat_ */
#undef TCL_TRASAT
//____________________________________________________________
float *TCL::trata(float *a, float *r__, int m, int n)
{
// trata.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRATA .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
/* Local variables */
int i__, j, ia, mn, ir, iat;
double sum;
/* Parameter adjuTments */
--r__; --a;
/* Function Body */
mn = m * n;
ir = 0;
for (i__ = 1; i__ <= m; ++i__) {
for (j = 1; j <= i__; ++j) {
ia = i__;
iat = j;
sum = 0.;
do {
sum += a[ia] * a[iat];
ia += m;
iat += m;
} while (ia <= mn);
++ir;
r__[ir] = sum;
}
}
++r__;
return r__;
} /* trata_ */
//____________________________________________________________
// trats.F -- translated by f2c (version 19970219).
float *TCL::trats(float *a, float *s, float *b, int m, int n)
{
/* Local variables */
int inds, i__, j, k, ia, ib, is;
double sum;
/* CERN PROGLIB# F112 TRATS .VERSION KERNFOR 4.15 861204 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--b; --s; --a;
/* Function Body */
ib = 0; inds = 0; i__ = 0;
do {
inds += i__;
ib = i__ + 1;
for (j = 1; j <= m; ++j) {
ia = j;
is = inds;
sum = 0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
sum += a[ia] * s[is];
ia += m;
++k;
} while (k < n);
b[ib] = sum;
ib += n;
}
++i__;
} while (i__ < n);
++b;
return b;
} /* trats_ */
//____________________________________________________________
// tratsa.F -- translated by f2c (version 19970219).
/* Subroutine */float *TCL::tratsa(float *a, float *s, float *r__, int m, int n)
{
/* Local variables */
int imax, i__, j, k;
int ia, mn, ir, is, iaa, ind;
double sum;
/* CERN PROGLIB# F112 TRATSA .VERSION KERNFOR 4.15 861204 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--r__; --s; --a;
/* Function Body */
imax = (m * m + m) / 2;
vzero(&r__[1], imax);
mn = m * n;
ind = 0;
i__ = 0;
do {
ind += i__;
ir = 0;
for (j = 1; j <= m; ++j) {
is = ind;
ia = j;
sum = 0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
sum += s[is] * a[ia];
ia += m;
++k;
} while (k < n);
iaa = i__ * m;
for (k = 1; k <= j; ++k) {
++iaa;
++ir;
r__[ir] += sum * a[iaa];
}
}
++i__;
} while (i__ < n);
++r__;
return r__;
} /* tratsa_ */
//____________________________________________________________
// trchlu.F -- translated by f2c (version 19970219).
float *TCL::trchlu(float *a, float *b, int n)
{
/* Local variables */
int ipiv, kpiv, i__, j;
double r__, dc;
int id, kd;
double sum;
/* CERN PROGLIB# F112 TRCHLU .VERSION KERNFOR 4.16 870601 */
/* ORIG. 18/12/74 W.HART */
/* Parameter adjuTments */
--b; --a;
/* Function Body */
ipiv = 0;
i__ = 0;
do {
++i__;
ipiv += i__;
kpiv = ipiv;
r__ = a[ipiv];
for (j = i__; j <= n; ++j) {
sum = 0.;
if (i__ == 1) goto L40;
if (r__ == 0.) goto L42;
id = ipiv - i__ + 1;
kd = kpiv - i__ + 1;
do {
sum += b[kd] * b[id];
++kd; ++id;
} while (id < ipiv);
L40:
sum = a[kpiv] - sum;
L42:
if (j != i__) b[kpiv] = sum * r__;
else {
dc = TMath::Sqrt(sum);
b[kpiv] = dc;
if (r__ > 0.) r__ = 1. / dc;
}
kpiv += j;
}
} while (i__ < n);
++b;
return b;
} /* trchlu_ */
//____________________________________________________________
// trchul.F -- translated by f2c (version 19970219).
/* Subroutine */float *TCL::trchul(float *a, float *b, int n)
{
/* Local variables */
int ipiv, kpiv, i__;
double r__;
int nTep;
double dc;
int id, kd;
double sum;
/* CERN PROGLIB# F112 TRCHUL .VERSION KERNFOR 4.16 870601 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--b; --a;
/* Function Body */
kpiv = (n * n + n) / 2;
i__ = n;
do {
ipiv = kpiv;
r__ = a[ipiv];
do {
sum = 0.;
if (i__ == n) goto L40;
if (r__ == 0.) goto L42;
id = ipiv;
kd = kpiv;
nTep = i__;
do {
kd += nTep;
id += nTep;
++nTep;
sum += b[id] * b[kd];
} while (nTep < n);
L40:
sum = a[kpiv] - sum;
L42:
if (kpiv < ipiv) b[kpiv] = sum * r__;
else {
dc = TMath::Sqrt(sum);
b[kpiv] = dc;
if (r__ > 0.) r__ = 1. / dc;
}
--kpiv;
} while (kpiv > ipiv - i__);
--i__;
} while (i__ > 0);
++b;
return b;
} /* trchul_ */
//____________________________________________________________
/* Subroutine */float *TCL::trinv(float *t, float *s, int n)
{
// trinv.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRINV .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
int lhor, ipiv, lver, j;
double sum = 0;
double r__ = 0;
int mx, ndTep, ind;
/* Parameter adjuTments */
--s; --t;
/* Function Body */
mx = (n * n + n) / 2;
ipiv = mx;
int i = n;
do {
r__ = 0.;
if (t[ipiv] > 0.) r__ = 1. / t[ipiv];
s[ipiv] = r__;
ndTep = n;
ind = mx - n + i;
while (ind != ipiv) {
sum = 0.;
if (r__ != 0.) {
lhor = ipiv;
lver = ind;
j = i;
do {
lhor += j;
++lver;
sum += t[lhor] * s[lver];
++j;
} while (lhor < ind);
}
s[ind] = -sum * r__;
--ndTep;
ind -= ndTep;
}
ipiv -= i;
--i;
} while (i > 0);
++s;
return s;
} /* trinv_ */
//____________________________________________________________
// trla.F -- translated by f2c (version 19970219).
/* Subroutine */float *TCL::trla(float *u, float *a, float *b, int m, int n)
{
int ipiv, ia, ib, iu;
double sum;
/* CERN PROGLIB# F112 TRLA .VERSION KERNFOR 4.15 861204 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--b; --a; --u;
/* Function Body */
ib = m * n;
ipiv = (m * m + m) / 2;
do {
do {
ia = ib;
iu = ipiv;
sum = 0.;
do {
sum += a[ia] * u[iu];
--iu;
ia -= n;
} while (ia > 0);
b[ib] = sum;
--ib;
} while (ia > 1 - n);
ipiv = iu;
} while (iu > 0);
++b;
return b;
} /* trla_ */
//____________________________________________________________
/* trlta.F -- translated by f2c (version 19970219).
// Subroutine */float *TCL::trlta(float *u, float *a, float *b, int m, int n)
{
int ipiv, mxpn, i__, nTep, ia, ib, iu, mx;
double sum;
/* CERN PROGLIB# F112 TRLTA .VERSION KERNFOR 4.15 861204 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--b; --a; --u;
/* Function Body */
ipiv = 0;
mx = m * n;
mxpn = mx + n;
ib = 0;
i__ = 0;
do {
++i__;
ipiv += i__;
do {
iu = ipiv;
nTep = i__;
++ib;
ia = ib;
sum = 0.;
do {
sum += a[ia] * u[iu];
ia += n;
iu += nTep;
++nTep;
} while (ia <= mx);
b[ib] = sum;
} while (ia < mxpn);
} while (i__ < m);
++b;
return b;
} /* trlta_ */
//____________________________________________________________
float *TCL::trpck(float *s, float *u, int n)
{
// trpck.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRPCK .VERSION KERNFOR 2.08 741218 */
// ORIG. 18/12/74 WH */
int i__, ia, ind, ipiv;
/* Parameter adjuTments */
--u; --s;
/* Function Body */
ia = 0;
ind = 0;
ipiv = 0;
for (i__ = 1; i__ <= n; ++i__) {
ipiv += i__;
do {
++ia;
++ind;
u[ind] = s[ia];
} while (ind < ipiv);
ia = ia + n - i__;
}
++u;
return u;
} /* trpck_ */
//____________________________________________________________
float *TCL::trqsq(float *q, float *s, float *r__, int m)
{
// trqsq.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRQSQ .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
int indq, inds, imax, i__, j, k, l;
int iq, ir, is, iqq;
double sum;
/* Parameter adjuTments */
--r__; --s; --q;
/* Function Body */
imax = (m * m + m) / 2;
vzero(&r__[1], imax);
inds = 0;
i__ = 0;
do {
inds += i__;
ir = 0;
indq = 0;
j = 0;
do {
indq += j;
is = inds;
iq = indq;
sum = (float)0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
if (k > j) iq += k;
else ++iq;
sum += s[is] * q[iq];
++k;
} while (k < m);
iqq = inds;
l = 0;
do {
++ir;
if (l > i__) iqq += l;
else ++iqq;
r__[ir] += q[iqq] * sum;
++l;
} while (l <= j);
++j;
} while (j < m);
++i__;
} while (i__ < m);
++r__;
return r__;
} /* trqsq_ */
//____________________________________________________________
float *TCL::trsa(float *s, float *a, float *b, int m, int n)
{
// trsa.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRSA .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
/* Local variables */
int inds, i__, j, k, ia, ib, is;
double sum;
/* Parameter adjuTments */
--b; --a; --s;
/* Function Body */
inds = 0;
ib = 0;
i__ = 0;
do {
inds += i__;
for (j = 1; j <= n; ++j) {
ia = j;
is = inds;
sum = 0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
sum += s[is] * a[ia];
ia += n;
++k;
} while (k < m);
++ib;
b[ib] = sum;
}
++i__;
} while (i__ < m);
++b;
return b;
} /* trsa_ */
//____________________________________________________________
/* Subroutine */float *TCL::trsinv(float *g, float *gi, int n)
{
// trsinv.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRSINV .VERSION KERNFOR 2.08 741218
// ORIG. 18/12/74 WH */
/* Function Body */
trchlu(g, gi, n);
trinv(gi, gi, n);
return trsmul(gi, gi, n);
} /* trsinv_ */
//____________________________________________________________
/* Subroutine */float *TCL::trsmlu(float *u, float *s, int n)
{
// trsmlu.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRSMLU .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
/* Local variables */
int lhor, lver, i__, k, l, ind;
double sum;
/* Parameter adjuTments */
--s; --u;
/* Function Body */
ind = (n * n + n) / 2;
for (i__ = 1; i__ <= n; ++i__) {
lver = ind;
for (k = i__; k <= n; ++k,--ind) {
lhor = ind; sum = 0.;
for (l = k; l <= n; ++l,--lver,--lhor)
sum += u[lver] * u[lhor];
s[ind] = sum;
}
}
++s;
return s;
} /* trsmlu_ */
//____________________________________________________________
/* Subroutine */float *TCL::trsmul(float *g, float *gi, int n)
{
// trsmul.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRSMUL .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
/* Local variables */
int lhor, lver, lpiv, i__, j, k, ind;
double sum;
/* Parameter adjuTments */
--gi; --g;
/* Function Body */
ind = 1;
lpiv = 0;
for (i__ = 1; i__ <= n; ++i__) {
lpiv += i__;
for (j = 1; j <= i__; ++j,++ind) {
lver = lpiv;
lhor = ind;
sum = 0.;
for (k = i__; k <= n; lhor += k,lver += k,++k)
sum += g[lver] * g[lhor];
gi[ind] = sum;
}
}
++gi;
return gi;
} /* trsmul_ */
//____________________________________________________________
float *TCL::trupck(float *u, float *s, int m)
{
// trupck.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRUPCK .VERSION KERNFOR 2.08 741218
// ORIG. 18/12/74 WH
int i__, im, is, iu, iv, ih, m2;
/* Parameter adjuTments */
--s; --u;
/* Function Body */
m2 = m * m;
is = m2;
iu = (m2 + m) / 2;
i__ = m - 1;
do {
im = i__ * m;
do {
s[is] = u[iu];
--is;
--iu;
} while (is > im);
is = is - m + i__;
--i__;
} while (i__ >= 0);
is = 1;
do {
iv = is;
ih = is;
while (1) {
iv += m;
++ih;
if (iv > m2) break;
s[ih] = s[iv];
}
is = is + m + 1;
} while (is < m2);
++s;
return s;
} /* trupck_ */
//____________________________________________________________
/* trsat.F -- translated by f2c (version 19970219).
// Subroutine */ float *TCL::trsat(float *s, float *a, float *b, int m, int n)
{
/* Local variables */
int inds, i__, j, k, ia, ib, mn, is;
double sum;
/* CERN PROGLIB# F112 TRSAT .VERSION KERNFOR 4.15 861204 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--b; --a; --s;
/* Function Body */
mn = m * n;
inds = 0;
ib = 0;
i__ = 0;
do {
inds += i__;
ia = 0;
for (j = 1; j <= n; ++j) {
is = inds;
sum = 0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
++ia;
sum += s[is] * a[ia];
++k;
} while (k < m);
++ib;
b[ib] = sum;
}
++i__;
} while (i__ < m);
++b;
return b;
} /* trsat_ */
// ------ double
//____________________________________________________________
// trata.F -- translated by f2c (version 19970219).
double *TCL::trata(double *a, double *r__, int m, int n)
{
/* Local variables */
int i__, j, ia, mn, ir, iat;
double sum;
/* CERN PROGLIB# F112 TRATA .VERSION KERNFOR 4.15 861204 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--r__; --a;
/* Function Body */
mn = m * n;
ir = 0;
for (i__ = 1; i__ <= m; ++i__) {
for (j = 1; j <= i__; ++j) {
ia = i__;
iat = j;
sum = (double)0.;
do {
sum += a[ia] * a[iat];
ia += m;
iat += m;
} while (ia <= mn);
++ir;
r__[ir] = sum;
}
}
return 0;
} /* trata_ */
//____________________________________________________________
// trats.F -- translated by f2c (version 19970219).
double *TCL::trats(double *a, double *s, double *b, int m, int n)
{
/* Local variables */
int inds, i__, j, k, ia, ib, is;
double sum;
/* CERN PROGLIB# F112 TRATS .VERSION KERNFOR 4.15 861204 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--b; --s; --a;
/* Function Body */
ib = 0; inds = 0; i__ = 0;
do {
inds += i__;
ib = i__ + 1;
for (j = 1; j <= m; ++j) {
ia = j;
is = inds;
sum = (double)0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
sum += a[ia] * s[is];
ia += m;
++k;
} while (k < n);
b[ib] = sum;
ib += n;
}
++i__;
} while (i__ < n);
return 0;
} /* trats_ */
//____________________________________________________________
// tratsa.F -- translated by f2c (version 19970219).
/* Subroutine */double *TCL::tratsa(double *a, double *s, double *r__, int m, int n)
{
/* Local variables */
int imax, i__, j, k;
int ia, mn, ir, is, iaa, ind;
double sum;
/* CERN PROGLIB# F112 TRATSA .VERSION KERNFOR 4.15 861204 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--r__; --s; --a;
/* Function Body */
imax = (m * m + m) / 2;
vzero(&r__[1], imax);
mn = m * n;
ind = 0;
i__ = 0;
do {
ind += i__;
ir = 0;
for (j = 1; j <= m; ++j) {
is = ind;
ia = j;
sum = (double)0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
sum += s[is] * a[ia];
ia += m;
++k;
} while (k < n);
iaa = i__ * m;
for (k = 1; k <= j; ++k) {
++iaa;
++ir;
r__[ir] += sum * a[iaa];
}
}
++i__;
} while (i__ < n);
return 0;
} /* tratsa_ */
//____________________________________________________________
double *TCL::trchlu(double *a, double *b, int n)
{
// trchlu.F -- translated by f2c (version 19970219).
/* Local variables */
int ipiv, kpiv, i__, j;
double r__, dc;
int id, kd;
double sum;
/* CERN PROGLIB# F112 TRCHLU .VERSION KERNFOR 4.16 870601 */
/* ORIG. 18/12/74 W.HART */
/* Parameter adjuTments */
--b; --a;
/* Function Body */
ipiv = 0;
i__ = 0;
do {
++i__;
ipiv += i__;
kpiv = ipiv;
r__ = a[ipiv];
for (j = i__; j <= n; ++j) {
sum = 0.;
if (i__ == 1) goto L40;
if (r__ == 0.) goto L42;
id = ipiv - i__ + 1;
kd = kpiv - i__ + 1;
do {
sum += b[kd] * b[id];
++kd; ++id;
} while (id < ipiv);
L40:
sum = a[kpiv] - sum;
L42:
if (j != i__) b[kpiv] = sum * r__;
else {
dc = TMath::Sqrt(sum);
b[kpiv] = dc;
if (r__ > 0.) r__ = (double)1. / dc;
}
kpiv += j;
}
} while (i__ < n);
return 0;
} /* trchlu_ */
//____________________________________________________________
// trchul.F -- translated by f2c (version 19970219).
double *TCL::trchul(double *a, double *b, int n)
{
/* Local variables */
int ipiv, kpiv, i__;
double r__;
int nTep;
double dc;
int id, kd;
double sum;
/* CERN PROGLIB# F112 TRCHUL .VERSION KERNFOR 4.16 870601 */
/* ORIG. 18/12/74 WH */
/* Parameter adjuTments */
--b; --a;
/* Function Body */
kpiv = (n * n + n) / 2;
i__ = n;
do {
ipiv = kpiv;
r__ = a[ipiv];
do {
sum = 0.;
if (i__ == n) goto L40;
if (r__ == (double)0.) goto L42;
id = ipiv;
kd = kpiv;
nTep = i__;
do {
kd += nTep;
id += nTep;
++nTep;
sum += b[id] * b[kd];
} while (nTep < n);
L40:
sum = a[kpiv] - sum;
L42:
if (kpiv < ipiv) b[kpiv] = sum * r__;
else {
dc = TMath::Sqrt(sum);
b[kpiv] = dc;
if (r__ > (double)0.) r__ = (double)1. / dc;
}
--kpiv;
} while (kpiv > ipiv - i__);
--i__;
} while (i__ > 0);
return 0;
} /* trchul_ */
//____________________________________________________________
double *TCL::trinv(double *t, double *s, int n)
{
// trinv.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRINV .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
//
int lhor, ipiv, lver, j;
double r__;
int mx, ndTep, ind;
double sum;
/* Parameter adjuTments */
--s; --t;
/* Function Body */
mx = (n * n + n) / 2;
ipiv = mx;
int i = n;
do {
r__ = 0.;
if (t[ipiv] > 0.) r__ = (double)1. / t[ipiv];
s[ipiv] = r__;
ndTep = n;
ind = mx - n + i;
while (ind != ipiv) {
sum = 0.;
if (r__ != 0.) {
lhor = ipiv;
lver = ind;
j = i;
do {
lhor += j;
++lver;
sum += t[lhor] * s[lver];
++j;
} while (lhor < ind);
}
s[ind] = -sum * r__;
--ndTep;
ind -= ndTep;
}
ipiv -= i;
--i;
} while (i > 0);
return 0;
} /* trinv_ */
//____________________________________________________________
/* Subroutine */double *TCL::trla(double *u, double *a, double *b, int m, int n)
{
//
// trla.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRLA .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
//
int ipiv, ia, ib, iu;
double sum;
/* Parameter adjuTments */
--b; --a; --u;
/* Function Body */
ib = m * n;
ipiv = (m * m + m) / 2;
do {
do {
ia = ib;
iu = ipiv;
sum = 0.;
do {
sum += a[ia] * u[iu];
--iu;
ia -= n;
} while (ia > 0);
b[ib] = sum;
--ib;
} while (ia > 1 - n);
ipiv = iu;
} while (iu > 0);
return 0;
} /* trla_ */
//____________________________________________________________
double *TCL::trlta(double *u, double *a, double *b, int m, int n)
{
// trlta.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRLTA .VERSION KERNFOR 4.15 861204
// ORIG. 18/12/74 WH
int ipiv, mxpn, i__, nTep, ia, ib, iu, mx;
double sum;
/* Parameter adjuTments */
--b; --a; --u;
/* Function Body */
ipiv = 0;
mx = m * n;
mxpn = mx + n;
ib = 0;
i__ = 0;
do {
++i__;
ipiv += i__;
do {
iu = ipiv;
nTep = i__;
++ib;
ia = ib;
sum = 0.;
do {
sum += a[ia] * u[iu];
ia += n;
iu += nTep;
++nTep;
} while (ia <= mx);
b[ib] = sum;
} while (ia < mxpn);
} while (i__ < m);
return 0;
} /* trlta_ */
//____________________________________________________________
/* Subroutine */double *TCL::trpck(double *s, double *u, int n)
{
// trpck.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRPCK .VERSION KERNFOR 2.08 741218 */
// ORIG. 18/12/74 WH */
int i__, ia, ind, ipiv;
/* Parameter adjuTments */
--u; --s;
/* Function Body */
ia = 0;
ind = 0;
ipiv = 0;
for (i__ = 1; i__ <= n; ++i__) {
ipiv += i__;
do {
++ia;
++ind;
u[ind] = s[ia];
} while (ind < ipiv);
ia = ia + n - i__;
}
return 0;
} /* trpck_ */
//____________________________________________________________
double *TCL::trqsq(double *q, double *s, double *r__, int m)
{
// trqsq.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRQSQ .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
int indq, inds, imax, i__, j, k, l;
int iq, ir, is, iqq;
double sum;
/* Parameter adjuTments */
--r__; --s; --q;
/* Function Body */
imax = (m * m + m) / 2;
vzero(&r__[1], imax);
inds = 0;
i__ = 0;
do {
inds += i__;
ir = 0;
indq = 0;
j = 0;
do {
indq += j;
is = inds;
iq = indq;
sum = 0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
if (k > j) iq += k;
else ++iq;
sum += s[is] * q[iq];
++k;
} while (k < m);
iqq = inds;
l = 0;
do {
++ir;
if (l > i__) iqq += l;
else ++iqq;
r__[ir] += q[iqq] * sum;
++l;
} while (l <= j);
++j;
} while (j < m);
++i__;
} while (i__ < m);
return 0;
} /* trqsq_ */
//____________________________________________________________
double *TCL::trsa(double *s, double *a, double *b, int m, int n)
{
// trsa.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRSA .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
/* Local variables */
int inds, i__, j, k, ia, ib, is;
double sum;
/* Parameter adjuTments */
--b; --a; --s;
/* Function Body */
inds = 0;
ib = 0;
i__ = 0;
do {
inds += i__;
for (j = 1; j <= n; ++j) {
ia = j;
is = inds;
sum = 0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
sum += s[is] * a[ia];
ia += n;
++k;
} while (k < m);
++ib;
b[ib] = sum;
}
++i__;
} while (i__ < m);
return 0;
} /* trsa_ */
//____________________________________________________________
/* Subroutine */double *TCL::trsinv(double *g, double *gi, int n)
{
// trsinv.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRSINV .VERSION KERNFOR 2.08 741218
// ORIG. 18/12/74 WH */
/* Function Body */
trchlu(g, gi, n);
trinv(gi, gi, n);
trsmul(gi, gi, n);
return 0;
} /* trsinv_ */
//____________________________________________________________
/* Subroutine */double *TCL::trsmlu(double *u, double *s, int n)
{
// trsmlu.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRSMLU .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
/* Local variables */
int lhor, lver, i__, k, l, ind;
double sum;
/* Parameter adjuTments */
--s; --u;
/* Function Body */
ind = (n * n + n) / 2;
for (i__ = 1; i__ <= n; ++i__) {
lver = ind;
for (k = i__; k <= n; ++k,--ind) {
lhor = ind; sum = 0.;
for (l = k; l <= n; ++l,--lver,--lhor)
sum += u[lver] * u[lhor];
s[ind] = sum;
}
}
return 0;
} /* trsmlu_ */
//____________________________________________________________
/* Subroutine */double *TCL::trsmul(double *g, double *gi, int n)
{
// trsmul.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRSMUL .VERSION KERNFOR 4.15 861204 */
// ORIG. 18/12/74 WH */
/* Local variables */
int lhor, lver, lpiv, i__, j, k, ind;
double sum;
/* Parameter adjuTments */
--gi; --g;
/* Function Body */
ind = 1;
lpiv = 0;
for (i__ = 1; i__ <= n; ++i__) {
lpiv += i__;
for (j = 1; j <= i__; ++j,++ind) {
lver = lpiv;
lhor = ind;
sum = 0.;
for (k = i__; k <= n;lhor += k,lver += k,++k)
sum += g[lver] * g[lhor];
gi[ind] = sum;
}
}
return 0;
} /* trsmul_ */
//____________________________________________________________
/* Subroutine */double *TCL::trupck(double *u, double *s, int m)
{
// trupck.F -- translated by f2c (version 19970219).
// CERN PROGLIB# F112 TRUPCK .VERSION KERNFOR 2.08 741218
// ORIG. 18/12/74 WH
int i__, im, is, iu, iv, ih, m2;
/* Parameter adjuTments */
--s; --u;
/* Function Body */
m2 = m * m;
is = m2;
iu = (m2 + m) / 2;
i__ = m - 1;
do {
im = i__ * m;
do {
s[is] = u[iu];
--is;
--iu;
} while (is > im);
is = is - m + i__;
--i__;
} while (i__ >= 0);
is = 1;
do {
iv = is;
ih = is;
while (1) {
iv += m;
++ih;
if (iv > m2) break;
s[ih] = s[iv];
}
is = is + m + 1;
} while (is < m2);
return 0;
} /* trupck_ */
//____________________________________________________________
double *TCL::trsat(double *s, double *a, double *b, int m, int n)
{
// trsat.F -- translated by f2c (version 19970219)
// CERN PROGLIB# F112 TRSAT .VERSION KERNFOR 4.15 861204
// ORIG. 18/12/74 WH
/* Local variables */
int inds, i__, j, k, ia, ib, mn, is;
double sum;
/* Parameter adjuTments */
--b; --a; --s;
/* Function Body */
mn = m * n;
inds = 0;
ib = 0;
i__ = 0;
do {
inds += i__;
ia = 0;
for (j = 1; j <= n; ++j) {
is = inds;
sum = 0.;
k = 0;
do {
if (k > i__) is += k;
else ++is;
++ia;
sum += s[is] * a[ia];
++k;
} while (k < m);
++ib;
b[ib] = sum;
}
++i__;
} while (i__ < m);
return 0;
} /* trsat_ */
// ------------ Victor Perevoztchikov's addition
//_____________________________________________________________________________
float *TCL::trsequ(float *smx, int m, float *b, int n)
{
// Linear Equations, Matrix Inversion
// trsequ solves the matrix equation
//
// SMX*x = B
//
// which represents a system of m simultaneous linear equations with n right-hand sides:
// SMX is an unpacked symmetric matrix (all elements) (m x m)
// B is an unpacked matrix of right-hand sides (n x m)
//
float *mem = new float[(m*(m+1))/2+m];
float *v = mem;
float *s = v+m;
if (!b) n=0;
TCL::trpck (smx,s ,m);
TCL::trsinv(s ,s, m);
for (int i=0;i<n;i++) {
TCL::trsa (s ,b+i*m, v, m, 1);
TCL::ucopy (v ,b+i*m, m);}
TCL::trupck(s ,smx, m);
delete [] mem;
return b;
}
//_____________________________________________________________________________
double *TCL::trsequ(double *smx, int m, double *b, int n)
{
// Linear Equations, Matrix Inversion
// trsequ solves the matrix equation
//
// SMX*x = B
//
// which represents a system of m simultaneous linear equations with n right-hand sides:
// SMX is an unpacked symmetric matrix (all elements) (m x m)
// B is an unpacked matrix of right-hand sides (n x m)
//
double *mem = new double[(m*(m+1))/2+m];
double *v = mem;
double *s = v+m;
if (!b) n=0;
TCL::trpck (smx,s ,m);
TCL::trsinv(s ,s, m);
for (int i=0;i<n;i++) {
TCL::trsa (s ,b+i*m, v, m, 1);
TCL::ucopy (v ,b+i*m, m);}
TCL::trupck(s ,smx, m);
delete [] mem;
return b;
}