HHOLDER.CPP Source File Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members
HHOLDER.CPP

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$ hholder.cpp                                   QR decomposition
00003 
00004 // Copyright (C) 1991,2,3,4: R B Davies
00005 
00006 #define WANT_MATH
00007 
00008 #include "include.h"
00009 
00010 #include "newmatap.h"
00011 
00012 #ifdef use_namespace
00013 namespace NEWMAT {
00014 #endif
00015 
00016     /*************************** QR decompositions ***************************/
00017 
00018     inline Real square(Real x) { return x*x; }
00019 
00020     void QRZT(Matrix& X, LowerTriangularMatrix& L) {
00021         Tracer et("QZT(1)");
00022         int n = X.Ncols(); int s = X.Nrows(); L.ReSize(s);
00023         Real* xi = X.Store(); int k;
00024         for (int i=0; i<s; i++) {
00025             Real sum = 0.0;
00026             Real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
00027             sum = sqrt(sum);
00028             L.element(i,i) = sum;
00029             if (sum==0.0) Throw(SingularException(L));
00030             Real* xj0=xi0; k=n; while(k--) { *xj0++ /= sum; }
00031             for (int j=i+1; j<s; j++) {
00032                 sum=0.0;
00033                 xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
00034                 xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
00035                 L.element(j,i) = sum;
00036             }
00037         }
00038     }
00039 
00040     void QRZT(const Matrix& X, Matrix& Y, Matrix& M) {
00041         Tracer et("QRZT(2)");
00042         int n = X.Ncols(); int s = X.Nrows(); int t = Y.Nrows();
00043         if (Y.Ncols() != n)
00044         { Throw(ProgramException("Unequal row lengths",X,Y)); }
00045         M.ReSize(t,s);
00046         Real* xi = X.Store(); int k;
00047         for (int i=0; i<s; i++) {
00048             Real* xj0 = Y.Store(); Real* xi0 = xi;
00049             for (int j=0; j<t; j++) {
00050                 Real sum=0.0;
00051                 xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
00052                 xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
00053                 M.element(j,i) = sum;
00054             }
00055         }
00056     }
00057 
00058     /*
00059       void QRZ(Matrix& X, UpperTriangularMatrix& U)
00060       {
00061       Tracer et("QRZ(1)");
00062       int n = X.Nrows(); int s = X.Ncols(); U.ReSize(s);
00063       Real* xi0 = X.Store(); int k;
00064       for (int i=0; i<s; i++)
00065       {
00066       Real sum = 0.0;
00067       Real* xi = xi0; k=n; while(k--) { sum += square(*xi); xi+=s; }
00068       sum = sqrt(sum);
00069       U.element(i,i) = sum;
00070       if (sum==0.0) Throw(SingularException(U));
00071       Real* xj0=xi0; k=n; while(k--) { *xj0 /= sum; xj0+=s; }
00072       xj0 = xi0;
00073       for (int j=i+1; j<s; j++)
00074       {
00075       sum=0.0;
00076       xi=xi0; k=n; xj0++; Real* xj=xj0;
00077       while(k--) { sum += *xi * *xj; xi+=s; xj+=s; }
00078       xi=xi0; k=n; xj=xj0;
00079       while(k--) { *xj -= sum * *xi; xj+=s; xi+=s; }
00080       U.element(i,j) = sum;
00081       }
00082       xi0++;
00083       }
00084       }
00085     */
00086 
00087     void QRZ(Matrix& X, UpperTriangularMatrix& U) {
00088         Tracer et("QRZ(1)");
00089         int n = X.Nrows(); int s = X.Ncols(); U.ReSize(s); U = 0.0;
00090         Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
00091         int j, k; int J = s; int i = s;
00092         while (i--) {
00093             Real* xj0 = xi0; Real* xi = xi0; k = n;
00094             while(k--) {
00095                 u = u0; Real Xi = *xi; xi += s; Real* xj = xj0; xj0 += s;
00096                 j = J; while(j--) *u++ += Xi * *xj++;
00097             }
00098 
00099             Real sum = sqrt(*u0); *u0 = sum; u = u0+1;
00100             if (sum == 0.0) Throw(SingularException(U));
00101             int J1 = J-1; j = J1; while(j--) *u++ /= sum;
00102 
00103             xj0 = xi0; xi = xi0++; k = n; while(k--) {
00104                 u = u0+1; Real Xi = *xi; xi += s; Real* xj = xj0; xj0 += s;
00105                 Xi /= sum; *xj++ = Xi;
00106                 j = J1; while(j--) *xj++ -= *u++ * Xi;
00107 
00108             }
00109             u0 += J--;
00110         }
00111     }
00112 
00113     void QRZ(const Matrix& X, Matrix& Y, Matrix& M) {
00114         Tracer et("QRZ(2)");
00115         int n = X.Nrows(); int s = X.Ncols(); int t = Y.Ncols();
00116         if (Y.Nrows() != n)
00117         { Throw(ProgramException("Unequal column lengths",X,Y)); }
00118         M.ReSize(s,t); M = 0;Real* m0 = M.Store(); Real* m;
00119         Real* xi0 = X.Store();
00120         int j, k; int i = s;
00121         while (i--) {
00122             Real* xj0 = Y.Store(); Real* xi = xi0; k = n; while(k--) {
00123                 m = m0; Real Xi = *xi; xi += s; Real* xj = xj0; xj0 += t;
00124                 j = t; while(j--) *m++ += Xi * *xj++;
00125             }
00126 
00127             xj0 = Y.Store(); xi = xi0++; k = n; while(k--) {
00128                 m = m0; Real Xi = *xi; xi += s; Real* xj = xj0; xj0 += t;
00129                 j = t; while(j--) *xj++ -= *m++ * Xi;
00130             }
00131             m0 += t;
00132         }
00133     }
00134 
00135     /*
00136 
00137       void QRZ(const Matrix& X, Matrix& Y, Matrix& M)
00138       {
00139       Tracer et("QRZ(2)");
00140       int n = X.Nrows(); int s = X.Ncols(); int t = Y.Ncols();
00141       if (Y.Nrows() != n)
00142       { Throw(ProgramException("Unequal column lengths",X,Y)); }
00143       M.ReSize(s,t);
00144       Real* xi0 = X.Store(); int k;
00145       for (int i=0; i<s; i++)
00146       {
00147       Real* xj0 = Y.Store();
00148       for (int j=0; j<t; j++)
00149       {
00150       Real sum=0.0;
00151       Real* xi=xi0; Real* xj=xj0; k=n;
00152       while(k--) { sum += *xi * *xj; xi+=s; xj+=t; }
00153       xi=xi0; k=n; xj=xj0++;
00154       while(k--) { *xj -= sum * *xi; xj+=t; xi+=s; }
00155       M.element(i,j) = sum;
00156       }
00157       xi0++;
00158       }
00159       }
00160     */
00161 
00162 #ifdef use_namespace
00163 }
00164 #endif
Generated on Fri Aug 23 01:37:02 2002 for VirtualU by
1.2.17