home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / newmat03 / hholder.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-08  |  1.5 KB  |  56 lines
  1. //$$ hholder.cxx                       Householder triangularisation
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.hxx"
  8.  
  9. #include "newmatap.hxx"
  10.  
  11.  
  12. /********************** householder triangularisation *********************/
  13.  
  14. inline real square(real x) { return x*x; }
  15.  
  16. void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
  17. {
  18.    int n = X.Ncols(); int s = X.Nrows(); L.ReDimension(s);
  19.    real* xi = X.Store(); int k;
  20.    for (int i=0; i<s; i++)
  21.    {
  22.       real sum = 0.0;
  23.       real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
  24.       sum = sqrt(sum);
  25.       L.element(i,i) = sum;
  26.       real* xj0=xi0; k=n; while(k--) { *xj0++ /= sum; }
  27.       for (int j=i+1; j<s; j++)
  28.       {
  29.          sum=0.0;
  30.      xi=xi0; real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
  31.      xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
  32.      L.element(j,i) = sum;
  33.       }
  34.    }
  35. }
  36.  
  37. void HHDecompose(const Matrix& X, Matrix& Y, Matrix& M)
  38. {
  39.    int n = X.Ncols(); int s = X.Nrows(); int t = Y.Nrows();
  40.    if (Y.Ncols() != n) MatrixError("Incompatible dimensions in HHDecompose");
  41.    M.ReDimension(t,s);
  42.    real* xi = X.Store(); int k;
  43.    for (int i=0; i<s; i++)
  44.    {
  45.       real* xj0 = Y.Store(); real* xi0 = xi;
  46.       for (int j=0; j<t; j++)
  47.       {
  48.          real sum=0.0;
  49.          xi=xi0; real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
  50.      xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
  51.      M.element(j,i) = sum;
  52.       }
  53.    }
  54. }
  55.  
  56.