hholder.cpp
Go to the documentation of this file.00001
00002
00003
00004
00005 #ifndef WANT_MATH
00006 #define WANT_MATH
00007 #endif
00008
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
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
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
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
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
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