home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume21 / newmat02 / part03 / newmat8.cxx < prev    next >
C/C++ Source or Header  |  1991-08-01  |  7KB  |  310 lines
  1. //$$ newmat8.cxx         Advanced matrix functions, scalar functions
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #include "include.hxx"
  6.  
  7. #include "boolean.hxx"
  8.  
  9. typedef double real;                 // all floating point double
  10.  
  11. #include "newmat.hxx"
  12.  
  13. #include "newmatap.hxx"
  14.  
  15. #include "newmatrc.hxx"
  16.  
  17. //#define REPORT { static ExeCounter ExeCount(__LINE__,8); ExeCount++; }
  18.  
  19. #define REPORT {}
  20.  
  21.  
  22. /********************** householder triangularisation *********************/
  23.  
  24. inline real square(real x) { return x*x; }
  25.  
  26. void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
  27. {
  28.    int n = X.Ncols(); int s = X.Nrows(); L.ReDimension(s);
  29.    real* xi = X.Store(); int k;
  30.    for (int i=0; i<s; i++)
  31.    {
  32.       real sum = 0.0;
  33.       real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
  34.       sum = sqrt(sum);
  35.       L.element(i,i) = sum;
  36.       real* xj0=xi0; k=n; while(k--) { *xj0++ /= sum; }
  37.       for (int j=i+1; j<s; j++)
  38.       {
  39.          sum=0.0;
  40.      xi=xi0; real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
  41.      xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
  42.      L.element(j,i) = sum;
  43.       }
  44.    }
  45. }
  46.  
  47. void HHDecompose(const Matrix& X, Matrix& Y, Matrix& M)
  48. {
  49.    int n = X.Ncols(); int s = X.Nrows(); int t = Y.Nrows();
  50.    if (Y.Ncols() != n) MatrixError("Incompatible dimensions in HHDecompose");
  51.    M.ReDimension(t,s);
  52.    real* xi = X.Store(); int k;
  53.    for (int i=0; i<s; i++)
  54.    {
  55.       real* xj0 = Y.Store(); real* xi0 = xi;
  56.       for (int j=0; j<t; j++)
  57.       {
  58.          real sum=0.0;
  59.          xi=xi0; real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
  60.      xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
  61.      M.element(j,i) = sum;
  62.       }
  63.    }
  64. }
  65.  
  66. /********* Cholesky decomposition of a positive definite matrix *************/
  67.  
  68. // Suppose S is symmetrix and positive definite. Then there exists a unique
  69. // lower triangular matrix L such that L L.t() = S;
  70.  
  71. LowerTriangularMatrix Cholesky(const SymmetricMatrix& S)
  72. {
  73.    int nr = S.Nrows();
  74.    LowerTriangularMatrix T(nr);
  75.    real* s = S.Store(); real* t = T.Store(); real* tk = t;
  76.  
  77.    for (int i=0; i<nr; i++)
  78.    {
  79.       real* tj = t; real* ti = tk;
  80.       for (int j=0; j<=i; j++)
  81.       {
  82.      tk = ti; real sum = 0.0; int k = j;
  83.      while (k--) { sum += *tj++ * *tk++; }
  84.      if (i==j) *tk++ = sqrt(*s++ - sum);
  85.          else *tk++ = (*s++ - sum) / *tj++;
  86.       }
  87.    }
  88. #ifndef __ZTC__
  89.    T.Release();              // not wanted if routine avoids return init
  90. #endif
  91.    return T;
  92. }
  93.  
  94. /************************** LU transformation ****************************/
  95.  
  96. inline real fabs(real f) { return (f>0.0) ? f : -f; }
  97.  
  98. void CroutMatrix::ludcmp()
  99. // LU decomposition - from numerical recipes in C
  100. {
  101.    REPORT
  102.    int i,j;
  103.  
  104.    real* vv=new real [nrows]; MatrixErrorNoSpace(vv);
  105.    real* a;
  106.  
  107.    a=store;
  108.    for (i=0;i<nrows;i++)
  109.    {
  110.       real big=0.0;
  111.       j=nrows; while (j--) { real temp=fabs(*a++); if (temp > big) big=temp; }
  112.       if (big == 0.0) MatrixError("Singular matrix");
  113.       vv[i]=1.0/big;
  114.    }
  115.  
  116.    real* aj=store;
  117.    for (j=0;j<nrows;j++)
  118.    {
  119.       real* ai=store;
  120.       for (i=0;i<j;i++)
  121.       {
  122.          real sum=*(ai+j); real* aix=ai; real* ajx=aj;
  123.          int k=i; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  124.          *ajx = sum; ai += nrows;
  125.       }
  126.  
  127.       real big=0.0; int imax;
  128.       for (i=j;i<nrows;i++)
  129.       {
  130.          real sum=*(ai+j); real* aix=ai; real* ajx=aj;
  131.          int k=j; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  132.          *aix = sum; ai += nrows;
  133.          real dum=vv[i]*fabs(sum); if (dum >= big) { big=dum; imax=i; }
  134.       }
  135.  
  136.       if (j != imax)
  137.       {
  138.          real* amax=store+imax*nrows; real* ajx=store+j*nrows;
  139.          int k=nrows; while(k--) { real dum=*amax; *amax++=*ajx; *ajx++=dum; }
  140.          d=!d; vv[imax]=vv[j];
  141.       }
  142.  
  143.       indx[j]=imax; ai=aj+j*nrows;
  144.       if (*ai == 0.0) MatrixError("Singular matrix");
  145.       real dum=1.0/(*ai);
  146.       i=nrows-j; while (--i) { ai += nrows; *ai *= dum; }
  147.  
  148.       aj++;
  149.    }
  150.    delete [nrows] vv;
  151. }
  152.  
  153. void CroutMatrix::lubksb(real* B, int mini)
  154. {
  155.    REPORT
  156.    int i,j; int ii=-1; real* ai=store;
  157.  
  158.    for (i=0;i<nrows;i++)
  159.    {
  160.       int ip=indx[i]; real sum=B[ip]; B[ip]=B[i];
  161.       if (ii>=0)
  162.       {
  163.          real* aij=ai+ii; real* bj=B+ii; j=i-ii;
  164.          while (j--) { sum -= *aij++ * *bj; bj++; }
  165.       }
  166.       else if (sum) ii=i;
  167.       B[i]=sum; ai += nrows;
  168.    }
  169.  
  170.    for (i=nrows-1;i>=mini;i--)
  171.    {
  172.       real* bj=B+i; ai -= nrows; real* ajx=ai+i; real sum=*bj; bj++;
  173.       j=nrows-i; while(--j) { sum -= *(++ajx) * *bj; bj++; }
  174.       B[i]=sum/(*(ai+i));
  175.    }
  176. }
  177.  
  178. /****************************** scalar functions ****************************/
  179.  
  180. real GeneralMatrix::SumSquare()
  181. {
  182.    REPORT
  183.    real sum = 0.0; int i = storage; real* s = store;
  184.    while (i--) sum += square(*s++);
  185.    tDelete(); return sum;
  186. }
  187.  
  188. real SymmetricMatrix::SumSquare()
  189. {
  190.    REPORT
  191.    real sum1 = 0.0; real sum2 = 0.0; real* s = store; int nr = nrows;
  192.    for (int i = 0; i<nr; i++)
  193.    {
  194.       int j = i;
  195.       while (j--) sum2 += square(*s++);
  196.       sum1 += square(*s++);
  197.    }
  198.    tDelete(); return sum1 + 2.0 * sum2;
  199. }
  200.  
  201. real BaseMatrix::SumSquare()
  202. {
  203.    REPORT GeneralMatrix* gm = Evaluate();
  204.    real s = gm->SumSquare(); return s;
  205. }
  206.  
  207. real Matrix::Trace()
  208. {
  209.    REPORT
  210.    int i = nrows; int d = i+1;
  211.    if (i != ncols) MatrixError("trace of non-square matrix");
  212.    real sum = 0.0; real* s = store;
  213.    while (i--) { sum += *s; s += d; }
  214.    tDelete(); return sum;
  215. }
  216.  
  217. real DiagonalMatrix::Trace()
  218. {
  219.    REPORT
  220.    int i = nrows; real sum = 0.0; real* s = store;
  221.    while (i--) sum += *s++;
  222.    tDelete(); return sum;
  223. }
  224.  
  225. real SymmetricMatrix::Trace()
  226. {
  227.    REPORT
  228.    int i = nrows; real sum = 0.0; real* s = store; int j = 2;
  229.    while (i--) { sum += *s; s += j++; }
  230.    tDelete(); return sum;
  231. }
  232.  
  233. real LowerTriangularMatrix::Trace()
  234. {
  235.    REPORT
  236.    int i = nrows; real sum = 0.0; real* s = store; int j = 2;
  237.    while (i--) { sum += *s; s += j++; }
  238.    tDelete(); return sum;
  239. }
  240.  
  241. real UpperTriangularMatrix::Trace()
  242. {
  243.    REPORT
  244.    int i = nrows; real sum = 0.0; real* s = store;
  245.    while (i) { sum += *s; s += i--; }
  246.    tDelete(); return sum;
  247. }
  248.  
  249. real BaseMatrix::Trace()
  250. {
  251.    REPORT
  252.    GeneralMatrix* gm = Evaluate(MatrixType::Diag);
  253.    real sum = gm->Trace(); return sum;
  254. }
  255.  
  256. void LogAndSign::operator*=(real x)
  257. {
  258.    if (x > 0.0) { log_value += log(x); }
  259.    else if (x < 0.0) { log_value += log(-x); sign = -sign; }
  260.    else sign = 0;
  261. }
  262.  
  263. real LogAndSign::Value() { return sign * exp(log_value); }
  264.  
  265. LogAndSign DiagonalMatrix::LogDeterminant()
  266. {
  267.    REPORT
  268.    int i = nrows; LogAndSign sum; real* s = store;
  269.    while (i--) sum *= *s++;
  270.    tDelete(); return sum;
  271. }
  272.  
  273. LogAndSign LowerTriangularMatrix::LogDeterminant()
  274. {
  275.    REPORT
  276.    int i = nrows; LogAndSign sum; real* s = store; int j = 2;
  277.    while (i--) { sum *= *s; s += j++; }
  278.    tDelete(); return sum;
  279. }
  280.  
  281. LogAndSign UpperTriangularMatrix::LogDeterminant()
  282. {
  283.    REPORT
  284.    int i = nrows; LogAndSign sum; real* s = store;
  285.    while (i) { sum *= *s; s += i--; }
  286.    tDelete(); return sum;
  287. }
  288.  
  289. LogAndSign BaseMatrix::LogDeterminant()
  290. {
  291.    REPORT GeneralMatrix* gm = Evaluate();
  292.    LogAndSign sum = gm->LogDeterminant(); return sum;
  293. }
  294.  
  295. LogAndSign GeneralMatrix::LogDeterminant()
  296. {
  297.    REPORT
  298.    if (nrows != ncols) MatrixError("determinant of non-square matrix"); 
  299.    CroutMatrix C(*this); return C.LogDeterminant();
  300. }
  301.  
  302. LogAndSign CroutMatrix::LogDeterminant()
  303. {
  304.    REPORT
  305.    int i = nrows; int dd = i+1; LogAndSign sum; real* s = store;
  306.    while (i--) { sum *= *s; s += dd; }
  307.    if (!d) sum.ChangeSign(); return sum;
  308.  
  309. }
  310.