// @(#)root/star:$Name:  $:$Id: TCL.cxx,v 1.3 2003/01/27 20:41:36 brun 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.3 2003/01/27 20:41:36 brun Exp $
// $Log: TCL.cxx,v $
// Revision 1.3  2003/01/27 20:41:36  brun
// New version of the Table package by Valeri Fine.
// New classes TIndexTable TResponseIterator TResponseTable TTableMap
//
// Revision 1.3  2003/01/03 20:17:13  fisyak
// *** empty log message ***
//
// Revision 1.2  2002/11/11 12:02:40  brun

// Remove unused variable mn in several functions.

//
// Revision 1.1  2002/04/15 20:23:39  fine
// NEw naming schema for RootKErnel classes and a set of classes to back geometry OO
//
// 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, 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); 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, is; double sum; /* CERN PROGLIB# F112 TRSAT .VERSION KERNFOR 4.15 861204 */ /* ORIG. 18/12/74 WH */ /* Parameter adjuTments */ --b; --a; --s; /* Function Body */ 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, 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); 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, is; double sum; /* Parameter adjuTments */ --b; --a; --s; /* Function Body */ 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; }


ROOT page - Class index - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.