home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / n / newmat06.zip / HHOLDER.CXX < prev    next >
C/C++ Source or Header  |  1992-11-30  |  2KB  |  60 lines
  1. //$$ hholder.cxx                       Householder triangularisation
  2.  
  3. // Copyright (C) 1991,2: R B Davies
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.h"
  8.  
  9. #include "newmatap.h"
  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.    Tracer et("HHDecompose(1)");
  19.    int n = X.Ncols(); int s = X.Nrows(); L.ReDimension(s);
  20.    Real* xi = X.Store(); int k;
  21.    for (int i=0; i<s; i++)
  22.    {
  23.       Real sum = 0.0;
  24.       Real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
  25.       sum = sqrt(sum);
  26.       L.element(i,i) = sum;
  27.       if (sum==0.0) Throw(SingularException(L));
  28.       Real* xj0=xi0; k=n; while(k--) { *xj0++ /= sum; }
  29.       for (int j=i+1; j<s; j++)
  30.       {
  31.          sum=0.0;
  32.      xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
  33.      xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
  34.      L.element(j,i) = sum;
  35.       }
  36.    }
  37. }
  38.  
  39. void HHDecompose(const Matrix& X, Matrix& Y, Matrix& M)
  40. {
  41.    Tracer et("HHDecompose(1)");
  42.    int n = X.Ncols(); int s = X.Nrows(); int t = Y.Nrows();
  43.    if (Y.Ncols() != n)
  44.    { Throw(ProgramException("Unequal row lengths",X,Y)); }
  45.    M.ReDimension(t,s);
  46.    Real* xi = X.Store(); int k;
  47.    for (int i=0; i<s; i++)
  48.    {
  49.       Real* xj0 = Y.Store(); Real* xi0 = xi;
  50.       for (int j=0; j<t; j++)
  51.       {
  52.          Real sum=0.0;
  53.          xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
  54.      xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
  55.      M.element(j,i) = sum;
  56.       }
  57.    }
  58. }
  59.  
  60.