home *** CD-ROM | disk | FTP | other *** search
/ PC Media 23 / PC MEDIA CD23.iso / share / prog / newmat / hholder.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1995-01-11  |  4.2 KB  |  168 lines
  1. //$$ hholder.cpp                                   QR decomposition
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.h"
  8.  
  9. #include "newmatap.h"
  10.  
  11.  
  12. /*************************** QR decompositions ***************************/
  13.  
  14. inline Real square(Real x) { return x*x; }
  15.  
  16. void QRZT(Matrix& X, LowerTriangularMatrix& L)
  17. {
  18.     Tracer et("QZT(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 QRZT(const Matrix& X, Matrix& Y, Matrix& M)
  40. {
  41.     Tracer et("QRZT(2)");
  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. /*
  61. void QRZ(Matrix& X, UpperTriangularMatrix& U)
  62. {
  63.     Tracer et("QRZ(1)");
  64.     int n = X.Nrows(); int s = X.Ncols(); U.ReDimension(s);
  65.     Real* xi0 = X.Store(); int k;
  66.     for (int i=0; i<s; i++)
  67.     {
  68.         Real sum = 0.0;
  69.         Real* xi = xi0; k=n; while(k--) { sum += square(*xi); xi+=s; }
  70.         sum = sqrt(sum);
  71.         U.element(i,i) = sum;
  72.         if (sum==0.0) Throw(SingularException(U));
  73.         Real* xj0=xi0; k=n; while(k--) { *xj0 /= sum; xj0+=s; }
  74.         xj0 = xi0;
  75.         for (int j=i+1; j<s; j++)
  76.         {
  77.             sum=0.0;
  78.             xi=xi0; k=n; xj0++; Real* xj=xj0;
  79.             while(k--) { sum += *xi * *xj; xi+=s; xj+=s; }
  80.             xi=xi0; k=n; xj=xj0;
  81.             while(k--) { *xj -= sum * *xi; xj+=s; xi+=s; }
  82.             U.element(i,j) = sum;
  83.         }
  84.         xi0++;
  85.     }
  86. }
  87. */
  88.  
  89. void QRZ(Matrix& X, UpperTriangularMatrix& U)
  90. {
  91.     Tracer et("QRZ(1)");
  92.     int n = X.Nrows(); int s = X.Ncols(); U.ReDimension(s); U = 0.0;
  93.     Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
  94.     int j, k; int J = s; int i = s; while (i--)
  95.     {
  96.         Real* xj0 = xi0; Real* xi = xi0; k = n; while(k--)
  97.         {
  98.             u = u0; Real Xi = *xi; xi += s; Real* xj = xj0; xj0 += s;
  99.             j = J; while(j--) *u++ += Xi * *xj++;
  100.         }
  101.  
  102.         Real sum = sqrt(*u0); *u0 = sum; u = u0+1;
  103.         if (sum == 0.0) Throw(SingularException(U));
  104.         int J1 = J-1; j = J1; while(j--) *u++ /= sum;
  105.  
  106.         xj0 = xi0; xi = xi0++; k = n; while(k--)
  107.         {
  108.             u = u0+1; Real Xi = *xi; xi += s; Real* xj = xj0; xj0 += s;
  109.             Xi /= sum; *xj++ = Xi;
  110.             j = J1; while(j--) *xj++ -= *u++ * Xi;
  111.         }
  112.         u0 += J--;
  113.     }
  114. }
  115.  
  116. void QRZ(const Matrix& X, Matrix& Y, Matrix& M)
  117. {
  118.     Tracer et("QRZ(2)");
  119.     int n = X.Nrows(); int s = X.Ncols(); int t = Y.Ncols();
  120.     if (Y.Nrows() != n)
  121.     { Throw(ProgramException("Unequal column lengths",X,Y)); }
  122.     M.ReDimension(s,t); M = 0;Real* m0 = M.Store(); Real* m;
  123.     Real* xi0 = X.Store();
  124.     int j, k; int i = s;
  125.     while (i--)
  126.     {
  127.         Real* xj0 = Y.Store(); Real* xi = xi0; k = n; while(k--)
  128.         {
  129.             m = m0; Real Xi = *xi; xi += s; Real* xj = xj0; xj0 += t;
  130.             j = t; while(j--) *m++ += Xi * *xj++;
  131.         }
  132.  
  133.         xj0 = Y.Store(); xi = xi0++; k = n; while(k--)
  134.         {
  135.             m = m0; Real Xi = *xi; xi += s; Real* xj = xj0; xj0 += t;
  136.             j = t; while(j--) *xj++ -= *m++ * Xi;
  137.         }
  138.         m0 += t;
  139.     }
  140. }
  141.  
  142. /*
  143.  
  144. void QRZ(const Matrix& X, Matrix& Y, Matrix& M)
  145. {
  146.     Tracer et("QRZ(2)");
  147.     int n = X.Nrows(); int s = X.Ncols(); int t = Y.Ncols();
  148.     if (Y.Nrows() != n)
  149.     { Throw(ProgramException("Unequal column lengths",X,Y)); }
  150.     M.ReDimension(s,t);
  151.     Real* xi0 = X.Store(); int k;
  152.     for (int i=0; i<s; i++)
  153.     {
  154.         Real* xj0 = Y.Store();
  155.         for (int j=0; j<t; j++)
  156.         {
  157.             Real sum=0.0;
  158.             Real* xi=xi0; Real* xj=xj0; k=n;
  159.             while(k--) { sum += *xi * *xj; xi+=s; xj+=t; }
  160.             xi=xi0; k=n; xj=xj0++;
  161.             while(k--) { *xj -= sum * *xi; xj+=t; xi+=s; }
  162.             M.element(i,j) = sum;
  163.         }
  164.         xi0++;
  165.     }
  166. }
  167. */
  168.