Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

hholder.cpp

Go to the documentation of this file.
00001 //$$ hholder.cpp                                   QR decomposition
00002 
00003 // Copyright (C) 1991,2,3,4: R B Davies
00004 
00005 #ifndef WANT_MATH
00006 #define WANT_MATH
00007 #endif
00008 //#define WANT_STREAM
00009 
00010 #include "include.h"
00011 
00012 #include "newmatap.h"
00013 
00014 #ifdef use_namespace
00015 namespace NEWMAT {
00016 #endif
00017 
00018 #ifdef DO_REPORT
00019 #define REPORT { static ExeCounter ExeCount(__LINE__,16); ++ExeCount; }
00020 #else
00021 #define REPORT {}
00022 #endif
00023 
00024 
00025 /*************************** QR decompositions ***************************/
00026 
00027 inline Real square(Real x) { return x*x; }
00028 
00029 void QRZT(Matrix& X, LowerTriangularMatrix& L)
00030 {
00031    REPORT
00032    Tracer et("QRZT(1)");
00033    int n = X.Ncols(); int s = X.Nrows(); L.ReSize(s);
00034    if (n == 0 || s == 0) { L = 0.0; return; }
00035    Real* xi = X.Store(); int k;
00036    for (int i=0; i<s; i++)
00037    {
00038       Real sum = 0.0;
00039       Real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
00040       sum = sqrt(sum);
00041       if (sum == 0.0)
00042       {
00043          REPORT
00044          k=n; while(k--) { *xi0++ = 0.0; }
00045          for (int j=i; j<s; j++) L.element(j,i) = 0.0;
00046       }
00047       else
00048       {
00049          L.element(i,i) = sum;
00050          Real* xj0=xi0; k=n; while(k--) { *xj0++ /= sum; }
00051          for (int j=i+1; j<s; j++)
00052          {
00053             sum=0.0;
00054             xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
00055             xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
00056             L.element(j,i) = sum;
00057          }
00058       }
00059    }
00060 }
00061 
00062 void QRZT(const Matrix& X, Matrix& Y, Matrix& M)
00063 {
00064    REPORT
00065    Tracer et("QRZT(2)");
00066    int n = X.Ncols(); int s = X.Nrows(); int t = Y.Nrows();
00067    if (Y.Ncols() != n)
00068       { Throw(ProgramException("Unequal row lengths",X,Y)); }
00069    M.ReSize(t,s);
00070    Real* xi = X.Store(); int k;
00071    for (int i=0; i<s; i++)
00072    {
00073       Real* xj0 = Y.Store(); Real* xi0 = xi;
00074       for (int j=0; j<t; j++)
00075       {
00076          Real sum=0.0;
00077          xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
00078          xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
00079          M.element(j,i) = sum;
00080       }
00081    }
00082 }
00083 
00084 /*
00085 void QRZ(Matrix& X, UpperTriangularMatrix& U)
00086 {
00087   Tracer et("QRZ(1)");
00088   int n = X.Nrows(); int s = X.Ncols(); U.ReSize(s);
00089   Real* xi0 = X.Store(); int k;
00090   for (int i=0; i<s; i++)
00091   {
00092     Real sum = 0.0;
00093     Real* xi = xi0; k=n; while(k--) { sum += square(*xi); xi+=s; }
00094     sum = sqrt(sum);
00095     U.element(i,i) = sum;
00096     if (sum==0.0) Throw(SingularException(U));
00097     Real* xj0=xi0; k=n; while(k--) { *xj0 /= sum; xj0+=s; }
00098     xj0 = xi0;
00099     for (int j=i+1; j<s; j++)
00100     {
00101       sum=0.0;
00102       xi=xi0; k=n; xj0++; Real* xj=xj0;
00103       while(k--) { sum += *xi * *xj; xi+=s; xj+=s; }
00104       xi=xi0; k=n; xj=xj0;
00105       while(k--) { *xj -= sum * *xi; xj+=s; xi+=s; }
00106       U.element(i,j) = sum;
00107     }
00108     xi0++;
00109   }
00110 }
00111 */
00112 
00113 void QRZ(Matrix& X, UpperTriangularMatrix& U)
00114 {
00115    REPORT
00116    Tracer et("QRZ(1)");
00117    int n = X.Nrows(); int s = X.Ncols(); U.ReSize(s); U = 0.0;
00118    if (n == 0 || s == 0) return;
00119    Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
00120    int j, k; int J = s; int i = s;
00121    while (i--)
00122    {
00123       Real* xj0 = xi0; Real* xi = xi0; k = n;
00124       if (k) for (;;)
00125       {
00126          u = u0; Real Xi = *xi; Real* xj = xj0;
00127          j = J; while(j--) *u++ += Xi * *xj++;
00128          if (!(--k)) break;
00129          xi += s; xj0 += s;
00130       }
00131 
00132       Real sum = sqrt(*u0); *u0 = sum; u = u0+1;
00133       if (sum == 0.0)
00134       {
00135          REPORT
00136          j = J - 1; while(j--) *u++ = 0.0;
00137 
00138          xj0 = xi0++; k = n;
00139          if (k) for (;;)
00140          {
00141             *xj0 = 0.0;
00142             if (!(--k)) break;
00143             xj0 += s;
00144          }
00145          u0 += J--;
00146       }
00147       else
00148       {
00149          int J1 = J-1; j = J1; while(j--) *u++ /= sum;
00150 
00151          xj0 = xi0; xi = xi0++; k = n;
00152          if (k) for (;;)
00153          {
00154             u = u0+1; Real Xi = *xi; Real* xj = xj0;
00155             Xi /= sum; *xj++ = Xi;
00156             j = J1; while(j--) *xj++ -= *u++ * Xi;
00157             if (!(--k)) break;
00158             xi += s; xj0 += s;
00159          }
00160          u0 += J--;
00161       }
00162    }
00163 }
00164 
00165 void QRZ(const Matrix& X, Matrix& Y, Matrix& M)
00166 {
00167    REPORT
00168    Tracer et("QRZ(2)");
00169    int n = X.Nrows(); int s = X.Ncols(); int t = Y.Ncols();
00170    if (Y.Nrows() != n)
00171       { Throw(ProgramException("Unequal column lengths",X,Y)); }
00172    M.ReSize(s,t); M = 0;Real* m0 = M.Store(); Real* m;
00173    Real* xi0 = X.Store();
00174    int j, k; int i = s;
00175    while (i--)
00176    {
00177       Real* xj0 = Y.Store(); Real* xi = xi0; k = n;
00178       if (k) for (;;)
00179       {
00180          m = m0; Real Xi = *xi; Real* xj = xj0;
00181          j = t; while(j--) *m++ += Xi * *xj++;
00182          if (!(--k)) break;
00183          xi += s; xj0 += t;
00184       }
00185 
00186       xj0 = Y.Store(); xi = xi0++; k = n;
00187       if (k) for (;;)
00188       {
00189          m = m0; Real Xi = *xi; Real* xj = xj0;
00190          j = t; while(j--) *xj++ -= *m++ * Xi;
00191          if (!(--k)) break;
00192          xi += s; xj0 += t;
00193       }
00194       m0 += t;
00195    }
00196 }
00197 
00198 /*
00199 
00200 void QRZ(const Matrix& X, Matrix& Y, Matrix& M)
00201 {
00202   Tracer et("QRZ(2)");
00203   int n = X.Nrows(); int s = X.Ncols(); int t = Y.Ncols();
00204   if (Y.Nrows() != n)
00205   { Throw(ProgramException("Unequal column lengths",X,Y)); }
00206   M.ReSize(s,t);
00207   Real* xi0 = X.Store(); int k;
00208   for (int i=0; i<s; i++)
00209   {
00210     Real* xj0 = Y.Store();
00211     for (int j=0; j<t; j++)
00212     {
00213       Real sum=0.0;
00214       Real* xi=xi0; Real* xj=xj0; k=n;
00215       while(k--) { sum += *xi * *xj; xi+=s; xj+=t; }
00216       xi=xi0; k=n; xj=xj0++;
00217       while(k--) { *xj -= sum * *xi; xj+=t; xi+=s; }
00218       M.element(i,j) = sum;
00219     }
00220     xi0++;
00221   }
00222 }
00223 */
00224 
00225 void UpdateQRZT(Matrix& X, LowerTriangularMatrix& L)
00226 {
00227    REPORT
00228    Tracer et("UpdateQRZT");
00229    int n = X.Ncols(); int s = X.Nrows();
00230    if (s != L.Nrows())
00231       Throw(ProgramException("Incompatible dimensions",X,L)); 
00232    if (n == 0 || s == 0) return;
00233    Real* xi = X.Store(); int k;
00234    for (int i=0; i<s; i++)
00235    {
00236       Real r = L.element(i,i); 
00237       Real sum = 0.0;
00238       Real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
00239       sum = sqrt(sum + square(r));
00240       if (sum == 0.0)
00241       {
00242          REPORT
00243          k=n; while(k--) { *xi0++ = 0.0; }
00244          for (int j=i; j<s; j++) L.element(j,i) = 0.0;
00245       }
00246       else
00247       {
00248          Real frs = fabs(r) + sum;
00249          Real a0 = sqrt(frs / sum); Real alpha = a0 / frs;
00250          if (r <= 0) { REPORT L.element(i,i) = sum; alpha = -alpha; }
00251          else { REPORT L.element(i,i) = -sum; }
00252          Real* xj0=xi0; k=n; while(k--) { *xj0++ *= alpha; }
00253          for (int j=i+1; j<s; j++)
00254          {
00255             sum = 0.0;
00256             xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
00257             sum += a0 * L.element(j,i);
00258             xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
00259             L.element(j,i) -= sum * a0;
00260          }
00261       }
00262    }
00263 }
00264 
00265 void UpdateQRZ(Matrix& X, UpperTriangularMatrix& U)
00266 {
00267    REPORT
00268    Tracer et("UpdateQRZ");
00269    int n = X.Nrows(); int s = X.Ncols();
00270    if (s != U.Ncols())
00271       Throw(ProgramException("Incompatible dimensions",X,U));
00272    if (n == 0 || s == 0) return; 
00273    Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
00274    RowVector V(s); Real* v0 = V.Store(); Real* v; V = 0.0;
00275    int j, k; int J = s; int i = s;
00276    while (i--)
00277    {
00278       Real* xj0 = xi0; Real* xi = xi0; k = n;
00279       if (k) for (;;)
00280       {
00281          v = v0; Real Xi = *xi; Real* xj = xj0;
00282          j = J; while(j--) *v++ += Xi * *xj++;
00283          if (!(--k)) break;
00284          xi += s; xj0 += s;
00285       }
00286 
00287       Real r = *u0;
00288       Real sum = sqrt(*v0 + square(r));
00289       
00290       if (sum == 0.0)
00291       {
00292          REPORT
00293          u = u0; v = v0;
00294          j = J; while(j--) { *u++ = 0.0; *v++ = 0.0; }
00295          xj0 = xi0++; k = n;
00296          if (k) for (;;)
00297          {
00298             *xj0 = 0.0;
00299             if (!(--k)) break;
00300             xj0 += s;
00301          }
00302          u0 += J--;
00303       }
00304       else
00305       {
00306          Real frs = fabs(r) + sum;
00307          Real a0 = sqrt(frs / sum); Real alpha = a0 / frs;
00308          if (r <= 0) { REPORT alpha = -alpha; *u0 = sum; }
00309          else { REPORT *u0 = -sum; }
00310       
00311          j = J - 1; v = v0 + 1; u = u0 + 1;     
00312          while (j--)
00313             { *v = a0 * *u + alpha * *v; *u -= a0 * *v; ++v; ++u; }
00314 
00315          xj0 = xi0; xi = xi0++; k = n;
00316          if (k) for (;;)
00317          {
00318             v = v0 + 1; Real Xi = *xi; Real* xj = xj0;
00319             Xi *= alpha; *xj++ = Xi;
00320             j = J - 1; while(j--) *xj++ -= *v++ * Xi;
00321             if (!(--k)) break;
00322             xi += s; xj0 += s;
00323          }
00324          
00325          j = J; v = v0;
00326          while (j--) *v++ = 0.0;
00327          
00328          u0 += J--;
00329       }
00330    }
00331 }
00332 
00333 
00334 #ifdef use_namespace
00335 }
00336 #endif
00337 

newmat11b
Generated Mon May 9 04:54:17 2016 by Doxygen 1.6.3